2016年8月4日木曜日

C++による気象関数




 私にはC++は難しい、関数を作ってもファイルからデータをstringで読み込んで、charに変換して、さらにdoubleに変換してやっと計算ができるみたいです。
 ハードルが高すぎるのでサンプルプログラムを作るのはあきらめました。
 SSIは500hPaの飽和相当温位から850hPaの相当温位の差を求めればよいので作成しませんでした。
 SSIは温度差ですが、飽和相当温位と相当温位の差はCpをかければエネルギーの差になります。
 単位系は普通の「物理」にある単位系にしました。
 状態方程式PV=nRTRが窒素とヘリウムで値が変わらないようにしました。
 気象では1kgあたりの状態方程式を使いますのでRの値が違ってしまう不具合を避けました。
 いろいろ書きたいのですが・・
 まず、関数をみましょう。
 不具合があっても責任は取りません。個人の責任で使用するのはかまいません。
 
 

template <typename TYPE>
TYPE Abs(TYPE a);
//絶対値を返す関数

double e(double td);
//露点温度(td)から水蒸気圧(e)を求める

double Td(double rh,double t);
//相対湿度()と気温()から露点温度を求める

double g(double hight);
//高さ(hight)毎の重力加速度(g)を求める

double Tw(double temp,double td,double press);
//気温(temp)露点(td)気圧(press)から湿急温度(tw)を求める
double pt(double hight,double temp);
//高さ(hight)と気温(temp)から温位(pt)を求める cf気温と気圧から求めるられる温位より正確です
double ept(double hight,double temp,double td,double press);
//高さ(hight,気温(temp,露点温度()と気圧から相当温位(ept)を求める cf普通求める相当温位より正確です

int main(){
/*????
ゴメン
????*/
}

template <typename TYPE>
TYPE Abs(TYPE a)
{
return a<0 ? -a : a;
}

double e(double td){
static double const e0=6.11;//0℃の水の水蒸気圧
static double const LR0=5374.48884;// L/R0Lは昇華熱R0は水蒸気のガス定数 水蒸気は理想気体ではない
static double const ek=2.791233047;
static double const t0=273.15;//0℃の絶対温度
double LR=LR0-ek*(td-20);//水は理想気体ではないのでL/Rに温度依存性がある
double result=e0*exp(LR*(1/t0-1/(td+t0)));
return result;
}
//水蒸気圧表からekを決定している


double Td(double rh,double t)
{
double t1,t2,tr;
double* t3;
double d,a;
t3=&tr;
cout<<" rh : "<<rh<<" t : "<<t<<endl;
t1=t-100.0;//十分大きな差をとること
t2=t;
*t3=0.5*(t1+t2);
// cout<<"*t : "<<*t3<<endl;
a=rh*e(t)/100;
d=e(*t3)-a;
//cout<<"d : "<<d<<endl;
while(Abs(d)>0.0001){
if( d<0 ){
t1=*t3;
}else{
t2=*t3;
}
*t3=0.5*(t1+t2);
d=e(*t3)-a;
// cout<<"*t3 : "<<*t3<<" d : "<<d<<" a : "<<a<<endl;
}
return tr;
}
//水蒸気圧からTdを追い込んでいる この追い込みは別に関数化できると思うが・・;



double g(double hight)
{
static double const k=0.0000031;//一般気象学の表からkを求めた高度10km以上だと影響がでそう
double result =9.807-k*hight;
return result;
}
//gを固定すると20〜30kmで不自然に温位が高くなる

double Tw(double temp,double td,double press){
static double const L=51012;
static double const cp=29.093;
double t1=temp;
double t2=td;
double *t3;
double tr;
t3=&tr;
*t3=0.5*(t1+t2);
cout<<*t3<<endl;
double d;
double a=cp*(temp+273.15)+L*e(td)/press;
d=a-cp*(*t3+273.15)-L*e(*t3)/press;
cout<<"a : "<<a<<endl;
cout<<"d : "<<d<<endl;
while(Abs(d)>0.01){
// for(int i=0;i<10;++i){
cout<<"d : "<<d<<endl;
if( d<0 ){
t1=*t3;
// cout<<"t2 : "<<t2<<endl;
}else{
t2=*t3;
// cout<<"t1 : "<<t1<<endl;
}
*t3=0.5*(t1+t2);
//cout<<"a : "<<a<<endl;
d=a-cp*(*t3+273.15)-L*e(*t3)/press;
cout<<"*t3 : "<<*t3<<" L*e(*t3)/press : "<< L*e(*t3)/press<<endl;
}
return tr;
}
// pT+Le(Td)/P=pTw+Le(Tw)/PからTwを追い込んでいる 追い込みはTdと同じ方法

