It has been 500 days since the last update, the content of the article may be outdated.
前面两节的模拟都是干空气泡的情形，本节介绍对Bryan湿空气团的模拟，会涉及到相变和热力学的能量交换过程。Bryan & Fritsch, 2002 提出，用来作为湿空气团的模拟标准。在情景上和Robert问题一样，不过是湿空气。
以下内容来百度百科
可逆干绝热过程
可逆湿绝热过程凝结并释放出潜热 ，如果饱和气块继续上升且凝结物全部保留 在气块内，并与外界无热量交换 ；当气块下沉增温湿，这些凝结物又蒸发 ，使气块始终维持饱和状态，所耗的潜热与原来释放的潜热相等，沿逆过程后仍能回到原来的状态，这样的过程称为可逆湿绝热过程 ，又称为湿绝热过程 ，这里的“湿”表示在饱和的绝热气块内发生水相变化。
可逆的绝热过程是等熵过程。等熵过程的对立面是等温过程，在等温过程中，最大限度的热量被转移到了外界，使得系统温度恒定如常。由于在热力学中，温度与熵是一组共轭变量，等温过程和等熵过程也可以视为“共轭”的一对过程。
问题描述 本例将干燥模拟扩展到潮湿模拟，引入水汽和云，我们将模拟在饱和大气环境中上升的饱和气团。
模拟和结果 (env ) jihenghu@canopy:~/bryanproj/build [14:12] $ cmake ~/canoe -DTASK=bryan . $ make -j8 $ cd .. (env ) jihenghu@canopy:~/bryanproj [14:14] $ ln -s build/bin/bryan.release . $ cp build/bin/bryan.inp . $ ln -s build/bin/combine.py .
开始模拟
$ mpirun -n 4 ./bryan.release -i bryan.inp time=1.0000e+03 cycle=7734 Terminating on success cpu time used = 2.459308e+02 (s)
模拟了1000s的变化，使用python进行输出的合并
$ python combine.py -o test $ ncdump -h bryan-test-main.nc netcdf bryan-test-main { dimensions: time = UNLIMITED ; // (11 currently) x1 = 200 ; x1f = 201 ; x2 = 400 ; x2f = 401 ; x3 = 1 ; variables: float time(time) ; time:axis = "T" ; time:long_name = "time" ; time:units = "seconds" ; float x1(x1) ; x1:axis = "Z" ; x1:units = "meters" ; x1:long_name = "Z-coordinate at cell center" ; float x1f(x1f) ; x1f:axis = "Z" ; x1f:units = "meters" ; x1f:long_name = "Z-coordinate at cell face" ; float x2(x2) ; x2:axis = "Y" ; x2:units = "meters" ; x2:long_name = "Y-coordinate at cell center" ; float x2f(x2f) ; x2f:axis = "Y" ; x2f:units = "meters" ; x2f:long_name = "Y-coordinate at cell face" ; float x3(x3) ; x3:axis = "X" ; x3:units = "meters" ; x3:long_name = "X-coordinate at cell center" ; float rho(time, x1, x2, x3) ; rho:units = "kg/m^3" ; rho:long_name = "density" ; float press(time, x1, x2, x3) ; press:units = "pa" ; press:long_name = "pressure" ; float vel1(time, x1, x2, x3) ; vel1:units = "m/s" ; vel1:long_name = "velocity" ; float vel2(time, x1, x2, x3) ; vel2:units = "m/s" ; vel2:long_name = "velocity" ; float vel3(time, x1, x2, x3) ; vel3:units = "m/s" ; vel3:long_name = "velocity" ; float r0(time, x1, x2, x3) ; float vapor(time, x1, x2, x3) ; vapor:units = "kg/kg" ; vapor:long_name = "mass mixing ratio of vapor" ; float mse(time, x1, x2, x3) ; mse:units = "J/kg" ; mse:long_name = "moist static energy" ; float qtol(time, x1, x2, x3) ; float rh(time, x1, x2, x3) ; float temp(time, x1, x2, x3) ; temp:units = "K" ; temp:long_name = "temperature" ; float theta(time, x1, x2, x3) ; theta:units = "K" ; theta:long_name = "potential temperature" ; float theta_e(time, x1, x2, x3) ; float theta_v(time, x1, x2, x3) ; // global attributes: :Conventions = "COARDS" ; :history = "Fri Nov 24 14:34:47 2023: ncks -A ./bryan.out3.nc ./bryan.out2.nc" ; :history_of_appended_files = "Fri Nov 24 14:34:47 2023: Appended file ./bryan.out3.nc had no \"history\" attribute\n" , "" ; :NCO = "netCDF Operators version 4.8.1 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)" ; }
温度场、位温场的的演变：
源码解读 头文件 除了之前介绍的几个类库，这里新增了几个，这里使用注释的方式简介一下
#include <athena/bvals/bvals.hpp> #include <application/application.hpp> #include <application/exceptions.hpp> #include <index_map.hpp> #include <climath/core.h> #include <climath/interpolation.h> #include <climath/root.hpp> #include "bryan_vapor_functions.hpp"
输出诊断量 int iH2O, iH2Oc; Real p0, grav;
注册输出参数：temperature, potential temperature, virtual potential temperature, moist static energy, relatve humidity, total mixing ratio of water。
void MeshBlock::InitUserMeshBlockData (ParameterInput *pin) AllocateUserOutputVariables (7 ); SetUserOutputVariableName (0 , "temp" ); SetUserOutputVariableName (1 , "theta" ); SetUserOutputVariableName (2 , "theta_v" ); SetUserOutputVariableName (3 , "mse" ); SetUserOutputVariableName (4 , "theta_e" ); SetUserOutputVariableName (5 , "rh" ); SetUserOutputVariableName (6 , "qtol" ); }
Grid iterations to assign the value before ouput, Most of the diagnosis variables are members of Class Thermoddynamics, use a pinter to refer the instance.
void MeshBlock::UserWorkBeforeOutput (ParameterInput *pin) auto pthermo = Thermodynamics::GetInstance (); for (int k = ks; k <= ke; ++k) for (int j = js; j <= je; ++j) for (int i = is; i <= ie; ++i) { user_out_var (0 , k, j, i) = pthermo->GetTemp (this , k, j, i); user_out_var (1 , k, j, i) = pthermo->PotentialTemp (this , p0, k, j, i); user_out_var (2 , k, j, i) = user_out_var (1 , j, i) * pthermo->RovRd (this , k, j, i); user_out_var (3 , k, j, i) = pthermo->MoistStaticEnergy (this , grav * pcoord->x1v (i), k, j, i); user_out_var (4 , k, j, i) = pthermo->EquivalentPotentialTemp (this , p0, iH2O, k, j, i); user_out_var (5 , k, j, i) = pthermo->RelativeHumidity (this , iH2O, k, j, i); auto &&air = AirParcelHelper::gather_from_primitive (this , k, j, i); user_out_var (6 , k, j, i) = air.w[iH2O] + air.c[iH2Oc]; } }
AirParcel类定义在air_parcel.hppAirParcel类包含了数种介质类型的定义及其单位的换算函数。介质使用指针表示，一共6种：
Real *const w; Real *const c; Real *const q; Real *const x; Real *const t; Real const *d;
分别表示对应物质的浓度值，可以采用air.w[iH2O]的形式引用对应的水汽含量；enum class Type { MassFrac = 0, MassConc = 1, MoleFrac = 2, MoleConc = 3 };
AirParcel &ToMassFraction () ;AirParcel &ToMassConcentration () ;AirParcel &ToMoleFraction () ;AirParcel &ToMoleConcentration () ;
namespace:: AirParcelHelper定义了几个用于更新<void distribute_to_primitive()>、获取<AirParcel/AirColumn gather_from_conserved()>一个或一段气柱的AirParcel信息的函数。此处用于获取水含量信息。
初始化Mesh类 本个例依然不考虑额外的Forcing， 所以没有注册source function。
void Mesh::InitUserMeshData (ParameterInput *pin) auto pindex = IndexMap::GetInstance (); grav = -pin->GetReal ("hydro" , "grav_acc1" ); p0 = pin->GetReal ("problem" , "p0" ); iH2O = pindex->GetVaporId ("H2O" ); iH2Oc = pindex->GetCloudId ("H2O(c)" ); }
初始化MeshBlock类 依然是先获取热力学的实例指针*pthermo，以便对热力学参数机型操控；block中设定的初始参数: 地表温度，气团中心坐标，椭圆轴长，温度梯度参数，气团初始时刻的水汽混合比 vapor mixing ratio。
void MeshBlock::ProblemGenerator (ParameterInput *pin) auto pthermo = Thermodynamics::GetInstance (); Real Ps = p0; Real Ts = pin->GetReal ("problem" , "Ts" ); Real xc = pin->GetReal ("problem" , "xc" ); Real zc = pin->GetReal ("problem" , "zc" ); Real xr = pin->GetReal ("problem" , "xr" ); Real zr = pin->GetReal ("problem" , "zr" ); Real dT = pin->GetReal ("problem" , "dT" ); Real qt = pin->GetReal ("problem" , "qt" );
AirParcel air (AirParcel::Type::MassFrac) ; air.w[iH2O] = qt; air.c[iH2Oc] = 0. ; air.ToMoleFraction (); qt = air.w[iH2O];
遍历格点，（1）首先使用地表气压和地表温度插值出第一层（is）和地表(0)中间高度处的气温和气压，水汽，也就是 pcoord->dx1f(is) / 2.高度处的空气参数，注意is表示其实下标；
for (int k = ks; k <= ke; ++k) for (int j = js; j <= je; ++j) { air.w[iH2O] = qt; air.w[IPR] = Ps; air.w[IDN] = Ts; air.c[iH2Oc] = 0. ; pthermo->Extrapolate (&air, pcoord->dx1f (is) / 2. , Thermodynamics::Method::ReversibleAdiabat, grav); for (int i = is; i <= ie; ++i) { AirParcelHelper::distribute_to_conserved (this , k, j, i, air); pthermo->Extrapolate (&air, pcoord->dx1f (i), Thermodynamics::Method::ReversibleAdiabat, grav); } }
Extrapolate()函数的定义如下,此处用来构建一个满足固定大气模型的廓线：
void Thermodynamics::Extrapolate ( AirParcel * qfrac, Real dzORdlnp, Method method, Real grav = 0. , Real userp = 0. ) const Construct an 1d atmosphere Parameters [in,out] qfrac mole fraction representation of air parcel [in] dzORdlnp vertical grid spacing [in] method choose from [reversible, pseudo, dry, isothermal] [in] grav gravitational acceleration [in] userp user parameter to adjust the temperature gradient
这个怎么实现呢，实际上Extrapolate调用了rk4IntegrateLnp()函数进行微分方程的数值积分：
#include <algorithm> void Thermodynamics::rk4IntegrateLnp (AirParcel *air, Real dlnp, Method method, Real adlnTdlnP)
rk4IntegrateLnp()进一步调用了calDlnTDlnP()函数来计算温压关系：
Real Thermodynamics::calDlnTDlnP (AirParcel const & qfrac, Real latent[]) const { Real gammad = GetGammad (qfrac); Real q_gas = 1. ; for (int n = 0 ; n < NCLOUD; ++n) q_gas -= qfrac.c[n]; Real f_sig = 1. ; for (int n = 1 ; n <= NVAPOR; ++n) f_sig += qfrac.w[n] * (cp_ratio_mole_[n] - 1. ); for (int n = 0 ; n < NCLOUD; ++n) f_sig += qfrac.c[n] * (cp_ratio_mole_[1 + NVAPOR + n] - 1. ); Real cphat_ov_r = gammad / (gammad - 1. ) * f_sig / q_gas; Real xd = q_gas; for (int n = 1 ; n <= NVAPOR; ++n) xd -= qfrac.w[n]; Real c1 = 0. , c2 = 0. , c3 = 0. ; for (int iv = 1 ; iv <= NVAPOR; ++iv) { c1 += qfrac.w[iv] / xd * latent[iv]; c2 += qfrac.w[iv] / xd * latent[iv] * latent[iv]; c3 += qfrac.w[iv] / xd; } return (1. + c1) / (cphat_ov_r + (c2 + c1 * c1) / (1. + c3)); }
这个函数的来源与Cheng Li博士期间的工作：《Moist Adiabats with Multiple Condensing Species: A New Theory with Application to Giant-Planet Atmospheres》 ，其提出了适用于多种可凝结气体共存情形下的湿绝热温压关系（经典理论只适用于水蒸气存在的地球大气，不适用于巨行星的大气廓线构建），形式是这样的，大家感兴趣可以去研读：
Moist Adiabats T-P with Multiple Condensing Species
最后这部分是再初始场中叠加Guassian气团的温湿度异常场，然后计算水的相变平衡点。
for (int k = ks; k <= ke; ++k) for (int j = js; j <= je; ++j) for (int i = is; i <= ie; ++i) { Real x1 = pcoord->x1v (i); Real x2 = pcoord->x2v (j); Real L = sqrt (sqr ((x2 - xc) / xr) + sqr ((x1 - zc) / zr)); if (L < 1. ) { auto &&air = AirParcelHelper::gather_from_conserved (this , k, j, i); air.ToMoleFraction (); Real temp_v = air.w[IDN] * pthermo->RovRd (air); temp_v *= 1. + dT * sqr (cos (M_PI * L / 2. )) / 300. ; Real temp; int err = root (air.w[IDN], air.w[IDN] + dT, 1.E-8 , &temp, [&pthermo, &air, temp_v](Real temp) { air.w[IDN] = temp; auto rates = pthermo->TryEquilibriumTP_VaporCloud (air, iH2O); air.w[iH2O] += rates[0 ]; air.c[iH2Oc] += rates[1 ]; return temp * pthermo->RovRd (air) - temp_v; }); if (err) throw RuntimeError ("pgen" , "TVSolver doesn't converge" ); air.w[IDN] = temp; auto rates = pthermo->TryEquilibriumTP_VaporCloud (air, iH2O); air.w[iH2O] += rates[0 ]; air.c[iH2Oc] += rates[1 ]; AirParcelHelper::distribute_to_conserved (this , k, j, i, air); } } }
TODO