double pt(double hight,double temp){
static double const cp=29.093;//N2:28.01 78.088% O2:32.00 20.949% Ar:39.94 0.93% 単原子のArまで考慮した
static double const m=0.028957; //N2:28.01 78.088% O2:32.00 20.949% Ar:39.94 0.93% 単原子のArまで考慮した MKSに統一したほうがよいのだけど・・yoku matigaeru
static double const t0=273.15;//平成23年理科年表P390
//cout<<"m: "<<m<<" g: "<<g(hight)<<" T: "<<temp+t0<< endl;
double result=m*g(hight)*hight/cp + temp+t0;
return result;
}
//エネルギー CpT+mgh をCpで割っている普通の温位より物理的意味は明確


double ept(double hight,double temp,double td,double press){
static double const L=51012.0;//気象関係の値を採用
static double const cp=29.093;
double result=pt(hight,temp)+L*e(td)/press/cp;
return result;
}
//エネルギー CpT+mgh+Le(td)/P をCpで割っている普通の相当温位より物理的意味は明確
//普通の温位は1000Phaまで圧縮するエネルギーを無視するミスを犯している結果SSIがプラスでも不安定なんて変な観測結果になる。


1.重力加速度

 一般気象学 P23
 表2 諸物理量の各高度における値(米国標準大気モデル, 1976)
 (tama は対象を30000m程度までにした)

 高度(Km)  重力加速度(ms−2)
      0       9.807
                    5                            9.791
                  10                            9.776
                  15                            9.761
                  20                            9.745
                  25                            9.730
                  30                            9.715

 0mと20000mを基準とする
 g(h)=9.807−kh
 g(0)=9.807
 g(20000)=9.807−20000k=9.745
 結果 k=0.0000031

2.水蒸気圧
 Es=E0Exp{L/R)(/T0 ー1/T)}
 大気科学講座2 P14
 (注 Rは理想気体のガス状数とは異なる)
 Rは温度に依存する可能性があり
 L/Rが温度に依存すると考える

E06.11 hPa

T=20℃=257.15Kの場合 水蒸気圧125hPa
1.25=6.11Exp{L/R(/273.15-1/253.15)}
ln(1.25)-ln(6.11)=L/R(/273.15-1/253.15 )
1.586783222=L/R(-0.000289235)
L/R=5486.138163


T=20℃=297.15Kの場合 水蒸気圧23.39hPa
L/R=5374.488841


グラフにすると (温度は本来絶対温度で示すべきだった)
 











 




飽和水蒸気圧の値は 大気科学講座2 P13から

LTに依存するL=L(T)
T=20℃=297.15Kを基準として(LR)の計算式(近似式)を作ると

L/R=5374.48884-2.791233047*(T-297.15)
Es=E0Exp{L/R)(/T0 ー1/T)}











高層観測の蒸気圧の値は相対湿度から計算することになる
氷の飽和蒸気圧か水の飽和蒸気圧を使うで値がことなってしまう
軽井沢の地上データを見ると水の飽和蒸気圧を使っている。
天気相談所で確かめておいたほうがよさそうだが・・;
マッイイカ(^^



3.定圧比熱

一般気象学P52
定積比熱Cv717JKkg−1
定積比熱Cp1004JKkg−1


大気科学講座P17
定積比熱Cv716JKkg−1
定積比熱Cp1005JKkg−1

現代熱力学 イリヤ・プリゴジンP33
T=298.15K P=1atm
N2 Cp=29.12 Cv=20.74 (Jmol-1)
O2 Cp=29.36 Cv=20.95 (Jmol-1)
Ar Cp=20.79 Cv=12.47 (Jmol-1)
CO2 Cp=37.11 Cv=28.46 (Jmol-1)

一般気象学P13
地表付近の大気組成(80kmまでほとんど変わらない)
N2 分子量 28.01 存在比 78.088
O2     32.00     20.949
Ar      39.94      0.93
CO2     44.01      0.03

CO2等を無視すると
N2,O2,Arで99.957%








乾燥空気をN2,O2,Arと考えるとCp=29.093 Cv=20.707を採用
平均分子量は28.95712665g/mole
これで計算すると
1kgCpは1004.69
1kgCvは715.09
多分これが正しい値




4.融解熱

一般気象学P58
氷の融解       334 j/g→  6012 j/mole
液体の蒸発      2500 jg→ 45000 jmole

大気講座2 P12
液体の蒸発      2501 jg→ 45018 jmole
氷を水蒸気にする   2834 jg→ 51012 jmole


現代熱力学イリヤ・プリゴジン P131
P=100000Pa=1000hPa
氷の融解 T=273.15  6008 j/mole
水の蒸発 T=373.15 40656 j/mole

気象関係と水の蒸発が大きく違う

理科年表 平成23年
氷の融解  6010  jmole 飽和蒸気圧下
水の蒸発 44000 jmole 25度 P=23.75Torr=3165.875Pa=31.658775hPa
Pa31.66は25℃の飽和蒸気圧31.69に近い(多分誤差範囲)
水の蒸発 40700 jmole 100度 飽和蒸気圧下

氷の融解熱は共通している
水の蒸発圧力や温度に左右されるようだ
根拠ははっきりしないが気象関係の教科書に従ったほうが無難と考える
欲しいのは上空の冷たい環境で水蒸気が氷になる時のエネルギー(潜熱)

潜熱は水蒸気が上空で氷になる潜熱L=51012 j/moleを採用する

クラジュウス・クラペイロンの式(熱力学・統計力学 W.グライナー P65)
平衡である条件

T液体 = T気体
P液体 = P気体
μ液体 = μ気体
ギブスーデュエムの関係式(P59)から独立ではない
もし状態方程式が知られいればTPからμを求めることができる

μ液体(P,T) = μ気体(P,T) (3.12)

(3.12)が成り立つように温度をdTだけ変えると、平衡を保つためには蒸気圧もある値dPだけ変える必要がある。
そうした変化では
液体(P,T) = 気体(P,T)
が成り立っていないといけないギブスーデュエムの関係式
SdT ー VdP + Ndμ = 0
を液体と気体について考えると
液体(P,T) = ー(S液体/N液体)dT + (V液体/N液体)dP
気体(P,T) = ー(S気体/N気体)dT + (V気体/N気体)dP
s液体=(S液体/N液体)
v液体=(V液体/N液体)
s気体=(S気体/N気体)
v気体=(V気体/N気体)
とすれば
dP(v液体 ー v気体)=dT(s液体 ー s気体)
s液体 ー s気体 は液体同じ温度で全部気体した時のエントロピーの変化とみなせるから
s液体 ー s気体=⊿Q’(液体→気体)/T
Q’は液体を気体にするために必要な熱
dP/dT=⊿Q’(液体→気体)/T(v液体 ー v気体)〜ー⊿Q’(液体→気体)/T(v気体)
Q’’を気体を液体にすると出てくる熱とするとー⊿Q’=⊿Q’’だから
dP/dT〜⊿Q’’(液体→気体)/T(v気体)
1moleの水蒸気を考えると⊿Q’’は水にしたの熱量L

moleの水蒸気はPVRvTに従ったとして
dP/dT=LT)(PRvT
dP/P=(L/Rv)dT/T2
lnP(T)=(L/Rv)(ー1/T)const
lnP(273.15=0℃)=ln(6.11)=(L/Rv)(−1/273.15)+const
lnP(T)=(L/Rv)(ー1/T) + ln(6.11) ー (L/Rv)(−1/273.15)
ln(P(T)/6.11) = ー (L/Rv)(1Tー /273.15)
PT)6.11exp{ー (L/Rv)(1Tー /273.15)}
水蒸気の飽和蒸気圧はeesであらわされるようなので
e(T)=6.11exp{ー (L/Rv)(1Tー /273.15)}

Rvは水蒸気の気体定数で大気科学講座2P14でRv=0.4615 J−1K
水の分子量は18だからRv=8.307  JmoleKと気体定数(H23理科年表P367
)R8.314472に近い
Rvを使ってと大気科学講座2P13の氷の蒸気圧を使ってLを求めて見よう

Td=ー10 L51017 j/mole
Td=ー20 L51133 j/mole
Td=ー30 L51080 j/mole

Rvが一定である根拠もないしL51012にかなり近い
大気講座2 P12
液体の蒸発      2501 jg→ 45018 jmole
氷を水蒸気にする   2834 jg→ 51012 jmole

上空の飽和蒸気圧に対応する露点温度はマイナス

下層の水蒸気が上空に行って氷になると仮定してL=51012を採用して良さそう