202. 気象・地学・環境

2011年3月14日 (月)

地震に関する過去の科学記事(バックナンバー)

 この際ですから,地震関係の過去の科学記事を紹介しておきます。

(↑ははーん?なるほど手抜きだね。)

 まず,2006年11/4の「結晶内での弾性波(地震波)」です。

 そして,2009年9/8の「定量的地震学1」,9/15の「定量的地震学2」,9/23の「定量的地震学3」,9/29の「定量的地震学4」,11/17の「定量的地震学5」,11/19の「定量的地震学6」です。

 この後,このシリーズをなぜ途中で中断したのかという理由はよく覚えていませんが,4から5の間があいてるのは10月から介護ヘルパー2級をとるために学校に通っていた時期と丁度一致しています。

 さらには,,2009年12/10の「震源の探知(大森公式等)」があります。

 最後の小記事については丸々転載しておきます。

 (転載開始)↓

 2009年12月10日 「震源の探知(大森公式等)」

今日は,ちょっと手抜きのツナギで2006年11/14の記事 「結晶内での弾性波(地震波)」の続きとして,中学か高校の地学で習った簡単な知識を披瀝してお茶をにごします。

 年末,体調はよくも悪くもないですがヒマ人なりに予定がこみ合っていて貧乏なこともありフリ-マンをエンジョイというわけにはいきません。

 ではまず,記事の再掲から始めます。

 (※再掲)

 今日は等方的とは限らない一般の結晶内での弾性波,特に地震波について考察してみます。

弾性体の密度をρ,応力テンソルを{σjk},それを構成する部分の局所的速度を={uj}とすると,この弾性体が従うべき運動方程式は流体方程式と同じくρd2j/dt2=∂σjk/∂xkで表与えられます。

 

そして,歪み速度テンソル{ujk}をjk≡(∂uj/∂xk∂uk/∂xj)/2で定義すれば,これは{σjk}と同じく反対称の単なる回転を除いた対称テンソルです。これは6個の独立成分を持ちます。

 

線形弾性体近似では,フック(Hooke)の法則:σjk=Cjklmlmが成立するため,先の運動方程式はρd2j/dt2=Cjklm(∂2m/∂xk∂xl)となります。

jklmは,一般に81個の成分を持ちますが,{σjk}も{jk}も対称テンソルなのでjとk,lとmについて対称ですから,独立成分は6×6=36個となります。

  

また,歪みエネルギーWはΔW=σjkΔjkから求まり,対称2次形式W=(1/2)jklmjklmになるので,jklm,kとl,mの交換についても対称であり,結局,その独立成分は(30/2)+6=21個となります。

この方程式の平面波の解をj=u0jexp[i(kr-ωt)]として代入すると,ρω2j=Cjklmklmとなります。これを書き直すと(ρω2δjm-Cjklmkl)um0 となります。

 

この1次方程式が自明でない解を持つためには,3行3列の行列係数:(ρω2δjm-Cjklmkl) (j,m=1,2,3)の行列式がゼロ,すなわち,det(ρω2δjm-Cjklmkl)=0 が成立する必要があります。

 

そして,ω2を未知数としてこの方程式を解けば,ω2の3個の解が得られ,それらはの関数となります。

ところで,もしも結晶が等方性弾性体であれば,歪みエネルギーW=(1/2)jklmjklmは座標系の回転に対して不変なスカラーであるはずですが,jkの2次形式の形で得られる独立な不変スカラーはukk2jkjkのみです。

 

そこで,適当な係数λ,μを選んでW=(1/2)λukk2+μjkjkと書けるはずです。そこで,応力テンソル{σjk}={Cjklmlm}はσjk=λullδjk+2μujkと書けます。ここに弾性定数λ,μはラメ(Lame)の定数と呼ばれています。

このときには運動方程式もρd2/dt2(λ+μ)∇(∇)+μ∇2と,簡単になります。

 

先に挙げた1次方程式の係数の行列要素はρω2δjm-Cjklmkl(ρω2-μ2jm(λ+μ)jmとなります。

 

の向きをx軸の正の向きに取ると,k1=k=||,k2=k30 ですから,行列[(ρω2-μ2jm(λ+μ)jm]は対角成分が[ρω2(λ+2μ)k2]と2つの(ρω2-μk2)の対角行列になります。

 

したがって,det(ρω2δjm-Cjklmkl)=0 の解はVP2≡(ωP/k)2(λ+2μ)/ρ,VS2S/k)2=μ/ρとなります。

 

それぞれの速度に対応する平面波は,速度Pで波の伝播方向への振動である"縦波=P波"と,速度VSで波の伝播方向に垂直な方向の振動である"横波=S波"です。

 

しかし,もしも異方性の弾性体の場合ならω>0 の3つの解ω()はの関数として1次の同次式ではあっても,単なる1次関数になるとは限らず,弾性波は一般に分散性の波でもあると思われます。

 

そこで,伝播速度は"単一波の速度=位相速度"ではなく,群速度=∂ω/∂で与えられると考えられます。

ところで,σjkやujk[σ]=t112233233112)のように6次元の列ベクトルで表わすボイト(Vogit)の表記という表記法があります。

 

この表記で表現すると,フックの法則は[σ]=C[u]と簡単になります。ここでCは6行6列のテンソル行列です。

 

また,このとき,特に結晶が等方性弾性体ならCの成分のうちで独立なのはやはり2つだけで,先に述べたラメの定数はλ=C12,μ=C44=(C11-C12)/2と表わされます。

例えば結晶が立方体構造をしている立方晶系では独立成分は3つです。つまり,C11=C22=C33,C12=C23=C31,かつ4~6行の成分は4~6列しか成分のない対角行列であってC44=C55=C66であり,結局,この結晶構造を示すにはC11とC12とC44の3つだけあれば十分である,ということになります。

この条件ではx軸をの向きに取り,先のdet(ρω2δjm-Cjklmkl)= 0 でのCjklmをボイトの表記でCを表わせば,ξ=(ω/k)2,a=C11/ρ,b=C44/ρ,c=C12/ρとして(ξ-a)(ξ-b)2=0 ですから,解はξ=a,bとなります。

 

つまり速度をVP,VSとすると,VP2=C11/ρとVS2=C44/ρの2つの速度が得られます。これは等方性弾性体でのC44=μ,C11=λ+2μのケースと一致しています。

一方,六方晶系ではC11=C22で,C13=C23,また4~6は対角行列でC44=C55,C66=(C11-C12)/2ですから,結局のところC11,C12,C13,C33,C44の5つだけがあれば十分となります。

 

結果だけ書くと,VP2=C11/ρ,およびVS2=C44/ρ,(C11-C12)/(2ρ)となり,"横波=S波"には2つの速度があるので地震の観測ではS波のほうに2重の波が観測されると予測されます。

 

ただし,異方性の結晶では一般に波は完全にS波とP波になるように行列が対角形になるとは限らないので,これらの計算においては偶々波の進行方向が結晶の対称軸と一致したり直交していたりする特別なケースだけしか扱っていません。

 

一般的な扱いについては,暇があってその気になればまた記事にするかもしれません。

参考文献;ランダウ=リフシッツ 著「弾性理論」(東京図書),角谷典彦 著「連続体力学」(共立出版)

(再掲終了※)

 と書きましたが,ここからもっと身近な地震の話題へとつなげます。

等方性媒質を仮定すると,P{(λ+2μ)/ρ}1/2,VS=(μ/ρ)1/2ですからVP>VSです。

そこで,震源Oから地震を感知する場所:AまでP波,およびS波が到着するまでの時間を,それぞれTP,およびTSと書けばTP=∫OA(1/VP)ds,TS=∫OA(1/VS)dsです。

 

P>VSですから,TS>TPとなります。

つまり,まず縦揺れのP波だけが到着しその後に横揺れS波も到着して両方が重ね合わされた大きい揺れが始まるのですね。

特にOA間の到るところで媒質が同じ均等な弾性体であれば,この区間でP,VSが一定です。

 

そこで,このときにはOA間の波の進行距離(直線距離でなくてもよい)をdAOAdsとするとTPA/VP,TSA/VSです。

これから,引き算するとTS-TP(A/VSA/VP)=A(1/VS-1/VP)です。

 

それ故,TSP≡TS-TP,k≡(1/VS-1/VP)-1=VPS/(VP-VS)と定義すればA=kTSPという比例関係の公式を得ます。

 この式は発見者の大森房吉氏(1918)の名前を取って,大森公式(Ohmori's law)と呼ばれます。

 

    (下図は中越地震本震の地震計記録からです。)

 

     

こうした等方的で均質な媒質を伝わる場合,弱いP波だけが来て「あ,ひょっとして地震かな?」と思ってから本格的で大きな揺れがくるまでの初期微動の時間をTSPとするとその地点Aから震源Oまでの距離AはTSPに比例します。

初期微動時間が長いほど震源までの距離(地理的な距離だけでなく震源深さまでも含めた距離)が長いということが言えます。もしも震源が直下にあれば縦揺れは上下震動,横揺れは左右震動になります。

さて,標高hが違う3地点A,B,Cがありその座標がそれぞれ(xA,yA,hA),(xB,yB,hB) (xC,yC,hC)であるとき,上記のような公式に基づいて初期震動時間の観測から観測点から震源までの距離A,B,Cが全てわかったとします。

このとき,未知の震源Oの座標を(xO,yO,hO)と仮定すれば地震波が全て直進ならA2(xA-xO)2+(yA-yO)2+(hA-hO)2,B2(xB-xO)2+(yB-yO)2+(hB-hO)2,C2(xC-xO)2+(yC-yO)2+(hC-hO)2が成立します。

これは3つの未知数xO,yO,hOに対する3つの連立2次方程式ですから解くことが可能です。つまり,正確に震源までの距離を観測できる点が3個あれば原理的には震源の位置がわかります。これは3点法と呼ばれるものですね。

普通,震源の高さhOは海抜で測って負の数であり震源は地中または海中にあります。

大森公式A=kTSPは地震に対する公式で大森係数と呼ばれるk=VPS/(VP-VS)は8(km/s)程度の値ですが,この式は雷における音と光の2つのように速さの異なる2つの信号が届く現象では全く同じ公式が成立します。

ただし,雷の場合は光速がc~30万km/s,空気中の音速がv~340m/sでc>>vなのでk=cv/(c-v)~cv/c=vですから,事実上A~vTなのでほとんど光速は関係なしです。(Tは光が見えてから音が聞こえるまでの時間です。)

 

例えば,ピカッと光ってから雷の音が聴こえるまでの時間が1秒なら,自分から340m程度離れたところ(空中かもしれません)で"放電=落雷"があったと推測されます。

    

PS:>私のマノン・レスコーたちへ 

 病気とか事故に会っていなければどこで誰と遊んでいてもいいけれど,連絡つかない状況だととにかく心配です。大丈夫かなあ。。。。)

 ("マゾ, ロリコン, 準デブ専, メガネフェチ=変態"のTOSHIより。。)

(T.ウッズも聖人君子ではなく適当にストレスを解消しているらしいので安心しました。あれだけの収入を得る能力があれば大勢の家族を養えますね。

 性欲のある健康な男だし,大有名人で品行方正キャラが売り物なら日本だと風俗通いもままならないので女遊びは囲い込むしか仕方ないとして,もしもプロテスタントならマドンナの養子のように大勢を扶養するのはむしろ義務かな?)

PS2:ドクちゃんって日本が好きなんだ。。??

 国ではなくて故郷(くに)が好きなんでしょ?

 http://news.aimu-net.com/read.cgi/wildplus/1257045443/

PS3:「陳情政治」からいきなり「非陳情政治」へと革命的に移行するのではなく,途中で1人(小沢さん)に権力が集中するシステムを経る。。

 共産革命でも過渡期ではプロレタリア独裁(現実には中国では毛沢東,ロシアではレーニンの1人独裁)です。それはそれでマヌーバとしてありなのでしょうが独裁者がプラトン的(プラトニック)に正しい哲人でいるうちはいいけど,うまく移行できないと"元の木阿弥"か,もっと悪くなります。

 小沢氏が中国に大挙して出かけていって米国にブラフかけて「ポッポ政権」を裏で後押ししてるのかも知れないけど,うまくいくのかねえ。基地問題。。

 朝っぱらから女子アナが「吉宗は何将軍?」と聞かれて「暴れん坊将軍?」と答えていました。それでもいいんだろうけど,ニュース番組なので一応「米(コメ)将軍だろ?」とツッコミたくなりました。。。アハハ

 ↑(転載終了)

| | コメント (0) | トラックバック (0)

2010年12月20日 (月)

水滴の成長と蒸発(2)

 さて水滴の成長と蒸発の続きです。

 

 色々としあわせな状態にかまけてはいますが,基本的にやることは変わっていません。

 

 水滴が媒質を含む希薄溶液から成る場合について記述します。 

§3.希薄溶液の水滴の凝結と蒸発

①定常状態(平衡状態)の熱と質量の関係

 

 蒸発の潜熱をLe(単位質量当たり)とすると,凝結によって単位時間当たり解放される熱量はLe(dm/dt)です。

そこで,この潜熱を熱源とする=0での熱方程式は(dq/dt)­0=-∫r=ahdS+Le(dm/dt)­0 (31)となります。ただし,aは水滴半径,qは水滴の熱量です。

この場合,熱平衡(定常)の条件:(dq/dt)­0=0はLe(dm/dt)­0=∫r=ahdS,またはVを水滴の体積として凝結による潜熱の解放が均質に生じるとすれば,∇h=Le(dm/dt)­0/V=∇hです。

 あるいは,h=-k∇T(25)より

r=ahrdS=-∫r=aa(∂T/∂r)dS=-4πaka(T-Ta)

なので,水滴表面温度Taが一定条件下では,

 

e(dm/dt)­0=∫r=ahdS=4πaka(Ta-T)

 

と表現できます。

 

ただし,前述したようにkaは水滴の熱伝導度です。

最後の式は,(蒸発に必要な熱量)=-Le(dm/dt)­0=4πaka(T-Ta),または(凝結によって解放される熱量)=Le(dm/dt)­0==4πaka(Ta-T)=(伝導により流出する熱量)を示しています。

このLe(dm/dt)­0=4πaka(Ta-T)なる式は,(dm/dt)­0=(4πaka/Le)(Ta-T)(32),Ta=T+{Le/(4πaka)}(33)と書き直すことができます。

故に,蒸発水滴:(dm/dt)­0<0 の温度はまわりよりも低く,凝結水滴:(dm/dt)­0>0 の温度はまわりよりも高いといえます。

Clausius-Clapeyronの公式

 

 さて2相(液相,気相)平衡状態の飽和水蒸気esat(T)の温度保存式を求めます。

 熱力学によれば,平衡の条件は液相,気相のGibbs自由エネルギーについてGw=Gv(34),またはμw=μvが成立することです。

 

 Gw,Gvは1モル当たりの水,水蒸気のGibbs自由エネルギーです。

 

 また,μwvは1モル当たりの水,水蒸気の化学ポテンシャルですが,Gw=μw,Gv=μvであり,同じものを別記号で表わしているだけです。

 温度TがT+dT,圧力pがp+dpになっても(34)の関係が維持されるためには,その際のGw,Gvの増加分dGw,dGvについてもdGw=dGv(35)が成立することが必要です。

 ところで,熱平衡状態ではdu=Tds-pdvでG=u+pv-Tsですから,dG=-sdT+vdpです。

 

 ただし,u,sはそれぞれ1モル当たりの内部エネルギー,エントロピー,vは1モル当たりの体積です。

 したがって,dGw=dGv (35)は-swdT+vwdp=-svdT+vvdpを意味します。

 

 つまりdp/dT=(sv-sw)/(vv-vw)(36)を意味します。

 ところが,T(sv-sw)=Lew (37)ですから,これは

  dp/dT=Lew/{T(vv-vw)}とも書けます。

同じ1モルでは水の体積は水蒸気の体積よりはるかに小さいので,この式の右辺で水蒸気の体積vvに比して水の体積vwを無視します。

  

左辺のpp=esat(T),右辺のvvにvv=RT/esatを代入すると,

  

(1/esat)(desat/dT)=Lew/(RT2)(38)を得ます。

 

これが,有名なClausius-Clapeyronの公式(クラウジウス・クラペイロンの公式)です。

 

念のため,改めて追記するとLeは単位質量当たりの蒸発の潜熱,Mwは水の分子量です。

 

     (↓下図はネット検索で入手した図の転載です。)

  

  

(注):エンタルピー:h≡u+pvを導入するとdu=Tds-pdvよりdh=Tds+vdpです。

 

 2相平衡にあるときの圧力pは,dp/dT=(sv-sw)/(vv-vw)(36)を満たしますから,dh=Tds+vdp=Tds+v(dp/dT)dTと書けます。

そこで,hv-hw=T(sv-sw)+∫v(dp/dT)dTですが,相変化が定温(dT=0)の下で移行する場合の潜熱はエンタルピーの変化として表現されます。

 

つまり,hv-hw=T(sv-sw)=Lewです。

したがって,p=esatに対するClausius-Clapeyronの公式:

(1/esat)(desat/dT)=Lew/(RT2)(38)は,

 

(1/esat)(desat/dT)=d(lnesat)=(hv-hw)/(RT2)とも表現されることがわかります。(注終わり)※

 蒸発の潜熱:LeがTに依らず一定値を取ると仮定して,

 (1/esat)(desat/dT)=Lew/(RT2)(38)を積分すると,

 esat(Ta)=esat(T)exp[{Lew/(RT)}(Ta-T)/Ta](39)

 が得られます。

Kelvinの公式,Raoultの規則

 

 半径aの純水滴と湿潤空気の混合系で,気相と液相が平衡にあるためにはμw=μv(40),かつpw=p+2σ/a(41)が成立する必要があります。

 

 ただし,σは表面張力です。

 

(このブログのバックナンバーで表面張力について言及しているのは2009年4/28の「水の波(2)(浅水波,深水波,表面張力波)」だけですね。

  

 一般の球でない水面の場合にはpw=p+σ/a1+σ/a2です。)

 

 ここで,p,Tの変化の下で,μw/T=μv/Tの関係が維持されるためにはd(μw/T)=d(μv/T)(42)が成立することが必要です。

 

 ところが(∂μw/∂pw)T=vw,[∂(μw/T)/∂T]pw=-hw/T2より,

d(μw/T)=(-hw/T2)dT+(vw/T)dpwです。

 

 ただしhはエンタルピー:h≡u+pvです。

  

 一方,水蒸気のみの1成分系では,Gv=μvであり,

混合系ではGxa=xaa+xvv(xa+xv=1)です。

ここで,aの添字aは(乾燥)空気(air)のa,

vのvは水蒸気(vapor)のvです。

そして,水蒸気の分圧はe=xvpですから,Gv=μv=μv+(T)+RTlneよりμv=μv+(T)+RTlnp+RTlnxvです。

 

これは,μv=μv(T,p,xv)=μv0(T,p)+RTlnxvと書けます。

 

μv0(T,p)は空気が純水蒸気(xv=1)で,その気圧がpのときの化学ポテンシャル(Gibbs自由エネルギー)です。

(∂μv/∂p)T=vv,[∂(μv/T)/∂T]pw=-hv/T2より,d(μv/T)=(-hv/T2)dT+(vv/T)dp+Rd(lnxv)ですから,

 

条件d(μw/T)=d(μv/T)(42)は,(-hw/T2)dT+(vw/T)dpw=(-hv/T2)dT+(vv/T)dp+Rd(lnxv)と書けます。

また,pw=p+2σ/a(41)よりdpw=dp+2d(σ/a)でありxv=e/pですから,{-(hv-hw)/T2}dT+{(vv-vw)/T2}dp-(2vw/T)d(σ/a)+Rd{ln(e/p)}=0 です。

 

つまり,(-Lew/T2)dT+{(vv-vw)/T2}dp-(2vw/T)d(σ/a)+Rd{ln(e/p)}=0 (43)が成立することが必要です。

ここでのeはe=esat(T)を意味しますが,今の混合系では定温,定圧下で半径aの変化に対してe=esat(T)が変化すると予想されるため,e=esat(T)を,改めてe=ea sat(T)と書きます。

定温,定圧下ではdT=dp=0なので,(43)は-(2vw/T)d(σ/a)+Rd{ln(easat/p)}=0(44)となります。

 

p=一定ですから,ln{easat(T)}=2vwσ/(RTa)+const.より,

a sat(T)=esat(T)exp{2Mwσ/(RTρwa)}(45) を得ます。

 

これがKelvinの公式です。

 一方,溶解液に対するRaoult(ラウール)の規則から,e'satをこの溶液の溶媒(水)の飽和蒸気圧とし,mw,ms,mを溶媒,溶質,(溶媒+溶質)の質量,Mw,Msを溶媒,溶質の分子量とすると,

 

 iをvan'tHoff factorとして, 

 e'sat/esat=nw/(nw+ins)=mws/(mws+imsw)

 =[1+imsw/{(m-ms)Ms}]-1(46)と書けます。

(45),(46)よりesatをesatと書けば,e'asat(T)=

sat(T)exp{2vwσ/(RTρwa)}[1+imsw/{(m-ms)Ms}]-1(47)

です。

④溶解液の蒸発,凝結

 

 ここまでの公式を以下のように要約します。 

すなわち,a(da/dt)0={Dvw/(ρwR)}{e/T-e'asat(Ta)/Ta}(7)',a(da/dt)­0={ka/(Leρw)}(Ta-T)(32)',

および,

 

sat(Ta)=esat(T)exp[{Mwe/(RT)}(Ta-T)/Ta](39)', e'asat(Ta)=esat(Ta)exp{2Mwσ/(RTaρwa)}[1+imsw/{(m-ms)Ms}]-1 (47)'です。

そして,(32)'をTa=T(1+δ),

δ≡{Leρw/(ka)}a(da/dt)­0 (48)と表現すれば,

 

(7)',(39)',(47)',(48)によって,

 

a(da/dt)0={Dvw/(ρwR)}(e/T-esat(T)exp[Mweδ/{RT(1+δ)}+2Mwσ/{RTρwa(1+δ)}][1+imsw/{(m-ms)Ms}]-1) (49)が得られます。

δ={Leρw/(ka)}a(da/dt)­0ですが,δ<<1なので

exp[Mweδ/{RT(1+δ)}~ 1+Mweδ/(RT),

かつexp[2Mwσ/{RTρwa(1+δ)}]~ exp{2Mwσ/(RTρwa)}

と近似できます。

さらに,

 

exp{2Mwσ/(RTρwa)}[1+imsw/{(m-ms)Ms}]-1~ 1+y(50),

y≡2Mwσ/(RTρwa)-im/{Ms(m-ms)}(51)と置くことで,

近似式:e'a sat(Ta)=e sat(T)(1+y)(52)

 

を得ます。

以上から,

 

a(da/dt)0={Dvw/(ρwR)}{esat(T)/T}][e/esat(T)-{1/(1+δ)}{1+Mweδ/(RT)}(1+y)]です。

さらに,δの2次以上とδyを無視すると,

 

a(da/dt)0={Dvwsat(T)/(ρwRT)}[e/esat(T)-(1+y)-δ{Lew/(RT)-1}]

 

です。

そこで,δ={Leρw/(ka)}a(da/dt)­0より,

a(da/dt)0wRT/{Dvwsat(T)}+{Leρw/(Ta)}{Lew/(RT)-1}]

~{e-esat(T)(1+y)}/esat(T),

 

あるいは,

a(da/dt)­0~[{e-e'sat(T)}/esat(T)]/[{Leρw/(Ta)}{Lew/(RT)-1}+ρwRT/{Dvwsat(T)}](53)

を得ます。

 ただし,e'sat(T)=esat(T)exp{2Mwσ/(RTρwa)}[1+imsw/{(m-ms)Ms}]-1(54)です。もちろん,m=4πa3ρ/3です。

 ここまでは,もっぱら(dm/dt)­0,a(da/dt)­0の静止水滴に関する話でしたが,運動水滴の場合には単にDv→ Shv/2,ka→ Nuaとすればいいだけです。

内容を思い出すのにてこずりましたが,ノートはここで終わっているので,この項目についてはこれで終わりにします。

  

PS:一昨年の暮れ頃には観音様の御出現,今年は天使の御出現と,ここのところ浮世離れしたことが続いています。

 

翌日(21日)も天使がそばにいてとてもよかったです。

 

夜には今年最後の手話講習(初心コース)に出席,終わって21時頃から椎名町駅近くで応用,上級コース合同の忘年会でした。

 

初対面のろう(あ)の人との手話もはずみましたが,翌日が早いので夜中の零時半頃には帰宅してそのまま就寝しました。

 

今日(22日)は,仕事をお休みして朝から眼底出血していた右目を含む両目眼底の総合的検査のため,予約していた帝京病院眼科診療に行ってきます。

 

しかし,今は右目も左目と同じくらい見えるようですから,医者の予想通り幸い自然に出血が止まって治癒しているようで,レーザー手術もしないで済みそうです。

| | コメント (0) | トラックバック (0)

2010年12月12日 (日)

水滴の成長と蒸発(1)

 前の一連の記事のプルーム上昇式について書いていた1983年8月の覚書きノートには続きがあります。それは水滴の成長と蒸発という項目です。

これも1000ページ近い"英文のANLモデル(アルゴンヌ国立研究所模型;後のSACTI)のドキュメントを全部翻訳して内容を完全に把握する"という当時のジョブの一環として勉強したものです。

冷却塔(cooling tower)からの温排気排出の影響には,水蒸気による白煙発生の他に,それが周辺大気で冷却され凝結した水滴の飛散による濃霧,また冷却水が海水のように塩分などを含む場合には飛散水滴による農作物等への塩害etc.も考えられるからです。

 

そこで,この項目についても覚書きとして記述しておきます。

 

(当時は,仕事の都合上,タイムリミットもあって,とにかく当座の計算アルゴリズム等に必要な式に到達するのが目的だったので,途中のプラス,マイナスなど細かいところはアバウトだったりテキトーだったりなところもあるので,あとで修正が必要なのです。)

以下は本文です。

§1.拡散流束による水滴質量の増減

 水蒸気の密度をρv,その密度流束(局所流束密度)をvとします。

この局所流束密度vは,移流風速(水蒸気を含む空気の平均流速)をとすると,そのゆらぎ(fluctuation)に相当する拡散流束を加えて,

 

v=ρv-Dv∇ρv と表現されます。

 

ただし,Dvは拡散係数です。

この流束を用いた水蒸気の質量保存を示す連続の方程式は,

∂ρv/∂t+∇v=0 or ∂ρv/∂t=-∇vです。

 

そして,v=ρv-Dv∇ρvによりv=∇(ρv-Dv∇ρv)です。

 

しかし,この水蒸気を含む空気は非圧縮(または,高々Boussinesq近似)を仮定しているため,∇=0 を満たしますから,

 

これは,v∇ρv-∇(Dv∇ρv) と書けます。

さらに,拡散係数Dv空間位置によらない定数であると仮定すると,

v∇ρv-Dv2ρvです。

結局,上記連続の方程式:∂ρv/∂t=-∇vは,

 

∂ρv/∂t+∇ρv=Dv2ρv (1) を意味します。

 

これは,Fick型の(移流)拡散方程式(diffusion equation)と呼ばれる方程式です。

 一方wを液体の水の密度として半径aの球水滴の質量をmとするとm=∫ρwdV=4πa3ρw/3です。

水滴の中心を原点:Oとするような極座標で考えると,r>aの領域では水蒸気に対して,上述の拡散方程式:∂ρv/∂t+∇ρv=Dv2ρv (1)が成立しますが,r<aの球内は(1)が成立しない特異領域と考えられます。

 そして,空気は水滴表面から内部に浸入することはないので,表面r=aでは速度の表面への垂直成分:ur=0 であるべきですが,現実の水滴表面ではより強い粘着条件:=0 が満たされていると思われます。

 

 つまり,表面r=aでは=0 なる境界条件が満たされると考えます。

 水滴の質量保存から,明らかにdm/dt=-∫r=avrdSですが,

 上記境界条件を要求すると,r=aでは

 v=ρv-Dv∇ρv=-Dv∇ρvなので,

 

 dm/dt=Dvr=a(∂ρv/∂r)r=adS (2)です。

 特に,表面だけでなく全空間で=0 を満たす,水滴に対して静止した空気中でのdm/dtを(dm/dt)0と書けば,場は中心対称ですから,

 

 (dm/dt)0=4πa2v(dρv/dr)r=a (3)と書けます。

他方,=0 なら∂ρv/∂t=Dv2ρvです。さらにt→ ∞ の極限では拡散場は定常(∂ρv/∂t=0)になり,しかも球対称ですから,

 

2ρv=r-12(rρv(r)/dr2=0 です。 

これを,r=aでρv=ρv(a);r=∞でρv=ρの境界条件下で解けば

 

ρv(r)=ρ+a{ρv(a)-ρ}/r (4)を得ます。

これにより,(dρv/dr)r=a={ρ-ρv(a)}/a (5)ですから,

 

=0の静止大気中での水滴の質量保存の式は, 

(dm/dt)0=4πa2v(dρv/dr)r=a=4πaDv-ρv(a)}

と書けます。

水蒸気が理想気体なら,それが従う状態方程式:e=ρvRT/Mwから,

ρ=ew/(RT),ρv(a)=eaw/(RTa)を代入して,

 

(dm/dt)0=(4πaDvw/R)(e/T-ea/Ta) (6)

 

です。

または,m=4πa3ρw/3なる式により,

 

方程式:a(da/dt)0={Dvw/(ρwR)}(e/T-ea/Ta) (7)

 

を得ます。

ただし,aは水滴表面(r=a)の温度,eは水蒸気圧,Rはモル気体定数(R~ 8.31J/(molK)),Mwは水の分子量(Mw~ 18g/mol)です。

しかし,水滴半径aが分子半径と比較してそれほど大きくない場合には,空気や水蒸気を連続体とみなす巨視的流体力学の適用範囲の外に出ると考えられるため,必要な小水滴に対する補正をします。

 

以下,Fuchs(1959)によって与えられた方法に従います。

 拡散方程式とその解は水滴表面からの水蒸気の反跳長さ(recoiling length):Δvを超えないなら正しくないと考えます。

すなわち,r≧a+Δvではρv(r)=ρ+A/r(8)が成立しますが,

a≦r≦a+Δの境界領域では水蒸気流束を構成する分子はMaxwell-Boltzmannの速度分布に従う粒子として気体運動論的に扱うことにします。

水分子の質量をmwとし,分子の平均速度を=(vx,vy,vz),その大きさをv≡||,Boltzmann定数をkB≡R/NA(NAはAvogadro数)とすると,温度がTのときの分子の平均速さは,

  

<v>=(8kBT)1/2/(πmw)1/2=(8RT)1/2/(πMw)1/2 (9)

 

で与えられます。

(注):<v>=∫vexp{-mw2/(2kBT)}d3/∫exp{-mw2/(2kBT)}d303exp{-mw2/(2kBT)}dv/∫02exp{-mw2/(2kBT)}dv です。

ところが,∫0exp(-αv2)dv=(1/2)(π/α)1/2ですから,両辺をαで微分することで,∫02exp(-αv2)dv=(1/4)π1/23/2を得ます。

 

さらに,これをαで微分すると,∫04exp(-αv2)dv=(3/8)π1/25/2です。

一方,∫03exp(-αv2)dv=(1/2)∫0uexp(-αu)du=(1/2)(1/α2)より,∫03exp(-αv2)dv/∫02exp(-αv2)dv=2/(πα)1/2です。

 

<v>=∫03exp{-mw2/(2kBT)}dv/∫02exp{-mw2/(2kBT)}dvの右辺は,上記のα=mw/(2kBT)の場合に当たりますから,

 <v>=(8kBT)1/2/(πmw)1/2=(8RT)1/2/(πMw)1/2です。

一方,二乗平均速度を評価するのなら,

<v21/2=[∫v2exp{-mw2/(2kBT)}d3]1/2[∫04exp{-mw2/(2kBT)}dv/∫02exp{-mw2/(2kBT)}dv]1/2

={3/(2α)}1/2=(3kBT/mw)1/2(3RT/Mw)1/2です。

 

後者はエネルギーの等分配則<mv2/2>=(3/2)kBTから予想される結果に一致しています。

 

(8/π)1/231/2 には大差がないので,<v>と<v21/2は同じオーダーの量です。(注終わり)※

 さて,nを分子濃度,つまり単位体積当たりの分子数とします。

分子の従うMaxwell-Boltzman分布を,

"=(vx,vy,vz)と+d=(vx+dvx,vy+dvy,vz+dvz)を対角線とする速度空間の微小直方体内に,単位体積当たりF(v2)dvxdvydvz個の分子があること"で表現すれば,

 

∫F(v2)dvxdvydvz=n,∫vF(v2)dvxdvydvz/n=<v>

 

です。

すると,速度空間では速度+dの間にあり,かつ座標空間の任意の向きの平面の微小立体角:dΩ=cosθdθdφに衝突する単位時間,単位面積当たりの分子数は,

 

{vsinθdΩ/(4π)}F(v2)dvxdvydvz

=vF(v2)dvxdvydvzsinθcosθdθdφ/(4π)

 

で与えられます。

ただし,この平面がとなす角がθであるように極座標軸を取ってます。

平面の一方の側の単位面積に衝突する単位時間当たりの分子の総数は,これを角度θについて 0 からπ/2まで,φについて 0 から2πまで,さらにvについて 0 から∞まで積分した値です。

0π/2sinθcosθdθ∫0dφ=πですから,この値:[∫vF(v2)dvxdvydvz0π/2sinθcosθdθ∫0π/2dφ]/(4π)は,(1/4)∫vF(v2)dvxdvydvz=(1/4)n<v>に等しいことがわかります。

したがって,(分子数密度流束)=(任意の面に衝突する単位時間,単位面積当たりの分子数)は, 

 

(1/4)n<v> (10) と書けます。

 

右辺の因子:<v>は先に書いたように,

<v>=(8kBT)1/2/(πmw)1/2=(8RT)1/2/(πMw)1/2 (9)

で与えられます。

また,水蒸気の温度を改めてTの代わりにTgと書くと,水蒸気はe=ρwBg/mw=nkBgなる状態方程式に従います。

  

故に,水滴の周囲の水蒸気の分子数密度流束は

=(1/4)n<v>=(1/4){e/(kBg)}(8kBg)1/2/(πmw)1/2

=e/(2πwBg)1/2です。

他方,水滴表面に付着する水蒸気の分子数密度流束をwと書き,平衡状態にあるときのwの水滴の周囲の分子数密度流束wに対する比を凝結定数と呼んでαcとします。

 

すると,=αcです。

そうすると,水滴表面の分子数密度流束は,

=αce/(2πmwBa)1/2 (13)

で与えられることになります。

 

実験によれば,αcの値は,αc≒0.01~ 0.07です。

ただし,水蒸気分子の平均速さは,<v>=(8RTg)1/2/(πMw)1/2でも

分子速度の分布は等方的で,そのベクトル平均はゼロと考えられるので

分子数密度流束wの向きはあらゆる方向について対称です。

水滴表面の水蒸気圧:eが飽和水蒸気圧:esatに等しいときに"定常状態=蒸発平衡"に達するとすれば,

 

平衡時の水滴表面での飽和分子数流束:wsatは, 

sat=αcsat/(2πmwBa)1/2 (14) です。

  

ただし,Taは水滴表面の絶対温度です。

それ故,流入する分子数流束と流出する分子数流束は互いに独立であるとしてTg~ Ta =Tと置けば,水滴に流入する総流束密度の大きさは

 

sat=w-wsat=αc(e-esat)/(2πmwBT)1/2(15)

 

と書けます。

この(15)式は,e<esatなら水滴が蒸発し,e≧esatなら蒸発しないことを意味しています。

さて,水滴中心からの距離rにおける水蒸気の分子数密度をNv(r),密度をρv(r)とすればρv(r)=mwv(r)です。

したがって,r=a+Δにおける水蒸気が,流出入する水蒸気を代表するとすれば,その水滴表面へと向かう凝結流束密度は,

αcwv(a+Δ)<vv>/4=αc<vv>ρv(a+Δ)/4 です。

他方,水滴表面r=aで蒸発に向かう蒸発流束密度は,

αc<vv>ρv(a)/4です。

そこで,表面での総蒸発流束v(a)は,

 

v(a)=4πa2v(a)=πa2αc<vv>{ρv(a)-ρv(a+Δc)} (16)

 

で与えられることがわかります。

一方,冒頭で書いたように,巨視的流体力学の見方では,

r=aでの流速密度はv=-Dv∇ρvなので,総流束は

v(r)=4πar2v(r)=-4πr2(dρv/dr)Dv

です。

ところで,今の気体運動論的考察では,

r≧a+Δρv(r)=ρ+A/r (8)が成立する,

と仮定しているので,

 

r=a+Δにおいてはρv(a+Δv)=ρ+A/(a+Δv)です。

 

故に,(dρv/dr)r=a+Δv=-A/(a+Δc)2ですから,

v(a+Δv)=-4π(a+Δv)2(dρv/dr)r=a+Δvv=4πDvA (17)

を得ます。

"定常状態=蒸発平衡"では,Jv(a+Δv)=Jv(a) (18)であるべきで,

ρv(a+Δv)=ρ+A/(a+Δv)ですから(16),(17)より,

4πDvA=πa2αc<vv>{ρv(a)-ρv(a+Δc)}

=πa2αc<vv>{ρv(a)-ρ-A/(a+Δv)}

が成立します。

したがって,

  

[a2αc<vv>/(a+Δv)+4Dv]A=a2αc<vv>{ρv(a)-ρ},

および<vv>=(8kaa)1/2/(πmw)1/2=(8RTa)1/2/(πMw)1/2 (9)

 

が成立します。

以上から,

 

A=a{ρv(a)-ρ}/[a/(a+Δv)+4Dv/(aαc<vv>)]

=a{ρv(a)-ρ}/[a/(a+Δv)+Dv{2πMw/(RTa)}1/2/(aαc)] (19) を得ます。

これは,冒頭で述べた巨視的流体力学の考察:

(4)ρv(r)=ρ+a{ρv(a)-ρ}/r,

および(dm/dt)0=4πaDv-ρv(a)}から導かれた式,

 

(6)(dm/dt)0=(4πaDvw /R)(e/T-ea/Ta),

または(7)a(da/dt)0={Dvw /(ρwR)}(e/T-ea/Ta)

 

を次のように補正することに相当します。

すなわち,(4)はa{ρv(a)-ρ}→ A=a{ρv(a)-ρ}/[a/(a+Δv)+Dv{2πMw/(RTa)}1/2/(aαc)]として,

 

ρv(r)=ρ+a{ρv(a)-ρ}/r→ (4)'ρv(r)=ρ+A/r

と補正します。

そこで,(6),および(7)は拡散係数をDv→ Dv'=Dv/[a/(a+Δv)+Dv{2πMw/(RTa)}1/2/(aαc)](20)として,

 

それぞれ(dm/dt)0=(4πaDv'Mw/R)(e/T-ea/Ta)(6)',

およびa(da/dt)0={Dv'Mw/(ρwR)}(e/T-ea/Ta)(7)'

に補正します。 

これで小水滴に対して補正された質量保存式を得ました。

 次に,水蒸気を含む空気に対して運動している水滴に対する補正を考察します。

 

 これは,水滴が運動する代わりに静止水滴に対して空気が運動する,つまり≠0 の等価な問題として設定できます。

この場合も,r>aにおける拡散は定常(∂ρv/∂t=0)に達しているとすれば,質量保存の拡散方程式は∇ρv=Dv2ρvです。

 

水滴表面での流束密度がv=-Dv∇ρvなので,

やはりdm/dt=-∫r=avrdS=∫r=av(∂ρv/∂r)dSです。

 ここで,平均通風係数と呼ばれる量fvを,

v≡dm/dt/(dm/dt)0(21)で,

 

 平均輸送係数kvをkv≡-∫r=av(∂ρv/∂r)dS/[4πa2v(a)-ρ}] (22) で定義します。

 

 すると,(2)'と(21),(22)からkv=Dvv/aです。

さらに,無次元のSherwood数Shを次式で定義します。

 

すなわち,Sh≡2kva/Dv=2fv (23)です。

 

こう定義すると,dm/dt=Shv2πa(ρ-ρsat) (24)です。

=0 のときにはfv=1 or Sh=2です。

 

ただし,これらは単に定義した,または言葉を置き換えただけで何ら新しいことを述べていません。文献を読む際の便宜だけでしょう。

§2.水滴中の熱の移動

 kを媒質中の熱伝導度,Tを絶対温度とすると,Fourierによれば熱流束ベクトル:hh=-k∇T (25)で与えられます。

そして,一般的な流体中のエネルギー保存の式は,粘性による散逸が無視できる場合には流体単位質量あたりのエントロピーをsとしてρT(∂s/∂t+∇s)=-∇hと書けます。ρは流体の密度です。

今の場合のように,媒質流体(水,空気)が非圧縮と見なせるときには

=0,かつTds=cpdT(cpは単位質量当たりの定圧比熱)なので,

 

これはρcp(∂T/∂t+∇T)=κ∇2

or ∂T/∂t+∇T=κ∇2T(26) と変形されます。

ただし,κ≡k/(ρcp)でκは熱拡散係数と呼ばれる量です。

 h=-k∇T(25)より,水滴の熱量をqとし,乾燥空気の熱伝導度をkaとすると,

  

 dq/dt=-∫r=ahrdS=∫r=aa(∂T/∂r)dS (27)

 

 です。

特に,0(水滴に対して水蒸気静止)のときのdq/dtを(dq/dt)0と書きます。

 

=0 の熱伝導方程式:∂T/∂t=κ∇2Tのt→∞での定常状態(∂T/∂t=0)の∇2T=r-12(rT)/dr2=0 の解:T=T+a(T-Ta)/rから,(dq/dt)0=4πaka(T-Ta)(28)です。

aが小さい小水滴のときには,先に拡散係数Dhの代わりにDh'を用いるよう補正したのと同じく,熱伝導度kaの代わりにka'≡ka/[a/(a+ΔT)+ka{2πMw/(RTa)}1/2/(aαT+ρcp)](29)を用います。

 

なお,ΔTTは拡散での反跳長さΔv,凝結定数αcに相当する熱伝導(熱拡散)での量です。

また,水滴に対して水蒸気が運動するときには,Nusselt数と呼ばれる数NuをNu≡2dq/dt/(dq/dt)0で定義します。

 

(dq/dt)0=4πaka(T-Ta) (28)からdq/dt=Nua 2πa(T-Ta)(30)と書けます。=0 のときにはNu=2です。

今日はここまでですが,まだ続きます。

なお,もちろんオリジナルではなく当時,参考文献を見たはずなのですが,それをノートに書いてないので今は不明です。

確か1983年当時に勤めていた会社の蔵書だったと記憶してます。

 

PS:新しい試みとしてワードのルビ(ふりがな)を記号としてココログにコピーしてみたのですが取り合えずのところは失敗でした。おかげで文章を校正する時間を浪費しました。

  

貧乏ひまなしで師走は日曜日でも時間がありません。バスに遅れるので続きはあとで。

 

PS:民主党だけドタバタしてるように見えるのは与党だから仕方ないとはいえ,こんなとき野党は何のためにあるのか?責任追及だけやるしかないのか?そもそも代議士とは何をするのが仕事なのか?

 

政界再編はいずれ必要なんだろうけど,今やられてもそれが終わる頃には私自身の寿命などは終わってるだろうなあ。。。

 

イヤ,またまたマスコミや評論家の多くと同じく,自分のこと抜きで他人を揶揄しちまったか。

 

なんとか法とかの法律に賛成反対とか,国や公共団体の迷惑行為,大きなお世話はやめて欲しいけど,基本的に政治のお世話には期待してないアナーキーを自認している自分。。

 

自分の世話さえままならず,社会的にも無力で有象無象の自分がいまさら何の能書きをたれてるのか。お笑いだ。。。

 

PS2:機械文明や都会のコンビニエンスに依存し切っている身はいまさらどうしようもないという気分もあるけれど,もしもここ10年くらいの間,PCやネットがなかったら,既に廃人になってたかも。。 

 

今のような状態でも,自分が生き生きと生きてるように感じることが多々あるのは悲しいことかも知れないけど,血の通ってないPCやネット依存症に身を委ねて合理化しているおかげでしょうネ。

| | コメント (0) | トラックバック (0)

2010年12月 5日 (日)

プルーム上昇のモデル式(3)

 プルーム上昇式の残りです。このシリーズはこれで終わりです。 

§5.逆転層貫通の条件

プルームが逆転層(inversion layer)を貫通する条件を考察します。

(補足):逆転層とは通常大気と逆に高度zが上昇するにつれて気温が増加していくような層のことです。

 

 断熱冷却されながら上昇していく周辺大気より暖かいプルームが上昇しても,こうした層があると自身よりも暖かい周辺大気に遭遇してプルームが浮力を失なうために,反射されて層を貫通できない可能性があります。

こうした事情は,2006年9/25の記事「ベナール対流の安定性とレイリー数」の序論では次のように書いています。

層の底z=0 では温度がT=T0 天井でT=T1 として定常熱伝導の方程式∇2T= 0 の解はT(z)=T0-βz=T0{(T1-T0)/H}z です。

ここでβ>0 ,つまりT0>T1としないと地球上の重力場の中では対流が起きません。そればかりでなく対流が起きるためにはβの値が断熱温度減率(湿度によって変化)よりも大きくないと対流は起きません。(例えばランダウ「流体力学1」を参照)

 

(補足終わり)※

 

さて,本題に戻ります。

 

 z=ziに微小厚さΔziの鋭い逆転層があって,大気はこの高さz=zi(と厚さΔziの極く薄い層)を除けば中立,つまり層の上下では共にs=(g/Ta)(∂θa/∂z)=0 であるとします。

そして逆転層より下:z<ziの中立状態の温位はθa=θであり,逆転層より上のz>ziの同じく中立状態の温位はθa=θ+Δθiであるとします。

(ⅰ)鉛直浮力プルームの場合

z=ziの層を通過できたと仮定したときにθp≧θa,つまりθp≧θ+Δθiが満足される,すなわちθ'≡θp-θと置いたときにθ'≧Δθiが満足されるならプルームは逆転層を貫通できるはずです。

何故なら,プルームが逆転層を貫通する仮想変位を行ったとき,受ける浮力は上向き,つまりg(θ'-Δθi)/Ta≧0 なので,それは下には押し戻されないからです。

(ただし,θ'≧Δθiは十分条件であって必要条件ではありません。つまりθ'<Δθiで浮力が負の向きであっても上向きの運動量が十分であれば貫通します。)

 微小厚さ:Δziの逆転層に対して,逆転層の強さ:biを,(32)bi≡gΔθia=∫Δzi(g/θa)(∂θa/∂z)dz=∫Δzisdzによって定義します。

また,改めてθ'≡θp-θに対する浮力流束:Fzと容積流束:Vを定義し直しておきます。Fz≡∫P(gρpθ'w'/Ta)dxdy/(πρa),V=∫Pρpw'dxdy/(πρa)です。

 こう定義すると,プルーム断面上でのθ'の平均値は,<θ'>=(Ta/g)(Fz/V)で与えられます。

 

 すぐ上で述べた逆転層貫通の(十分)条件θ'≧Δθiはプルーム断面について,<θ'>≧Δθi,つまり(Ta/g)(Fz/V)≧Δθi⇔ (Fz/V)≧(gΔθi/Ta)と解釈されます。

一方,bi=gΔθiaですから,θa~Taより結局逆転層の貫通条件は(Fz/V)≧biと表わせることがわかります。

ところで,中立大気ではs=(g/Ta)(∂θa/∂z)=0 ですから方程式の1つ:(22)dFz/dt=-s(wV)からFz=F (一定)です。

  

そこで,z=zi~zi+ΔziでもFzi=F(一定)と考えられます。

すると,(23)d(wV)/dt=Fzは,d(wV)/dt=F(一定)となりますから,これもすぐに解けてwV=Ft+Fmを得ます。

それ故,最後の(24)dV/dt=2α(wV)3/2/Vはd(V2/2)/dt=2α(Ft+Fm)3/2となります。

 

t=0でV=0 の点源なら,解はV2/2=(4α/5)F-1{(Ft+Fm)5/2-Fm5/2},すなわちV=(8α/5)1/2-1/2{(Ft+Fm)5/2-Fm5/2}1/2と書けます。

 今の場合は浮力プルーム(Fm=0)を想定しているので,wV=Ft+Fm=Ftであり,V=(8α/5)1/23/45/4です。

 

 したがってw=Ft/V=(8α/5)-1/21/4-1/4を得ます。

 これが,t=tiでz=ziに到達するとすれば,zi=∫0tiwdt=(8α/5)-1/2(4/3)F1/4i3/4ですから,ti=(8α/5)2/3(3/4)4/3-1/3i4/3です。

 

 そこで,Vi=(8α/5)1/2-1/4i5/4=(8α/5)4/3(3/4)5/31/3i5/3を得ますから,結局Vi=(4・35α4/54)1/31/3i5/3です。

 そして,実験によればVi~ 0.0717F1/3i5/3 ⇔ α~ 0.125です。

 以上から,浮力プルームの貫通条件:(Fzi/Vi)≧biは,{54/(4・35α4)}1/32/3i-5/3≧biと書けることがわかります。

したがって,無風時(calm)に中立状態で浮力プルームが逆転層(z=zi)を貫通する条件として,(32)zi≦(4・35α4/54)-1/52/5i-3/5が得られました。

(ⅰ)折れ曲がりプルームの場合

 この場合も,浮力プルームについては(Fzi/Vi)≧biを逆転層貫通の条件であるとします。

  

 周辺大気は中立(s=0)なのでFzi=F(一定)です。

そして,前記事で書いたように,(21)dV/dz=2β(uV)1/2のz=0(t=0)でV=0 の解は,V=u(βz)2ですからVi=u(βzi)2より,Fzi/Vi={F/(uβ2)}zi-2です。

 そこで,折れ曲がり浮力プルームの逆転層貫通条件は,F/(β2u)}zi-2≧biと表現できます。

したがって,有風時に中立状態で浮力プルームが逆転層(z=zi)を貫通する条件として(33)F≧β2ubii2~ 0.25ubii2, or (33)'zi≦β-1{F/(ubi)}1/2~ 2.0{F/(ubi)}1/2が得られました。

(ⅲ)鉛直プルームで浮力がゼロのジェット・プルーム,および浮力流束は不足だが十分な運動量の鉛直流束を持つ場合

 初期浮力がゼロのジェット(F=0)の場合を考えると,s=0 なら(22)dFz/dt=-s(wV)=0 よりFz=F=0 です。

 

 そこで,z<ziではFz=0 (一定)ですがz=ziで浮力流束:-bii=-(∫Δzisdz)Viを獲得して,z>ziではFz=-bii(一定)と考えることができます。

 逆転層:z=ziより上ではプルーム内部の温位θpが外気のそれθ+Δθiよりも小さく,エントレインメントが無視できてV=Vi(一定)と仮定すれば次のように書けるはずです。

z<zi(t<ti)ではFz=0 (一定),wV=Fm(一定),dV/dt==2α(wV)3/2/Vであり,一方,z>zi(t>ti)ではFz=-bii(一定),wV=-bii(t-ti)+Fm,V=Vi(一定)です。

そこで,t>tiではw=-bi(t-ti)+Fm/Viより,w=0 となるのはtmax=ti+Fm/(bii)のときです。

 

このときの高さ:z=zmaxは,zmax=zi+∫titmaxwdt=zi-(bi/2)(tmax-ti)2+(Fm/V)(tmax-ti)で与えられます。よってzmax-zi=Fm2/(2bii2)です。

さらに,ViをFm,ziで表わすためt<tiでのVを求めます。

 

従う方程式はdV/dt=2α(wV)3/2/V or d(V2/2)/dt=2α(wV)3/2=2αFm3/2です。

 

t=0 でV=V0≡Fm/w0の解は,V2/2=2αFm3/2t+(Fm/w0)2, or V={αFm3/2t+(Fm/w0)2}1/2で与えられます。

そして,w=Fm/V=よりzi=∫0ti(Fm/V)dt=2Fm(Vi-V0)/(4αFm3/2)ですからVi=2αFm1/2i+V0です。

点源ならV00 なので,Vi=2αFm1/2iより,zmax-zi=Fm2/(2bii2)=(1/2)(Fm/bi)(2α)-2i-2です。

 

よって,(34)zmax/zi=1+(1/8)α-2(Fm/bi)zi-3を得ます。

このケースでは,zmax3ziを貫通条件として採用します。そして,そして,α~ 0.125としてその条件を求めることにします。

 

無風時(calm)に中立状態でジェット・プルームが逆転層(z=zi)を貫通する条件として,(35)Fm≧16α2ii-3 ~ 0.25bii-3, or (35)'zi≦(4α)2/3(Fm/bi)1/3~ 1.59(Fm/bi)1/3が得られました。

もしも,わずかの浮力Fを持つプルームならz<zi(t<ti)ではFz=F(一定),wV=Ft+Fm,d(V2/2)/dt=2α(Ft+Fm)3/2であり,一方,z>zi(t>ti)ではFz=F-bii(一定),wV=-bii(t-ti)+Ft+Fm,V=Vi(一定)です。

t>tiでは,w=-bi(t-ti)+Ft/Vi+Fm/Viより,w=0 となるのはt=tmax=(biii+Fm)/(bii-F)のときです。

 

このときの高さ:z=zmaxは,zmax-zi=∫titmaxwdt=(Fti+Fm)2/{2Vi(bii-F)}で与えられます。

ここで,Vi=(8α/5)1/2-1/2{(Fti+Fm)5/2-Fm5/2}1/2,かつ(Fti+Fm)2={5Vi2F/(8α)+Fm5/2}4/5なのでzmax-zi={5Vi2F/(8α)+Fm5/2}4/3/{2Vi(bii-F)}なる式を得ます。

 以下,この場合もzmax≧3ziを貫通条件として中立状態大気中の逆転層(z=zi)を貫通する条件式が得られます。

 何故か,ここまででノートが終わってるので,この項はこれでおしまいにします。

参考文献:G.A.Briggs “Plume Rise",U.S.Atomic Energy Commision Division of Technical Information(1969)

 

PS1:土曜日4日午後に,4~6日に文京区春日の文京シビックセンターで開催されていた「障害者ふれあいの集い」に行ってきました。

 

  

  

 

 11/30に出席した手話講習会での情報によると,豊島区や北区でも同時に開催されているらしいです。

 

 実は私もつたない作品="牛乳などの紙パックから再生紙を作り折り紙などを切り貼りしてパウチしたコースター"を何枚か出品していて,1枚30円也で販売されてたりしていました。

   

   

  

  

  

 

PS2:鈴木宗男さん。がんばれ。

 

 一緒にしたら迷惑かもしれないけど高知のスクールバス運転手の件といい,理不尽だと思いますね。

 

 ちまたでウワサの市川海老蔵氏とそのケンカ相手も世論に負けないでくださいネ。。

 

 有名人といえども酒癖はしょうがない面ありますね。

 

 不特定の乗客を乗せるタクシー運転手と同じく(尤もこちらは遊びじゃなく仕事ですが),常識があろうとなかろうと外での個人的飲酒はキケンと背中合わせなのは.大人なんだから覚悟の上でしょう。

 

 男でも女でも,普段は物静かだったりリッパに見える方でも,酒の席では基本的に無礼講です。

 

 酔っ払うと自慢話やホラ話のテンコ盛り,酒乱,泣き上戸,怒り上戸,裸になるなどなど,性格が変わる(実は本性?)ような方々,いっぱい見てきています。本当は羨ましい限りです。

 

 本気で酒でストレス解消しようとに飲みにきているのなら,その解消法はヒトさまざまです。まあ,そうした愚行のお相手も込みで飲み屋の代金を払ってるのだと思っています。

 

 有名人とはいえ,かつての草薙つよポンのケースと同じく,普段の姿と重ねて批判されるのはカワイソーですね。

  

 私だって他人の一方的なストレス解消を迷惑だと感じて,ケンカは弱いけれど思うだけなら「こいつ,どついたろか?」と考えることもしばしばですが。。。

| | コメント (0) | トラックバック (0)

2010年12月 1日 (水)

プルーム上昇のモデル式(2)

 今日もプルーム上昇式の続きです。

§2.理論の単純化,エントレインメントの仮定

 これまでの議論では,未知数Fz,w,Vに対して2つの方程式しか与えられておらず,問題が閉じていません。

 

 これは,プルームの明確な定義が与えられていないため,プルームの大きさを与える式が無いことに起因すると考えられます。

 そこで,連続の方程式∇(ρpp)=0 を用いてプルームの大きさを規定する,ある方程式を与えます。

 プルーム半径をrとして,改めてプルームの平均進行速度の大きさをvsとします。さらに,エントレインメント(entrainment)速度,つまり,プルーム境界から乱流により流入する外気の平均速さをveとします。

 そして,プルームの行路長さ:sに対してds≡vsdt,またはdt≡ds/vsによってdtを定義しておきます。

 すると,質量の保存式から明らかにΔ(ρpπr2s)=ρa2πrveΔs,または(d/ds)(ρp2sa)=2rveです。こうした関係があることをエントレインメントの仮定といいます。

 特に,鉛直プルーム:s=z,dz=wdtの場合にはV=ρpπr2w/(πρa)=ρp2w/ρaなのでエントレインメントの仮定はdV/dz=2rveを意味します。

 また,(d/ds)(ρp2sa)=2rveよりρaはsに依らないとして(r2a){d(ρps)/ds}+(2rρpsa)(dr/ds)=2rveです。

さらに,dt=ds/vsよりd/ds=(1/vs)(d/dt)を用いればdr/dt=ρaep-{r/(2ρp)}{d(ρps)/ds}を得ます。

ここでTaylorによる自然な仮定:ve∝wを採用し,鉛直プルームに対するエントレインメント係数αを導入してve=αwとします。

また,再びV=ρpπr2w/(πρa)=ρp2w/ρaよりr2=ρaV/(ρpw)~V/wと近似できるという考察から,逆にプルーム半径:rをr≡(V/w)1/2によって定義します。

この定義によれば,V=r2wですが,Δ(ρpπr2s)=ρa2πrveΔs,または(d/ds)(ρp2sa)=2rveの関係は変わらず,ρp~ρaの近似で考えると依然としてエントレインメントの仮定はdV/dz=2rveです。

元々,これも仮定ですから以下ではdV/dz=2rveをエントレインメントの仮定として採用します。この右辺にve=αw,r=(V/w)1/2を代入するとdV/dz=2α(wV)1/2を得ます。

以上から,(抵抗力)を無視すると次のような簡単な基本方程式系が得られました。

(4)dFz/dz=-sV,(12)d(wV)/dz=Fz/w,

および,(13)dV/dz=2α(wV)1/2 です。

一応,未知量V,Fz,wの意味を再び述べておきます。

 

Vは鉛直容積流束,Fzは鉛直方向の浮力(gθ'/T)の流束,wはプルーム面の鉛直速度成分です。

  

プルームの行路長さと混同するおそれのある同じ変数sを用いていますが,(4)式のsは,もちろんs≡(g/Ta)(∂θa/∂z)で定義される安定度パラメータです。

 

そして,Morton,Taylor,Turner(1956)によれば,エントレインメント係数αの具体的な値の候補はα=0.093です。一方,Briggs(1968)によればα=0.075です。また,Ricou,Spalding(1961)によればジェットに対しα=0.080です。

 §3.折れ曲がりプルーム(bent-over plume)

 進行方向が完全に鉛直z方向でs=(一定)のプルーム断面が水平面内にある鉛直プルームとは異なる状態のプルームを折れ曲がりプルームといいます。

折れ曲がりプルームの場合,積分面をプルームと交わる鉛直面として,この面上ではs=(g/T)(∂θa/∂z),およびa=uは一定であると仮定します。

前と同様,∫pθ'p)dydz=-(∂θa/∂z)∫ρpw'dydzですがプルームの外ではθ'=0,w'=0 とすれば,これは(d/dx)∫Pρpθ'(u+u')dydz=-(∂θa/∂z)∫Pρpw'dydzを意味します。

そこで,浮力流束を(5)Fz≡∫P(gρpθ'w'/Ta)dxdy/(πρa)の代わりに(14)Fz≡∫P{gρpθ'(u+u')/Ta}dydz/(πρa)で定義すれば,(15)dFz/dx=-sMと表わすことができます。

ただし,再びsはs=(g/Ta)(∂θa/∂z)でありMはM≡∫Pρpw'dydz/(πρa)で定義される量です。

この場合も排出源がhot-sourceなら,浮力流束Fzの初期値Fは{(単位時間にプルーム積分面を通過する総排出熱量)×(g/Ta)/cp}/(πρa)と考えられるので,やはりF=gQH/(πρapa)で与えられます。

 そして,運動量に対する方程式は,(8)d(V)/dz=(Fz/w)a(dV/dz)+(抵抗力)の証明の際に遭遇した式:(d/dn)(∫PρppnpidS)={1+(dr/dn)2}1/2(∫Cρpppi)+δi3P(gρpθ'/Ta)dSにおいて進行方向に置くと得られます。

 d/dnをd/dxに変え,p=(u+u')i+v'+w',d={1+(dr/dx)2}-1/2(-dz+dy)+{1+(dr/dx)2}-1/2(dr/dx)dlとして,必要なi=3の場合のみを書き下します。

 

 それは,(d/dx){∫Pρp(u+u')w'dS}=∫P(gρpθ'/Ta)dS+{1+(dr/dx)2}1/2(∫Cρppw'd)=∫P(gρpθ'/Ta)dS-Cρpv'w'dz+∫Cρpw'2dy+∫Cρp(u+u')w'(dr/dx)dlです。

 u'<<uなのでu+u'~ uと近似して有効運動量流束をuMeff;Meff≡{∫Pρpw'dS-∫dx∫Cρpw'(dr/dx)dl}/(πρa)と定義しu',v',w'の2次以上の項を無視します。

 

 すると,Fz={u∫P(gρpθ'/Ta)dS}/(πρa)より,結局運動方程式として(16)u(dMeff/dx)=Fz/uが得られます。

 鉛直運動量の水平方向流束uMeffの初期値はプルーム面を通過する鉛直運動量ρpwπr2/(πρa)=ρpwr2aの総流束(鉛直方向流束)ですから,その初期値はFm=(ρ0002a)×w0=ρ00202a=ρ000aです。

そして,Briggs(1972),Richards(1963)によれば有効運動量流束uMeffの運動量流束uM;M=∫Pρpw'dydz/(πρa)に対する比率の評価値はMeff/M ~ 2.3です。

zをプルーム軸心の高さ,dt≡dx/uとしてw≡dz/dt=udz/dxと定義します。wはプルーム軸心の上昇速度を与え,解(x,z),or z=z(x)はプルームの軌跡を示します。

ここで,折れ曲がりプルームにおけるエントレインメント係数をβとします。すなわちve=βwです。有効プルーム半径をrとすると,Meff=πr2ρpw/(πρa)=wρp2aです。

折れ曲がりプルームの場合のエントレインメントの仮定はd(ρpπr2u)/dx=2πrρae,つまりd(ρp2u/ρa)/dx=2rveと表わせます。

そこで,ud(Meff/w)/dx=2rve=2rβwです。

 

r≡(Meff/w)1/2と定義すれば,(17)ud(Meff/w)/dx=2β(wMeff)1/2を得ます。

さらにwV≡uMeffによってVを定義すれば,Vは水平方向の容積流束でありr=(V/u)1/2です。これとu/dx=w/dzなる関係を用いて基本方程式系:(15),(16),(17)を書き直せば次のようになります。

すなわち,(18)dFz/dz=-s'V,(19)s'=s(M/Meff),(20)d(wV)/dz=Fz/w,(21)dV/dz=2β(uV)1/2です。

(21):dV/dz=2β(uV)1/2は,z=0 でV=0 という点源(point source)の初期条件に対しては容易に解けて,(V/u)1/2=βz,つまりr=βzです。

Richards(1963),TVA-Report(1967)によれば折れ曲がりプルームのエントレインメント係数はβ~ 0.5です。

§4.プルーム上昇高度の算出

 具体的にプルーム上昇式を求めるに際して,さらに次のような2つの仮定を与えます。

(仮定Ⅰ):鉛直プルームを満足する状態を無風状態(静穏:calm)に対応させ,折れ曲がりプルームを有風時の状態に対応させる。

(仮定Ⅱ):初期時には鉛直運動量か浮力流束のいずれか一方が支配的である。前者はジェットと呼んでF=0 とし,後者を浮力プルームと呼んでFm=0 とする。

 さて,まずdx=udt,dz=wdtにより,全ての方程式を独立変数tの式に変換します。

鉛直プルームについては,(4),(12),(13)より

 

(22)dFz/dt=-s(wV),(23)d(wV)/dt=Fz,(24)dV/dt=2α(wV)3/2/Vです。

また,折れ曲がりプルームについては,(18),(20),(21)より

 

(25)dFz/dt=-s'(wV),(26)d(wV)/dt=Fz,(27)dV/dt=2βwV(u/V)1/2です。

 以下,順に安定大気:s>0,s'>0 (s,s'は一定値を仮定)の中でのプルーム上昇高さを求めます。

(ⅰ)鉛直プルーム(無風:calm)

 このケースの基本方程式系は(22)dFz/dt=-s(wV),(23)d(wV)/dt=Fz,(24)dV/dt=2α(wV)3/2/Vです。

(22),(23)よりd2(wV)/dt2=-s(wV)を得ます。これはよく知られている線形振動の方程式です。

 

そこで,ω≡s1/2としてこれをd2(wV)/dt2=-ω2(wV)と書き,初期条件:t=0 でwV=Fm,d(wV)/dt=Fz=Fの解を求めると,wV=Fmcosωt+(F/ω)sinωtを得ます。

鉛直運動量流束wVがゼロとなる最初のt=t1がプルーム上昇高さを与えると考えられます。

  

上昇高さはΔh=∫0t1wdtですから,仮定Ⅱの浮力プルーム(Fm=0)ではt1=π/ω,よりΔh=ω-10πwdλです。

 

同様にジェットプルーム(F=0)ではt1=π/(2ω)よりΔh=ω-10π/2wdλ(λ=ωt)です。

     浮力プルーム(Fm=0)の場合

 (V2/2)/dt=V(dV/dt)=2α(wV)3/2の右辺にwV=(F/ω)sinωtを代入してd(V2/2)/dt=2α(F/ω)3/2sin3/2ωtです。

故に,V2/2=2α(F/ω)3/2ω-1{∫ωtsin3/2λdλ+const.},すなわちV=2α1/23/4ω-5/4{∫ωtsin3/2νdν+const.}1/2です。

 w=wV/V=(F/ω)sinωt/Vより,Δh=ω-10πwdλ=(1/2)α-1/21/4ω-3/4[∫0πsinλdλ{∫0λsin3/2νdν+const.}-1/2]を得ます。

したがって,ω=s1/2より定数cをc≡(1/2)α-1/2[∫0πsinλdλ{∫0λsin3/2νdν+const.}-1/2]で定義すれば,Δh=cF1/4-3/8なる関係式が成立することになります。

 t=0 のとき,V=0 の点源(const.=0)ではMorton,Taylor,Turner(1956)の数値計算があり,その示唆するところによればα=0.132としてc~ 5.0ですから,浮力プルーム上昇高さの無風状態での評価式として(28)Δh=5.0F1/4-3/8を得ます。

     ジェットプルーム(F=0)

 wV=Fmcosωtよりd(V2/2)/dt=V(dV/dt)=2α(wV)3/2の右辺にwV=(F/ω)sinωtを代入すれば,d(V2/2)/dt=2αFm3/2cos3/2ωtです。

故に,V2/2=2αFm3/2ω-1{∫ωtcos3/2λdλ+const.},すなわちV=2α1/2m3/4ω-1/2{∫ωtcos3/2νdν+const.}1/2です。

 そして,w=wV/V=Fmcosωt/Vより,Δh=ω-10π/2wdλ=(1/2)α-1/2m1/4ω-1・2[∫0π/2cosλdλ{∫0λcos3/2νdν+const.}-1/2]です。

したがって,ω=s1/2より定数c0をc0≡(1/2)α-1/2[∫0π/2cosλdλ{∫0λcos3/2νdν+const.}-1/2]と定義すれば,Δh=c0m1/4-1/4なる関係の成立することがわかります。

 t=0 のとき,V=0 の点源(const.=0)では,μ=π/2-νと変数変換すれば,cosν=sinμより∫0λcos3/2νdν=-∫π/2π/2-λsin3/2μdμ,c0=(1/2)α-1/2[∫0π/2cosλdλ{∫π/2-λπ/2sin3/2μdμ}-1/2]=(1/2)α-1/2[∫0π/2sinσdσ{∫0σsin3/2μdμ}-1/2]となります。

 点源のとき,c=(1/2)α-1/2[∫0πsinλdλ{∫0λsin3/2νdν+const.}-1/2]に対してc0=(1/2)α-1/2[∫0π/2sinλdλ{∫0λsin3/2νdν+const.}-1/2]であり,ほとんどcと同じです。

cを求めたのと同じ数値計算によれば,α=0.132に対してc0~ 4.0ですから,無風状態でのジェットプルームの評価式として(29)Δh=4.0(Fm/s)1/4を得ます。

(ⅱ)折れ曲がりプルーム(有風時)

 このケースの基本方程式系は(25)dFz/dt=-s'(wV),(26)d(wV)/dt=Fz,(27)dV/dt=2βwV(u/V)1/2です。

(25),(26)よりd2(wV)/dt2=-s(wV)なのでω'=s'1/2としてwV=Fmcosω't+(F/ω')sinω'tです。

 

これを(27)から得られるd(2V3/2/3)/dt=2βu1/2wVの右辺に代入すると,d(2V3/2/3)/dt=2βu1/2{Fmcosω't+(F/ω')sinω't}です。

 故に,V=[3βu1/2{(Fm/ω')sinω't+(F/ω'2)(1-cosω't)}]2/3=[(3βu1/2/s'){(ω'Fmsinω't+F(1-cosω't))}2/3です。

 一方,(21)dV/dz=2β(uV)1/2のz=0 (t=0)でV=0 の解はV=u(βz)2ですから,上式のVと等置するとz=(V/u)1/2/β={3/(β2us')}1/3{(ω'Fmsinω't+F(1-cosω't))}1/3を得ます。

 三角関数の合成ω'Fmsinω't-Fcosω't=(ω'2m2+F2)1/2sin(ω't+δ)を用いると,zの最大値がzmax={3/(β2us')}1/3{F+(ω'2m2+F2)1/2}1/3で与えられることがわかります。

 そして,(19)s'=s(M/Meff)によりs'∝ sなので,有風時の浮力プルーム(Fm=0)の場合は,Δh=c1{F/(us)}1/3,ジェットプルーム(F=0)の場合は,Δh=c2(Fm/u)1/3-1/6と書けます。

 ただし,c1=[6/{β2(s'/s)}]1/3,c2=[3/{β2(s'/s)}]1/3ですから,例えばβ~ 0.5,(s'/s)~ 2.3を代入して計算するとc1~ 2.19,c2~ 1.73です。

一方,Briggsによれば,有風時の浮力プルームの場合の上昇式は,(30)Δh=2.4{F/(us)}1/3,有風時のジェットプルームの場合の上昇式は(31)Δh=1.5(Fm/u)1/3-1/6で与えられるらしいです。

今日はここまでです。この項目はまだ少し続きます。

参考文献:G.A.Briggs "Plume Rise",U.S.Atomic Energy Commision Division of Technical Information(1969)

| | コメント (0) | トラックバック (0)

2010年11月28日 (日)

プルーム上昇のモデル式(1)

 唐突ですが,過去の資料を整理していたら1983年(33歳当時)の仕事関係の覚書きノートを見つけたので,ブログに覚書きを残す一環としてこれも書いておこうと思います。

内容は,古典的なBriggsによるプルーム(plume)上昇式の紹介です。(←ノートにはプリューム上昇式と書いてありましたが,今はplumeはプルームと読むべきだと思っているので直しました。)

当時は,環境アセス会社の社員で丁度,冷却塔(cooling-tower)からの温排気のシミュレーションについて調べていた頃と思います。

 

※(この後,最終的にはANLモデル(Argonne国立研究所モデル)と格闘することになったのですが。。この会社をやめて10年以上後にはSACTIというプログラムとも出会いましたね。illinois大学関連かな?) 

ただし,当時は常識だった気象力学の知識は今ではかなり忘れているので,この記事内容だけで私自身がわかるように補足しています。

.プルーム上昇式の導出

§1.基礎理論(鉛直プルーム)

プルーム内の局所密度ρpは不変:∂ρp/∂t=0 とします。

 

すると質量保存の連続の方程式は(1)∇(ρpp)=0 となります。ただしp,およびρpはプルームガスの局所速度,および密度です。

また,局所プルーム要素の運動は断熱的とします。

 

断熱条件はプルームの温位をθpとして(2)dθp/dt=0 です。

 

温度がT,圧力がpの気体に対して温位はθ≡T(p00/p)(R/Cp)で定義されます。ただしRは気体定数:R≒8.31J/(molK),Cvはモル定積比熱,Cpはモル定圧比熱です。

 

vとCpの間にはMayerの関係式:Cp=Cv+Rが成立します。

dθp/dt=0 は中立大気中の運動では浮力が保存すること,移流以外には熱の放射吸収がないことを意味します。

(補足):ガス(気体)が希薄で理想気体近似ができるとすると,その状態方程式は気体分子の分子量をmとしてp=ρRT/mです。p,ρ,Tは,それぞれ気体の圧力,密度,絶対温度です。

断熱変化の場合には,c1を定数としてPoissonの式:p=c1ργが成立します。γは比熱比でγ≡Cp/Cvです。

p=ρRT/mの両辺の対数微分を取ると,dlnp=dinρ+dlnTですからdlnT=dlnp-dlnρですが,断熱変化ではp=c1ργよりlnp=γlnρ+lnc1が成立するため,dlnρ=dlnp/γです。

それ故,dlnT=(1-1/γ)dlnp=(R/Cp)dlnp,すなわちdln{T/p(R/Cp)}=0 ,またはT/p(R/Cp)=(一定)です。

そこで,ある定数c2を与えθ≡c2T/p(R/Cp)によって温位という量θを定義すると,断熱変化はdθ=0 を意味するので,断熱変化の場合の温位の時間変化率はdθ/dt=0 です。

特に,圧力pの値がp00のときθ=Tとなるように定数c2を定めれば,θ=T(p00/p)(R/Cp)という表現ができます。

さて,プルームの周辺大気について対流が生じない条件はdp/dz=-ρg(静力学平衡)ですが,これは対数表現ではdlnp/dz=-mg/(RT)です。ただし,鉛直上向きをz軸の正の向きとしています。

さらに,断熱変化:dθ=0 ならdlnp=(Cp/R)dlnTなので,熱の放射吸収のない断熱環境で対流が生じない条件はdlnT/dz=-mg/(CpT),またはdT/dz=-mg/Cpと書けます。

そこで,鉛直高さzと共にこの比率で温度が減衰すれば対流が生じず,しかも高さ方向の変化dzに対する温度変化は断熱変化なので,温位勾配はゼロ:dθ/dz=0 です。

 

この勾配:dT/dz=-mg/Cpの温度減率(気温減率)を断熱減率といい,この状態にある大気を中立状態といいます。

ところで,通常の地上付近成層圏の大気の体積組成(分子数組成)は分子量28の窒素N2が78.08%,分子量32の酸素O2が20.95%,その他はアルゴン,二酸化炭素,水蒸気などです。

そこで,地上付近の乾燥大気の分子量:m=md28.964です。(小倉義光「一般気象学」より)

また,N2,O2のような2原子分子では理論的にはモル定積比熱は(5/2)R,モル定圧比熱はCp=(7/2)Rです。

そこで,この大気の断熱減率を計算するとΓd≡mdg/Cp=28.964(g/mol)×9.8(m/s2)×2/{7×8.31J/(molK)}=9.8(K/km)です。

 

これは,100m上昇する毎に気温が約1度下がるような状態です。この気温減率:Γd=mdg/Cp=9.8(K/km)を乾燥断熱減率といいます。

一般には,大気は湿度ゼロの乾燥空気ではなく水分を含んでいます。

 

容積Vの中に質量Mdの乾燥空気と質量Mvの水蒸気(vapor)が混在しているとき,この空気の密度はρ=(Md+Mv)/V=ρd+ρvです。

そして,pdを乾燥大気の分圧,eを水蒸気圧とすると,乾燥大気,水蒸気の状態方程式はそれぞれpd=ρdRT/md,e=ρvRT/mvです。

 

全大気圧pはp=pd+e=ρdRT/md+ρvRT/mvですから,p=ρRT{1+(mv/md)(ρvd)}/{1+(ρvd)},またはp=ρRT*/md;T*≡T{1+(mv/md)(ρvd)}/{1+(ρvd)}と書けます。

ここで,T*は仮温度(virtual temperature)と呼ばれる量で具体的にmd=28.964(g/mol),mv=18(g/mol)を代入すると,水蒸気の混合比:q≡ρvd=(md/mv){e/(p-e)}に対してT*≒T(1+1.61q)/(1+q)=T(1+0.61r)と表現されます。

ただし,r≡q/(1+q)=ρv/(ρv+ρd)は比湿(specific humidity)と呼ばれる量で,式から明らかにr≦1です。

 

乾燥大気の温度Tに対する安定条件:dT/dz=-mdg/Cp=-Γdは,Tを仮温度T*に置き換えると近似的にそのまま成立すると考えれば,dT*/dz=-mdg/Cp=-Γdなる式が安定条件に対応します。

 

左辺のT*にT*=T(1+0.61r)を代入すると,(1+0.61r)(dT/dz)=-Γdです。そこで,不飽和で凝結のない大気の断熱減率ΓはΓ=Γd/(1+0.61r)で与えられ,Γ≦Γdです。

さらに,飽和(saturate)して凝結水滴があり潜熱(latent heat)の存在によって温度降下が鈍くなる場合の湿潤断熱減率をΓwとします。

 

飽和水蒸気圧をesと書き,単位質量当たりの液体の水が水蒸気になるのに必要な潜熱(蒸発熱)をLとすれば,esと温度Tの関係は理論的にはClausius-Clapeyronの式:des/dT=Les/(RT2)で表わされます。

 

(上式は,lnes=-L/(RT)+const.またはes=Aexp{-L/(RT)}を意味しますが,2006年7/30の記事「気液平衡の統計力学」ではこれより精密な修正式:lnes=lnT-L/(RT)-1+const.またはes=ATexp{-L/(RT)}を与えています。)

  

例えばT=100℃ではL=597(cal/g)ですが常温付近でのLはTの変動に鈍感なのでLがTに無関係と近似してこれを積分すれば,ln(es/es0)=(mvL/R)(1/T0-1/T)です。

  

これはesがTだけに依存して周囲の大気圧pには依らない形の関係式です。

 

さて,乾燥空気と混合比qの水蒸気が混合した空気塊を考えます。

 

まず,単一気体の断熱変化は,0=dlnθ=(Cp/R)dlnT-dlnpで表わされます。これはエネルギー収支を示す式と見るなら,0=RTdθ/θ=CpdT-RTdp/pと表現されます。

 

そこで,乾燥空気と水蒸気が混合しただけの大気において水蒸気が飽和に達していないとき,各成分の上記左辺の温位変化:RTdθ/θへの寄与は,それぞれCpdT-RT(dpd/pd)=0,qCpvdT-qRT(de/e)=0 です。

   

ここで,qはこの大気の乾燥空気に対する水蒸気の混合比でq=ρvd=(md/mv)(e/pd)=(md/mv){e/(p-e)}です。また,Cpvは水蒸気の定圧モル比熱です。このとき,もちろんpd=p-eです。

   

この空気塊が上昇してq=qsとなって飽和した後,さらに上昇してdqsの水蒸気が凝結したとすれば,空気塊はLdqsだけの潜熱を受けるので,変化が断熱なら水蒸気が飽和した空気塊のエネルギー収支の式は次のようになります。

  

すなわち,CpdT+qspvdT-RT{d(p-es)/(p-es)+qsdes/es}+Ldqs=0 となります。qsは飽和時の混合比です。

  

これに,des/dT=Les/(RT2)を代入すれば,-Ld(qs/T)=CpdT/T+qspvdT/T-Rd(p-es)/(p-es)を得ます。

  

d=ρdRT/md,es=ρvRT/mvですから,静力学平衡:dp/dz=-ρg, or d(pd+es)/dz=-(ρd+ρV)gは,d(pd+es)/dz=-g(mdd+mvs)/(RT),またはdp=-g{md(p-es)+mVs}/(RT)dzと表現できます。

 

そこで,d(p-es)=-g{md(p-es)+mVs}/(RT)dz-desを-Ld(qs/T)=(Cp+qspv)dT/T-Rd(p-es)/(p-es)に代入すると,-Ld(qs/T)=(Cp+qspv)dT/T+g{md+mVs/(p-es)}dz/T+Rdes/(p-es)です。

  

これから,dT/dz=-g{md+mVs(p-es)}/(Cp+qspv)-RT(des/dz)/{(p-es)(Cp+qspv)}+{LT/(Cp+qspv)}{d(qs/T)/dz}を得ます。

   

ところが,qs=(md/mv){es/(p-es)}よりes/(p-es)=mvs/mdです。

   

そこで,dqs/dz=(md/mv)[(des/dz)/(p-es)-{d(p-es)/dz}{es/(p-es)2}]=(md/mv)[(des/dz)/(p-es)-{qs/(p-es)}{d(p-es)/dz]と書けます。

    

つまり,(des/dz)/(p-es)=(mv/md)(dqs/dz)+{qs/(p-es)}{d(p-es)/dz}です。これらの式を用いて,できるだけesを消去します。

   

以上から,湿潤断熱減率は,Γw≡-dT/dz=Γd/(1+qspv/Cp)+mV2g/{mdp(1+qspv/Cp)}+mvRT(dqs/dz)/{mdp(1+qspv/Cp}+RT{qs/(p-es)}{d(p-es)/dz}/{Cp(1+qspv/Cp-LT{d(qs/T)/dz}/{Cp(1+qspv/Cp)}です。

 

この式は閉じていないので古い正野重方著の岩波全書の「気象力学」などには近似公式もあるようですがここでは主要課題ではないので割愛します。ただ,一般に湿潤断熱減率は一定値ではありません。

 

(なお,水蒸気の定圧比熱については,2008年1/6の記事「氷,水,水蒸気の比熱」,および2009年2/9の記事「水蒸気の比熱」などを参照してください。) (補足終わり) ※

さて,浮力以外の気圧傾度力は無視できると考えて外気圧pは近似的に静力学平衡:dp/dz=-ρgを満足しているとします。

 

また,大気を流体と見るとき,対象とする運動でのReynolds数:Reは非常に大きいため,1/Reに比例する分子粘性項も無視します。

そして,プルームガスの周辺大気(ambient air)の局所速度,圧力,密度,温度,温位を,それぞれa,paa,Taaとします。

特に,気圧については,∂pa/∂x=∂pa/∂y=0,∂pa/∂z=-ρagが成立しており,プルームガス圧については近似的にpp=paが成立するとします。

この場合,周辺大気の運動方程式は,ρaa/dt=-∇pa-ρa=0 です。ただし,,は,それぞれx,y,z軸の方向単位ベクトルを表わします。

この方程式の解は明らかにaconstで与えられます。

一方,プルームガスの満たす運動方程式はρpp/dt=-∇pp-ρp~(-∇pa-ρa)-(ρp-ρa)gです。

 

すなわち,この近似での運動方程式はρp(dp/dt)=(ρa-ρp)/g,またはdp/dt={(ρa-ρp)/ρp}gです。これらの右辺が浮力項(buoyancy term)です。

状態方程式:pa=ρaRTa/M,pp=ρpRTp/Mから,pp~paの近似ではρap ~Tp/Ta~θpaなので温位偏差θ'をθ'≡θp-θaと定義すれば,(ρa-ρp)/ρp~ (θp-θa)/θa=θ'/θaです。

周辺大気の温位θaの定義a≡Ta(p00/pa)(R/Cp)において特にp00~paなるパラメータを採用すればθa~Taなので,結局この近似(Boussinesq近似)での運動方程式は,(3)dp/dt=(gθ'/Ta)となります。

連続方程式が質量の保存を意味するのに対して,運動方程式は運動量の保存を意味します。

さらに,pa'と書き,特に速度偏差:'の鉛直z方向の成分(or第3成分)をw'とします。また,周辺大気の風速aは鉛直成分を持たないとして,その向きをx軸(or1軸)の正の向きに取ります。

すると,p≡u,vp3=v'3=w'と書けます。さらにθaa,uはzだけの関数でありθa=θa(z),ρa=ρa(z),u=u(z)のように表わせるとします。

そして水平プルーム断面上ではw'と浮力gθ'/Taが比例する:gθ'/Ta∝ w'と仮定します。比例係数をc(z)とするとgθ'/Ta=c(z)w'です。

定常状態が実現されていると仮定し,面積分することによってプルームを水平面で切った断面上での平均量に対する方程式を求めます。

ただし,プルーム断面の外では明らかにθ'=0,かつ'=0 で特にw'=0 です。

この結果は,まず浮力の効果についての方程式:(4)dFz/dz=-sVを与えます。ただし,Fzは浮力流束で,(5)Fz≡∫P(gρpθ'w'/Ta)dxdy/(πρa)で定義されます。

Vは鉛直容積流束で,(6)V≡∫Pρpw'dxdy/(πρa)です。

 

また,sは安定度パラメータと呼ばれる量で,(7)s≡(g/Ta)(∂θa/∂z)です。

(証明):(2)dθp/dt=0 より,∂θp/∂t+p∇θp=0ですから,定常状態:∂θp/∂t=0 ではρpp∇(θa+θ')= 0です。

そして,∂θa/∂x=∂θa/∂y=∂θ'/∂x=∂θ'/∂y=0 よりρpp∇(θa+θ')=-ρpw'(∂θa/∂z)です。

 

これに,(1)∇(ρpp)=0 にθ'を掛けたものを加えると∇(ρpθ'p)=-ρpw'(∂θa/∂z)を得ます。

これを全水平面で積分すると∫pθ'p)dxdy=-(∂θa/∂z)∫ρpw'dxdyですが,無限遠では明らかにθ'=0 で,プルーム断面を含む全水平面領域では全ての量はzに依らず一定ですから,左辺の(∂/∂z)も積分記号の外に出せます。

プルーム水平断面の外ではρpθ'w'=ρpw'=0 なので,上式は(d/dz)∫Pρpθ'w'dxdy=-(∂θa/∂z)∫Pρpw'dxdyとなります。

周辺大気の密度ρaは"上昇高度までの層当たり平均密度=(一定)"と考え,両辺に(g/Ta)/(πρa)を掛けると{∂(1/Ta)}{(d/dz)∫Pρpgθ'w'dxdy}/(πρa)=-(g/Ta)(∂θ/∂z){∫ρpw'dxdy}/(πρa)です。

 

よって,dFz /dz=-sV-(1/Ta)(dTa/dz)Fzを得ました。

ところが,dTa/dz~dθa/dz,Fz ~gV(<θ'>/θa)なので(1/Ta)(dTa/dz)Fz~ (sV)・(<θ'>/θa)<<(sV)ですから,右辺第1項:(sV)に比べて右辺第2項:(1/Ta)(dTa/dz)Fzを無視します。

そこで,結局dFz /dz=-sVを得ます。(証明終わり)

次に,運動量に対する方程式は,(8)d(V)/dz=(Fz/w)a(dV/dz)+(抵抗力)です。ただし,Vは(9)V≡∫Ppρpw'dxdy/(πρa)で定義される運動量流束です。

 

この式で定義されるベクトルをプルームの速度と考えます。wはのz成分です。

(証明):(1)∇(ρpp)=0 より,p・∇(ρpp)=0 ,つまりΣj=13pijppj)=0 です。

 

 また,(3)dp/dt=(gθ'/Ta)より,定常状態ではρp(p∇)p=(gρpθ'/Ta),つまりΣj=13ρppjjpi=(gρpθ'/Tai3です。

したがって,∇(ρpppi)=Σj=13vpi∂jppj)+Σj=13ρppjjpi=(gρpθ'/Tai3です。

 プルームの進行方向をとするとき,プルームの断面と側面で囲まれた領域VにGaussの定理を適用すると,∫V∇(ρpppi)dV=∫P+ρppnpidS-∫P-ρppnpidS-Δn{1+(dr/dn)2}1/2(∫Cρpppi)です。 

ただし,ここでのrはプルーム半径(plume-radius)を意味します。プルーム断面は等方的なので,これは断面の形を円と想定したときの半径です。

 

-Δndはプルーム側面の外向き法線の向きを持ちΔnは側面方向の厚みです。それ故,dは周の微小線素の長さを持つ内向き法線方向ベクトルです。

 故に,Δ∫PρppnpidS≡∫P+ρppnpidS-∫P-ρppnpidS=Δn{1+(dr/dn)2}1/2(∫Cρpppi)+δi3Δn∫P(gρpθ'/Ta)dSです。

 

 よって,(d/dn)(∫PρppnpidS)={1+(dr/dn)2}1/2(∫Cρpppi)+δi3P(gρpθ'/Ta)dSを得ます。

 一方,連続の方程式から(d/dn)(∫PρppndS)={1+(dr/dn)2}1/2(∫Cρppです。

 

 さらにvpi=vai+v'iですから,(d/dn)(∫PρppnpidS)=δi3P(gρpθ'/Ta)dS+vai(d/dn)(∫PρppndS)+{1+(dr/dn)2}1/2(∫Cρppv'i)を得ます。

 特にをz方向のに選べばvpn=w'より,(d/dz)(∫Pρpw'vpidS)=δi3P(gρpθ'/Ta)dS+vi(d/dz)(∫Pρpw'dS)+(抵抗力)です。

 

 ここで,(抵抗力)≡{1+(dr/dz)2}1/2(∫Cρppv'i)であり,d={1+(dr/dz)2}-1/2(-dy+dx)+{1+(dr/dz)2}-1/2(dr/dz)dlです。(dl=(dx2+dy2)1/2)

 一方,∫P(gρpθ'/Ta)dS=(1/w)∫P(g/Ta)(ρpθ'w)dS=(1/w)∫P(g/Ta)(ρpθ'w')dS+(1/w)∫P(g/Ta){ρpθ'(w-w')}dSです。

V≡∫Ppρpw'dS/(πρa),V≡∫Pρpw'dS/(πρa),およびvp3=w'よりw=∫Pρpw'2dS/∫Pρpw'dSですから,w-w'={∫Pρpw'2dS-w'∫Pρpw'dS}/{∫Pρpw'dS}です。

 仮定によりgθ'/Ta=c(z)w'を代入すると,(1/w)∫P(g/Ta){ρpθ'(w-w')}dS={c(z)/w}[{∫Pρpw'dS}{∫Pρpw'2dS}-{∫Pρpw'2dS}{∫Pρpw'dS}]/{∫Pρpw'dS}=0 です。

 したがって,∫P(gρpθ’/Ta)dS=(1/w)∫P(g/Ta)(ρpθ'w')dS故,(d/dz)(∫Pρpw'vpidS)=(1/w)δi3P(g/Ta)(ρpθ'w')dS+vi(d/dz)(∫Pρpw'dS)+(抵抗力),

 

 または,(d/dz)(∫Ppρpw'dS)=(1/w){∫P(g/Ta)(ρpθ'w')dS}a(d/dz)(∫Pρpw'dS)+(抵抗力)なる式を得ます。

 最後の表式は,まさにd(V)/dz=(Fz/w)a(dV/dz)+(抵抗力)そのものです。

 

※特にプルームが軸対称で完全に鉛直向きなら(抵抗力)={1+(dr/dz)2}1/2(∫Cρppv'i)=0でd(wV)/dz=Fz/wです。※(証明終わり)

 さて,dFz/dz=-sV,d(wV)/dz=Fz/wを解くために次のような初期条件を与えます。

まず,初期速度は吐出速度に等しいとします。すなわち,(10a)=w0です。そして初期運動量流束は,(10b)V=(wV)=Fm≡(ρ0a)w0202とします。ρ0は排出ガス密度,r0は煙突口の口径です。

さらに,初期の浮力流束は,(10c)Fz=F≡(1-ρ0a)gw002で与えられるとします。

  

つまり,初期容積流束をV0≡w002とすれば初期の総質量流束はπρ00=πρ0002ですから,定義により運動量流束はFm=πρ000/(πρa)=(ρ0a)w00です。

  

そして,初期浮力流束はFz=F=(1-ρ0a)gV0~{(T0-Ta)/T0}gV0となるわけです。

 

ところで,排出源がhot-sourceの場合の初期浮力流束は,(11)F=gQH/(πρapa)で与えられます。

 

ここで,QHは排出されるガスによって単位時間に排出される熱量(cal/sec)であり,cpは空気の単位質量当たりの熱容量(定圧比熱)(cal/(kg・K))です。

 

なぜなら,Fz=F=(gθ'/Ta)(πρ0002)/(πρa)={g/(πρaa)}(θ'πρ00)ですが,πρ00は総質量流束なので排出ガスによって単位時間に外気に付加される熱量はQH=cpT'πρ00 ~cpθ'πρ00に等しいからです。 

 

つまり,θ'πρ00=QH/cpによってF=gQH/(πρapa)が得られるわけです。

 

この場合,排熱強度QHと排出ガス密度ρ0の情報を得ることにより逆にV0=w002が求まります。

 

そこで,煙突口径r0のデータから吐出速度w0が算定されて,結局全ての初期情報を得ることができます。

  

ただし,後述するように3つの未知数Fz,w,Vに対し微分方程式は2つしかないので解を決めるにはもう1つ方程式が必要です。

今日はここまでですがまだ続きます。

参考文献:G.A.Briggs "Plume Rise",U.S.Atomic Energy Commision Division of Technical Information(1969),小倉義光「一般気象学」(東京大学出版会),小倉義光「気象力学通論」(東京大学出版会)

| | コメント (0) | トラックバック (0)

2010年9月19日 (日)

台風の進路(コリオリの力)(再々掲)

 丁度,台風シーズンです。

 そこで,私自身の息抜き(手抜き?)として,4年前に書いたテキストだけの解説記事:「台風の進路(コリオリの力) 」に図を入れて再掲してみます。

(※再掲記事):2006年8月16日(水)

 そろそろ台風が頻発する季節になりました。

 今日は,北半球では赤道付近で発生し,海域から多量の水分やエネルギーを吸収しながら発達して北上する台風が,なぜ右(東)の方に進路を変えていくのか?そして,なぜ上空から見て左巻き(反時計回り)の風が吹くのかという,ごくありふれた疑問について解説してみます。

 例として,ちょっと古いけど左回転しているLPレコードがあり,その上に"一寸法師"よりも小さい小人が乗っているという仮想的な状況を考えてみます。(左回転は仮定であって実際のLPレコードは裏から見ない限り右回転(時計回り)です。)

 レコードの中心は地球の北極に相当し,レコードの最大半径の部分は地球の赤道に相当します。 

 まず,レコードの回転している"最大半径=赤道"の上にいる小人が"レコードの中心=北極"めがけて真っ直ぐに小石を投げたとします。

 本人は真っ直ぐ中心に向かって投げたつもりでも,小石が手を離れた瞬間には慣性によって小石はレコードの回転スピードと同じ速度で右に向かう接線速度を持ちますから,実はそれは中心の方向に向かって真っ直ぐに飛ぶのではなくて,次第に右の方に逸れていくということになりますね。

 ↑ここで右というのは,小石を投げた小人にとっての右です。(わかっているとは思いますが念のため))

 次に逆に"中心=北極"の上に小人がいて今度は"最大半径=赤道"めがけてやはり小石を投げたとします。

 今度は北極で小人は自転しているかもしれませんが,スピンの回転半径はゼロなので,その慣性による小石の左右方向への速度はゼロですから確かに"真っ直ぐ"進むはずです。

 ところが,レコードの上,つまり北半球の地球上にいる人は"左から右=西から東"に回転しています。その人から見ると"上=北"から真っ直ぐ飛んでくる小石は"下=南"から見て"左に左に(西に西に)"逸れていくように見えます。

 逆に"小石を投げた方=北"から見ると,見かけ上はやはり右の方に逸れていくわけですね。

 小石を台風だとみなし地球の自転の角速度をΩとすると Ωは"360度=2πラジアン(rad)"を24時間で回転する角速度です。

 地球の半径をRとし,緯度をθとすると,そこでの回転半径はRcosθですから,回転の接線速度はΩRcosθです。

 したがって"赤道"での接線速度は最大速度"ΩR=時速1667km=秒速463m"ですが,日本付近の緯度:θ=35度での接線速度は"ΩRcosθ=秒速379m"で,日本付近では回転速度は赤道より"約2割=秒速80m"くらい減少しています。

 したがって,赤道付近で発生した台風は地球のまさつにより"ΩR=時速1667km=秒速463m"の地球にひきずられて慣性による右向きの速度を持っていて,その右向き速度は北上しても全く変わらないものです。

 しかし,地球自身の回転速度は緯度が上がるにつれて次第に小さくなるものですから,日本付近では1秒間に80mくらいの割合で,右(東)へ右へと逸れていくことになります。

 先にLPレコードの例で述べたように仮に北極で台風が発生して南下したとしてもそれは右に逸れていきます。

 実は北半球ではどこから投げた石も見かけ上,右に逸れていくわけです。

 例えばスナイパー:ゴルゴ13が1km遠方の標的を狙って狙撃しても,弾丸は僅かに右に逸れていくのでそれを勘定に入れて狙う必要があるわけです。

 もしも南半球なら逆に左に逸れるわけですね。

 こうした北半球で右にそれる現象は結局,遠心力などと同じく"見かけの力=慣性力"が働いていると考えることができて,それを発見者の名前にちなんでコリオリ(Coriolis)の力と言います。

 そして北半球での台風を考えると,台風ですから"中心=目の部分"の気圧が最低でまわりの気圧は目の部分のそれより高いわけです。

 風はどのように吹くか,というと水が高いところから低いところへと流れるように,風も気圧の高いところから低いところ目指してその気圧のスロープに沿って吹いていくわけです。

 もしもコリオリの力がなければ,風は"外周部から中心に向かって一直線に進む=落下していく"はずなのですが,コリオリの力によって気圧のスロープも右にねじれてしまっています。

 それ故,風は外周部から中心に向かっていくときに,右にフックして逸れていきながら,最後は中心の気圧最低の目に向かっていくことになり,そのために左巻き(反時計回り)になるのです。

 南半球での台風は逆に右巻き(時計回り)ですね。

 どこかの"バカな大学教授"は,風呂の水が排水口へと流れていき排水されるときに,北半球では左巻き(反時計回り)ですが,赤道を越えて南半球に入ったとたんに右巻き(時計回り)に変わる,などと主張したと聞きましたが,それは誤りですね。

 風呂の水程度の規模では地球自転の影響などは出てきません。

 たまたま排水口付近で左巻きの角運動量を持っていたら左巻きになり,逆なら右巻きになるというだけで,それはカオス現象,つまり偶然の産物でしかありません。

 しかし台風くらいの規模になると地球の自転がもろに効いてきます。

 遠心力の加速度は緯度θでΩ2Rcosθですから,最大でも重力の加速度の0.3%程度です。

 北極で体重100kg重の人が赤道で体重を測っても300g重くらいしか軽くはなりません。一方,コリオリの力の加速度は台風の北上の速度をvとして2Ωsinθ×vです。

 Ω=7,2×10-5/sですから,緯度θが35度で台風の北上速度が100mを10秒で走る程度の時速36km程度なら,加速度a=8.3×10-4m/s2であり,重力g=9.8m/s2の約0.01%程度です。

 そこで,コリオリの力の加速度は最大で重力加速度の0.3%程度しかない遠心力のさらに1/30程度にすぎませんが,台風程度の規模だとそれがかなり効いてきます。

 地球自転の証拠であるとされるフーコー(Foucault)の振り子をこのコリオリ力で説明することもできます。

 ニュートン(Isac Newton)は"慣性系の同等性=ガリレイの相対性原理"は認めても"回転系を含む非慣性系の同等性=マッハ(Mach)原理 → 一般相対性原理"を認めることをあきらめました。

 そして,彼が"絶対座標系=絶対空間"に固執せざるを得なかったのも,こうした"遠心力やコリオリ力の絶対性"を解消する道はない,という考えからだったという話もあります。

 こうした"見かけの力=慣性力"の扱いはとても悩ましいところがあります。

※以上,再掲記事です。

PS:この連休では,2008年11/23の記事「超弦理論(7)(弦の相互作用と頂点演算子)」を再編集するため,図1-5から図1-12まで都合8個も図を入れましたが,これは私には慣れない作業なので疲れました。出来映えはいかに?。。。 

| | コメント (0) | トラックバック (0)

2009年12月10日 (木)

震源の探知(大森公式等)

 今日は,ちょっと手抜きのツナギで2006年11/14の記事 「結晶内での弾性波(地震波)」の続きとして,中学か高校の地学で習った簡単な知識を披瀝してお茶をにごします。

 年末,体調はよくも悪くもないですがヒマ人なりに予定がこみ合っていて貧乏なこともありフリ-マンをエンジョイというわけにはいきません。

 ではまず,記事の再掲から始めます。

 (※再掲)

 今日は等方的とは限らない一般の結晶内での弾性波,特に地震波について考察してみます。 

弾性体の密度をρ,応力テンソルを{σjk},それを構成する部分の局所的速度を={uj}とすると,この弾性体が従うべき運動方程式は流体方程式と同じくρd2j/dt2=∂σjk/∂xkで与えられます。

 

そして,歪み速度テンソル{ujk}をjk≡(∂uj/∂xk∂uk/∂xj)/2で定義すれば,これは{σjk}と同じく反対称の単なる回転を除いた対称テンソルです。これは6個の独立成分を持ちます。

 

微小変位に対しては適切な線形弾性体近似では,フック(Hooke)の法則:σjk=Cjklmlmが成立するため,先の運動方程式はρd2j/dt2=Cjklm(∂2m/∂xk∂xl)となります。

jklmは,一般に81個の成分を持ちますが,{σjk}も{jk}も対称テンソルなのでjとk,lとmについて対称ですから,独立成分は6×6=36個となります。

  

また,歪みエネルギーWはΔW=σjkΔjkから求まり,対称2次形式W=(1/2)jklmjklmになるので,jklm,kとl,mの交換についても対称であり,結局,その独立成分は(30/2)+6=21個となります。

この方程式の平面波の解をj=u0jexp[i(kr-ωt)]として代入すると,ρω2j=Cjklmklmとなります。

 

これは書き直すと(ρω2δjm-Cjklmkl)um0 とも書けます。

 

この1次方程式が自明でない解を持つためには,3行3列の行列係数:(ρω2δjm-Cjklmkl) (j,m=1,2,3)の行列式がゼロ,

 

すなわち,det(ρω2δjm-Cjklmkl)=0 が成立する必要があります。

 

そして,ω2を未知数としてこの方程式を解けば,ω2の3個の解が得られ,それらはの関数となります。 

ところで,もしも結晶が等方性弾性体であれば歪みエネルギー:

W=(1/2)jklmjklmは座標系の回転に対して不変なスカラーであるはずですが,jkの2次形式の形で得られる独立な不変スカラーはukk2jkjkのみです。

 

そこで,適当な係数λ,μを選んでW=(1/2)λukk2+μjkjkと書けるはずです。

 

そこで,応力テンソル{σjk}={Cjklmlm}はσjk=λullδjk+2μujkと書けます。ここに弾性定数λ,μはラメ(Lame)の定数と呼ばれています。

このときには運動方程式も,

ρd2/dt2(λ+μ)∇(∇)+μ∇2

と簡単になります。

 

先に挙げた1次方程式の係数の行列要素は,

ρω2δjm-Cjklmkl(ρω2-μ2jm(λ+μ)jm

となります。

 

の向きをx軸の正の向きに取ると,k1=k=||,k2=k30 ですから,行列[(ρω2-μ2jm(λ+μ)jm]は対角成分が[ρω2(λ+2μ)k2]と2つの(ρω2-μk2)の対角行列になります。

 

したがって,det(ρω2δjm-Cjklmkl)=0 の解は,

P2≡(ωP/k)2(λ+2μ)/ρ,VS2S/k)2=μ

となります。

 

それぞれの速度に対応する平面波は,速度Pで伝播する波の伝播方向への振動である"縦波=P波"と,速度VSで伝播する波の伝播方向に垂直な方向の振動である"横波=S波"です。

 

しかし,もしも異方性の弾性体の場合ならω>0 の3つの解ω()はの関数として1次の同次式ではあっても,単なる1次関数になるとは限らず,弾性波は一般に分散性の波でもあると思われます。

 

そこで,伝播速度は"単一波の速度=位相速度"ではなく,群速度=∂ω/∂で与えられると考えられます。

ところで,σjkやujk[σ]=t112233233112)のように6次元の列ベクトルで表わすボイト(Vogit)の表記という表記法があります。

 

この表記で表現すると,フックの法則は[σ]=C[u]と簡単になります。ここでCは6行6列のテンソル行列です。

 

また,このとき,特に結晶が等方性弾性体ならCの成分のうちで独立なのはやはり2つだけです。

 

先に述べたラメの定数はλ=C12,およびμ=C44=(C11-C12)/2と表わされます。

例えば結晶が立方体構造をしている立方晶系では独立成分は3つです。

 

つまり,C11=C22=C33,C12=C23=C31,かつ4~6行の成分は4~6列しか成分のない対角行列であってC44=C55=C66です。

 

結局,この結晶構造を示すにはC11とC12とC44の3つだけあれば十分である,ということになります。

この条件ではx軸をの向きに取り,先のdet(ρω2δjm-Cjklmkl)= 0 でのCjklmをボイトの表記でCを表わせば,ξ=(ω/k)2,a=C11/ρ,b=C44/ρ,c=C12/ρとして(ξ-a)(ξ-b)2=0 です。

 

この解はξ=a,bとなります。

 

つまり速度をVP,VSとすると,VP2=C11/ρとVS2=C44/ρの2つの速度が得られます。

 

これは等方性弾性体でのC44=μ,C11=λ+2μのケースと一致します。

一方,六方晶系ではC11=C22で,C13=C23,また4~6は対角行列でC44=C55,C66=(C11-C12)/2ですから,結局のところC11,C12,C13,C33,C44の5つだけがあれば十分となります。

 

結果だけ書くと,VP2=C11/ρ,およびVS2=C44/ρ,(C11-C12)/(2ρ)となり,"横波=S波"には2つの速度があるので地震の観測ではS波のほうに二重の波が観測されると予測されます。

 

ただし,異方性の結晶では一般に波は完全にS波とP波になるように行列が対角形になるとは限らないので,これらの計算においては偶々波の進行方向が結晶の対称軸と一致したり直交していたりする特別なケースだけしか扱っていません。

 

一般的な扱いについては,暇があってその気になればまた記事にするかもしれません。 

参考文献;ランダウ=リフシッツ 著「弾性理論」(東京図書),角谷典彦 著「連続体力学」(共立出版)

(再掲終了※)

 と書きましたが,ここからもっと身近な地震の話題へとつなげます。

等方性媒質を仮定すると,P{(λ+2μ)/ρ}1/2,VS=(μ/ρ)1/2ですからVP>VSです。

そこで,震源Oから地震を感知する場所:AまでP波,およびS波が到着するまでの時間を,それぞれTP,およびTSと書けばTP=∫OA(1/VP)ds,TS=∫OA(1/VS)dsです。

 

P>VSですからTS>TPとなります。

つまり,まず縦揺れのP波だけが到着しその後に横揺れS波も到着して両方が重ね合わされた大きいゆれが始まるのですね。 

特にOA間の到るところで媒質が同じ均等な弾性体であれば,この区間でP,VSが一定です。

 

そこで,このときにはOA間の波の進行距離(直線距離でなくてもよい)をdAOAdsとするとTPA/VP,TSA/VSです。

これから,辺々を引き算すると,

 

S-TP(A/VSA/VP)=A(1/VS-1/VP) です。

 

それ故,TSP≡TS-TP,k≡(1/VS-1/VP)-1=VPS/(VP-VS)と定義すれば,A=kTSPという比例関係の公式を得ます。 

 この式は発見者の大森房吉氏(1918)の名前を取って,大森公式(Ohmori's law)と呼ばれます。 

 

 (下図は中越地震本震の地震計記録からです。)

 

     

こうした等方的で均質な媒質を伝わる場合,弱いP波だけが来て「あ,ひょっとして地震かな?」と思ってから本格的で大きな揺れがくるまでの初期微動の時間をTSPとするとその地点Aから震源Oまでの距離AはTSPに比例します。 

初期微動時間が長いほど震源までの距離(地理的な距離だけでなく震源深さまでも含めた距離)が長いということが言えます。もしも震源が直下にあれば縦揺れは上下震動,横揺れは左右震動になります。 

さて,標高hが違う3地点A,B,Cがありその座標がそれぞれ(xA,yA,hA),(xB,yB,hB) (xC,yC,hC)であるとき,上記のような公式に基づいて初期震動時間の観測から観測点から震源までの距離A,B,Cが全てわかったとします。 

このとき,未知の震源Oの座標を(xO,yO,hO)と仮定すれば,地震波が全て直進なら,

 

A2(xA-xO)2+(yA-yO)2+(hA-hO)2,

B2(xB-xO)2+(yB-yO)2+(hB-hO)2,

C2(xC-xO)2+(yC-yO)2+(hC-hO)2

 

が成立します。

これは3つの未知数xO,yO,hOに対する3つの連立2次方程式ですから解くことが可能です。

 

つまり,正確に震源までの距離を観測できる点が3個あれば原理的には震源の位置がわかります。これは3点法と呼ばれるものですね。

普通,震源の高さ:hOは海抜で測って負の数であり,震源は地中または海中にあります。

大森公式:A=kTSPは地震に対する公式で大森係数と呼ばれる

k=VPS/(VP-VS)は8(km/s)程度の値です。

 

しかし,パラメータkを色々と変えると,雷における音と光のように速さの異なる2つの信号が届く現象では全く同じ公式が成立します。

ただし,雷の場合は光速がc~30万km/s,空気中の音速がv~340m/sでc>>vなのでk=cv/(c-v)~cv/c=vですから,事実上A~vTなのでほとんど光速は関係無しです。(Tは光が見えてから音が聞こえるまでの時間です。)

 

例えば,ピカッと光ってから雷の音が聴こえるまでの時間が1秒なら,自分から340m程度離れたところ(空中かもしれません)で"放電=落雷"があったと推測されます。

   

PS:>私のマノン・レスコーたちへ 

 病気とか事故に会っていなければどこで誰と遊んでいてもいいけれど,連絡つかない状況だととにかく心配です。大丈夫かなあ。。。。

 ("マゾ, ロリコン, 準デブ専, メガネフェチ=変態"のTOSHIより。。)

(T.ウッズも聖人君子ではなく適当にストレスを解消しているらしいので安心しました。あれだけの収入を得る能力があれば大勢の家族を養えますね。

 性欲のある健康な男だし,大有名人で品行方正キャラが売り物なら日本だと風俗通いもままならないので女遊びは囲い込むしか仕方ないとして,もしもプロテスタントならマドンナの養子のように大勢を扶養するのはむしろ義務かな?)

PS2:ドクちゃんって日本が好きなんだ。。??

 国ではなくて故郷(くに)が好きなんでしょ?

 http://news.aimu-net.com/read.cgi/wildplus/1257045443/

PS3:「陳情政治」からいきなり「非陳情政治」へと革命的に移行するのではなく,途中で1人(小沢さん)に権力が集中するシステムを経る。。

 共産革命でも過渡期ではプロレタリア独裁(現実には中国では毛沢東,ロシアではレーニンの1人独裁)です。それはそれでマヌーバとしてありなのでしょうが独裁者がプラトン的(プラトニック)に正しい哲人でいるうちはいいけど,うまく移行できないと"元の木阿弥"か,もっと悪くなります。

 小沢氏が中国に大挙して出かけていって米国にブラフかけて「ポッポ政権」を裏で後押ししてるのかも知れないけど,うまくいくのかねえ。基地問題。。

 朝っぱらから女子アナが「吉宗は何将軍?」と聞かれて「暴れん坊将軍?」と答えていました。それでもいいんだろうけど,ニュース番組なので一応「米(コメ)将軍だろ?」とツッコミたくなりました。。。アハハ

クリックして投票してくだかなりさい。(1クリック=1投票です。1人1日1投票しかできません。クリックすると「人気blogランキング」に跳びます。)

にほんブログ村 科学ブログへ にほんブログ村 科学ブログ 物理学へクリックして投票してください。(ブログ村科学ブログランキング投票です。1クリック=1投票です。1人1日1投票しかできません。クリックするとブログ村の人気ランキング一覧のホ-ムページームに跳びます。)

iconヤーマン プラチナゲルマローラー 1日3分コロコロエステ!ローラー型プラチナ配合美顔器  

ブックオフオンライン 

お売りください。ブックオフオンラインのインターネット買取 展開へ! ▼コミック 尾田栄一郎 「ONE PIECE(52)」 icon ▼コミック 「ONE PIECE」をオトナ買い icon

三国志特集 ▼コミック 横山光輝 「三国志全巻セット」 icon 「三国志(文庫版)全巻セット」 icon  「三国志(ワイド版)全巻セット」 icon  ▼書籍 「三国志」/吉川英治 icon  「三国志」/北方謙三 icon  「三国志」/宮城谷昌光

iconオンライン書店 boople.com(ブープル) 

| | コメント (0) | トラックバック (0)

2009年11月30日 (月)

定量的地震学6

 つなぎで地震学の続きを書きます。

前回は内部面上での応力,または変位(歪み)の不連続性で表現される地震源(面源)の体源による等価物の存在を証明しました。

そして,断層面を横切る剪断作用に対する体積力の等価物として次のような表現が得られるという結論で終わりました。 

まず,応力の不連続性[]がゼロでない場合には,その寄与-∫-∞dτ∫Σ[Tp((ξ,τ),ν)]Gnp(ξ,t-τ;,0)}dΣξは,-∫-∞dτ∫V(∫Σ{[Tp((ξ,τ),ν)]δ3(ηξ)dΣξ}Gnp(ξ,t-τ;,0))dVηで表現できます。

それ故,Σ上の応力不連続性の体積力等価物は[T](η,τ)≡-∫Σ[Tp((ξ,τ),ν)]δ3(ηξ)dΣξで与えられます。

一方,変位の不連続性による寄与∫-∞dτ∫Σ[ui(ξ,τ)]Cijpq(ξj{∂Gnp(ξ,τ-t;,0)/∂ξq}は,-∫-∞dτ∫V (∫Σ[ui(ξ,τ)]Cijpq(ξj{∂δ3(ηξ)/∂ηq}dΣξ)Gnp(ξ,t-τ;,0)dVηで表現されます。

そこで,Σ上の変位不連続性の体積力等価物はfp[](η,τ)≡-∫Σ[ui(ξ,τ)]Cijpq(ξj{∂δ3(ηξ)/∂ηq}dΣξで与えられると考えられます。

ただし,鍵括弧[ ]の量はΣ上の各点での不連続性を表わす量です。

 

すなわち,[(ξ,τ)]≡(ξ,τ)|Σ+(ξ,τ)|Σ-,かつ[ui(ξ,τ)]≡ui(ξ,τ)|Σ+-ui(ξ,τ)|Σ-etc.と定義されています。

 

という内容のことを書きました。

こうすれば応力,および変位の不連続性はそれぞれ-∫-∞dτ∫Σ[Tp((ξ,τ),ν)]Gnp(ξ,t-τ;,0)}dΣξ=∫-∞dτ∫V[fp[T](η,τ)Gnp(,t-τ;η,0)]dVη,および∫-∞dτ∫Σ[ui(ξ,τ)]Cijpq(ξj{∂Gnp(ξ,τ-t;,0)/∂ξq}=∫-∞dτ∫V[fp[](η,τ)Gnp(,t-τ;η,0)]dVηと体源で表現できるわけです。

まず,これらについて前回,少し書き残したことを補足することから始めます。

 上記,変位不連続性の表現の被積分関数は,各pに対して,添字i,j,qのそれぞれが異なる27個の項を含んでいますが,以下では媒質に対称性があって2個か3個を除く全ての項がゼロであるような重要な例を取り上げる予定です。

 しかし,取り合えず上記の恒等式自体は一般的な非均質,非等方媒質に対して成立する式です。

 

 これらはV内の各点における体源応力が断層面の上の弾性媒質の性質のみに依存する表現になっていることは注目すべきことです。

そして,V内の断層運動は内部プロセス(内力のみの系)なので,全運動量と全角運動量は保存されなければなりません。

すなわち,全てのτに対し∫V[](η,τ)dVη=0,かつ全てのτとη0に対し∫V[(ηη0[](η,τ)dVη=0 であるべきです。

実際,fp[](η,τ)≡-∫Σ[ui(ξ,τ)]Cijpq(ξj{∂δ3(ηξ)/∂ηq}dΣξですから,ガウスの積分定理によって∫Vp[](η,τ)dVη=-∫VdVηΣdΣξ{[([ui(ξ,τ)]Cijpq(ξj{∂δ3(ηξ)/∂ηq}=∫ΣdΣξSextdSη[ui(ξ,τ)]Cijpq(ξjδ3(ηξ)です。

 

ところがSextとΣは全く共通点がないので,決してηξとは成り得ず,最右辺は確かに消えます。

一方,∫V[(ηη0[](η,τ)dVηの第m成分は∫Vεmnpn-η0n)fp[](η,τ)dVη=-∫ΣdΣξ[ui(ξ,τ)]Cijpq(ξjVdVηmnpn-η0n){∂δ3(ηξ)/∂ηq})=∫Σmqpijpq(ξj[ui(ξ,τ)])dΣξで与えられます。

 

ところがCijpq(ξ)ijqp(ξ),かつεmqp=εmqpなので,これの右辺もゼロになります。

こうして体積力の表現が内力である条件を満たす合理的なものであることを示すことができました。これが前回書き残した補足事項です。 

さて,野外の不連続性に等価な体積力の簡単な例として,丁度1点と1方向に加えられた体積力のケースを考えます。

以前の記事「定量的地震学1,3」では次のように瞬間体積力を与えました。

ξに対する1つの個別粒子に時刻t=τにn方向に瞬時的に加えられる1つの体積力(,t)があれば,その成分fi(,t)は空間位置を与えるのに3次元のディラック(Dirac)のデルタ関数,衝撃の時刻を与えるのに1次元のデルタ関数に比例してfi(,t)=Aδ3(ξ)δ(t-τ)δinと表現されます。

ただしAは衝撃の強さを与える定数です。fi3(ξ),δ(t-τ)の次元がそれぞれ[力/体積]=MLT-2/L3,1/L3,1/Tであることに着目するとδinは無次元なので,衝撃の強さAは正しく,”衝撃=力積”の物理的次元を有することがわかります。”

これは,応力成分の不連続性とみなすこともできます。

ここでは,体積力の応力の不連続性との等価性を見るために,x3を深さ方向とし,x3=hの1点(0,0,h)に加わる体積力ですがδ(t-τ)に比例する瞬時的な力ではなく,τ=0から定常的にx3=hの1点(0,0,h)にFの大きさの体積力とします。

すなわち,(η,τ)=δ(η1)δ(η2)δ(η3-h)θ(τ),≡(0,0,F)とします。ここでθ(τ)はHeaviside関数(階段関数)です。

これは,平面ξ3=h上の1点を横切る応力の不連続性:[(ξ,τ)]=(ξ,τ)|(ξ1,ξ2,h+)(ξ,τ)|(ξ1,ξ2,h-)=-δ(ξ1)δ(ξ2)θ(τ)と同一視できることがわかります。

これは,先に得られた表現:[T](η,τ)≡-∫Σ{[Tp((ξ,τ),ν)]δ3(ηξ)dΣξの右辺において,平面ξ3=hをΣとして,[Tp((ξ,τ),ν)]=-δ(ξ1)δ(ξ2)θ(τ)を代入すれば,[T](η,τ)=δ(η1)δ(η2)δ(η3-h)θ(τ)が確かに得られることからわかります。

断層スリップ(地滑り)によって引き起こされる地震波は,モーメントを帳消しにするようなある断層上の力の分布によって引き起こされる地震波と同じです。

 

問題としている断層スリップに対し,この力の分布は一意的ではありませんが,等方的な媒質内ではそれを常に複数の偶力源の表面分布として選択可能です。

この結論は,地震が単一の偶力源か複数の偶力源のいずれによってモデル化されるかという疑問に関連して長期間論じられた論争の見地からは皮肉な結論です。

単一の偶力源を擁護する人々は地震は断層上のスリップが原因であると強く信じていました。そして彼らは,これが単一の偶力源,すなわち断層の相対する側の運動に対応する2つの力に等価であると直線的に信じました。

しかし,弾性力学では経験上直線的アプローチは危険なことが多いようです。

複数の偶力源を擁護する人々のあるものは,地震が既存の剪断応力下での体積的崩壊であるに違いないと思っていました。 

今や,複数の偶力源に等価であると認識されている近年の断層理論は,莫大な距離で観測された放射パターンによる支持のみならず,震源に非常に近いところで得られたデータによる強い支持を得ています。 

断層Σは平面ξ3=0 上にあるものとします。すると,スリップ[]はξ3方向には成分を持たず,ベクトル[]は平面Σに平行と考えられます。 

さらに,平面ξ3=0 上でのスリップ[]の方向をξ1方向とします。すると,[2]=[3]=0 でありν1=ν2=0 ,ν3=1 です。 

そこで,等価体積力の表現:p[](η,τ)=-∫Σ[ui(ξ,τ)]Cijpq(ξj{∂δ3(ηξ)/∂ηq}dΣξは,fp[](η,τ)=-∫Σ[u1(ξ,τ)]C13pq(ξ){∂δ3(ηξ)/∂ηq}dξ1dξ2となります。  

 

既述のように,不均質でも等方的な物質では弾性係数CijpqはCijpq=λδijδpq+μ(δipδjq+δiqδjp)なる形をしています。独立定数λ,μはLame(ラメ)の定数と呼ばれます。

 

よって,今の場合の被積分関数のC13pq(ξ)は,C1313(ξ)=C1331(ξ)=μ(ξ)を除いて全てゼロです。

それ故,f1[](η,τ)=-∫Σμ(ξ)[u1(ξ,τ)]δ(η1-ξ1)δ(η2-ξ2){∂δ(η3)/∂η3}dξ1dξ2,f2[](η,τ)=0,f3[](η,τ)=-∫Σμ(ξ)[u1(ξ,τ)]{∂δ(η1-ξ1)/∂η1}δ(η2-ξ2)δ(η3)dξ1dξ2となります。

まず,f1[](η,τ)=-∫Σμ(ξ)[u1(ξ,τ)]δ(η1-ξ1)δ(η2-ξ2){∂δ(η3)/∂η3}dξ1dξ2を見ると,これは±η1方向内の力でη2方向に腕がありη3方向にモーメントを持つΣにわたる単偶力源を示しているとわかります。

実際,積分を実行すると,f1[](η,τ)=-μ(η)[u1(η,τ)]{∂δ(η3)/∂η3}です。

 

これは平面η3=0+上に寄与する点力とη3=0-上に寄与する反対向きの点力のペアであると考えられ,補足事項の最初の式で示したように,内力なのでf1[](η,τ)の合力成分は消えます。

一方,力のモーメントはこの力の成分だけでは消えず,この成分によるη2軸の回りのモーメントは∫Vη31[]dV=-∫Vη3μ(η)[u1(η,τ)]{∂δ(η3)/∂η3}dη1dη2dη3=∫Σμ(ξ)[u1(ξ,τ)]dΣξとなります。

Σにわたるスリップを平均すると,<u(τ)>≡∫Σ[u1(ξ,τ)]dΣξ/Sです。ただし,S≡∫ΣdΣξは断層の全面積です。

そこで,もしも断層域が均質,つまりΣの上でμ(ξ)=μ(一定)なら,η2軸の回りのf1[](η,τ)による全モーメントは,単にη2軸の正の向きに∫Vη31[]dV=∫Σμ(ξ)[u1(ξ,τ)]dΣξ=μS<u(τ)>で与えられます。

体積力には,また3軸成分:f3[](η,τ)=-∫Σμ(ξ)[u1(ξ,τ)]{∂δ(η1-ξ1)/∂η1}δ(η2-ξ2)δ(η3)dξ1dξ2=-∂{μ(η)[u1(η,τ)]/∂η1}δ(η3)があります。

この成分によるη2軸の回りの力のモーメントは-∫Vη13[]dV=∫Vη1{μ(η)[u1(η,τ)]/∂η1}δ(η3)dη1dη2dη3=-∫Σμ(ξ)[u1(ξ,τ)]dΣξです。

 

これについても,もしも媒質が均質でΣの上でμ(ξ)=μ(一定)なら,-μS<u(τ)>です。

よって,均質,不均質いずれにしても全モーメントはゼロです。これも既に,補足事項として一般的な形で証明した事実に一致しています。

途中ですが今日はここで終わります。 

参考文献:K.Aki,& P.G.Richards 「Quantitative Seismology(Theory and Method)」

  

 クリックして投票してください。(1クリック=1投票です。1人1日1投票しかできません。クリックすると「人気blogランキング」に跳びます。)

にほんブログ村 科学ブログへ にほんブログ村 科学ブログ 物理学へクリックして投票してください。(ブログ村科学ブログランキング投票です。1クリック=1投票です。1人1日1投票しかできません。クリックするとブログ村の人気ランキング一覧のホ-ムページームに跳びます。)

iconヤーマン プラチナゲルマローラー 1日3分コロコロエステ!ローラー型プラチナ配合美顔器  

ブックオフオンライン 

お売りください。ブックオフオンラインのインターネット買取 展開へ! ▼コミック 尾田栄一郎 「ONE PIECE(52)」 icon ▼コミック 「ONE PIECE」をオトナ買い icon 

三国志特集 ▼コミック 横山光輝 「三国志全巻セット」 icon 「三国志(文庫版)全巻セット」 icon  「三国志(ワイド版)全巻セット」 icon  ▼書籍 「三国志」/吉川英治 icon  「三国志」/北方謙三 icon  「三国志」/宮城谷昌光

iconオンライン書店 boople.com(ブープル) 

icon  iconicon  icon icon icon予約】セックス・アンド・ザ・シティ ザ・ムービー コレクターズ・エディション<限定盤>

| | コメント (0) | トラックバック (0)

2009年11月17日 (火)

定量的地震学5

 1月以上の間があきましたが地震学の続きです。 

 前回は曲線座標変換というやや幾何学的な話に寄り道しましたが,今日は第3章の地震源(震源)の表示に入ります。

 固体地球外部の地震源の例を上げると,風,大洋波浪,隕石衝突,ロケットの打ち上げ,大気中の爆発などがあります。その他,人々の歩行によってさえ地震が生じます。

こうした外部の地震源については,通常単純な地球表面に加えられた時刻変動応力の解析の枠組みの中で処理できます。

 また,多くの実用的目的を有する現象に起因するものや,その他にも火山の噴火,ガス漏れなど比較的規模の大きい爆発,粉砕などもありますが,これも外部ソースであり上記に含まれるでしょう。

 一方,一般的な地震や地下爆発など地球内部のソースに対しては,解析の枠組みはよりむずかしい展開を必要とします。

なぜなら,これまで論じてきた弾性運動を支配する方程式が固体地球内部の到るところで有効というわけではないからです。

この章では内部ソースについてのみ論じます。これは断層源と体源の2つのカテゴリーに分けられます。 

断層源は破砕された平面を通っての"断層=地滑り"のような内部面と関連した事象です。一方,体源は体積的なソース領域における爆発的膨張のような内部体積に関わる事象です。

内部ソースの数学的記述は古典的には異なる2つのラインに沿って追跡されてきました。

 

1つはソースを含む媒質のある要素へ加わる体積力によるライン,もう1つは変位,または歪みの幾つかの不連続性(すなわち,断裂する断層面や体源の表面を通るそれ)によるラインです。

しかし,2番目のアプローチは有効的に1番目のそれに組み込むことが可能です。すなわち,断層面を横切る単純な剪断作用については体積力の等価物が存在します。

 

以下では,こうした体積力の等価物について理論を幾分詳細に展開し,根本的に異なる力の系が正確に同じ変位の不連続性に同等であることを示すことから解析を始めます。 

それから後に,BarridgeとKnopoff(1964)に従って断層源の一般理論を展開し,最後に体源に関する理論を概説する予定です。

一般に,震動図に記載された運動は地震源の効果と伝播の効果の両方の結果です。

 

震動図を理解することを追求してきた主な理由の1つは地震運動から伝播の効果を分離することです。というのは地震源の効果は地球の内部構造に関する情報を携帯しているからです。

最近では,地殻プレートの運動を図示して地震源のプレートがどのように動かされるかを学ぶという目的のために,次震源(震源)のメカニズムが研究されているようです。

現在では,近くの断層の性質と局地的応力の分布に関する地質学的,地球物理学的データに基づいて工学用途の敷地における地震危害の予測を行なうという観点から,こうしたソース(source;震源)の理論が広く展開されているようです。

さて,初めに述べた通り,[内部面における表示定理]="応力と変位における不連続性に対する体積力(実体力)の恒等式(等価物)"を与えるという主題に移ります。

第2章で与えられた表示定理は,表面Sを体積V内部の隣り合う面に選べば震源において強力な助けに成り得ます。

ここでの考察の動機は,H.F.Reidの仕事に起因するものです。

 

彼の1906年のSan Francisco地震前後のSan-Andreas断層の研究は,地震運動が活性化した地質学的断層上の自然発生的地滑りにより放射された波動によるという一般的認識に導きました。

 

我々は,こうしたソースのメカニズムをすぐ後により詳細に論じる予定です。他方,動力学的過程や他のソースのメカニズムについてはずっと後の第15章で述べる予定です。

現在の関心としては,埋没した地滑りのプロセスとそれから放射された波が,既に得られた表示定理からどのようにして自然に解析され得るかということです。 

既に述べたように,表示定理は次の3つの異なる表現を有します。 

n(,t)=∫-∞dτ∫V[fi(ξ,τ)Gin(,t-τ;ξ,0)]dVξ+∫-∞dτ∫S[{Gin(ξ,t-τ;,0)Ti((ξ,τ),)-ui(ξ,τ)Cijkl(ξ)nj{∂Gkn(ξ,τ-t;,0)/∂ξl}}dSξ..(1)

n(,t)=∫-∞dτ∫V[fi(ξ,τ)Grigin(,t-τ;ξ,0)]dVξ-∫-∞dτ∫S[ui(ξ,τ)Cijkl(ξ)nj{∂Grigkn(,τ-t;ξ,0)/∂ξl}]dSξ..(2)

n(x,t)=∫-∞dτ∫V[fi(ξ,τ)Gfreein(,t-τ;ξ,0)]dVξ+∫-∞dτ∫S{Gfreein(,t-τ;ξ,0)Ti((ξ,τ),)}dSξ..(3) です。

そこで述べたように,これらは変位ベクトル(,t)がS上の変位に依存するのか →(2),それとも応力に依存するのか →(3),または両方に依存するのか →(1)という疑問への矛盾を示しているように見えます。

  

しかし,弾性媒質上では応力と変位が独立というわけではないので矛盾はありません。

さて,この表示定理における表面SがVの外部境界面だけでなく埋没した断層の相対する面であるような隣り合う内部面をも含むとする視点は地震源の理論にとって中心的な論点です。

すなわち,S≡Sext+Σに取ってみます。ただし,SextはVの外部表面でありΣ≡Σ++Σ-はV内部の1つの断層の相対する面です。

まず,表示定理で最も一般的な式(1)をun(,t)=∫-∞dτ∫V[fp(ξ,τ)Gnp(,t-τ;ξ,0)]dVξ-∫-∞dτ∫S[ui(ξ,τ)Cijpq(ξ)nj{∂Gnp(ξ,τ-t;,0)/∂ξq}-{Gnp(ξ,t-τ;,0)Tp((ξ,τ),)}]dSξと書き直します。

 Sextは地球全体の表面としてもいいのでそう考えます。

 

 グリーン関数Gがext上で斉次境界条件Ai+γji,j=0,Bi+γji,j=0 を満足していると考えられます。ただし,Ai(,t)≡Gim(,t;ξ11),Bi(,t)≡Gin(,t;ξ2,-τ2)です。

こうすれば,右辺のS=Sext+Σにおける表面積分のうち,外部地球の表面Sextの寄与は無視できます。 

そして,Σ-の外向き法線単位ベクトルをn=νとします。するとΣ+の外向き法線単位ベクトルは明らかにn=-νです。 

また,Vの内部での一般座標をηとし,ξは表面(断面)Σの上の一般座標のみを表わすとします。

さらに,鍵括弧[ ]の中の量はΣ+とΣ-の差を示すとします。例えば,[(ξ,τ)]≡(ξ,τ)|Σ+(ξ,τ)|Σ-と定義されます。

すると,表示定理はun(,t)=∫-∞dτ∫V[fp(η,τ)Gnp(,t-τ;η,0)]dVη+∫-∞dτ∫Σ[ui(ξ,τ)Cijpq(ξj{∂Gnp(ξ,τ-t;,0)/∂ξq}-{Gnp(ξ,t-τ;,0)Tp((ξ,τ),ν)}]dΣ となります。

これが実際的な計算において意味を持つためには,Σ上で何らかの境界条件を与える必要があります。

境界条件のに対する選択は破裂する断層面を横切る応力と変位の現実の性質に従わせる必要がありますが,Gに対する選択は有益な形になるように自由に選ぶことが可能です。

 

変位については,両断層面の上でゼロではないスリップ[]が存在するはずですが,応力については連続性から[(,ν)]=0 と考えられます。

また,Σ上でのGの性質の定義を確立する最も簡単で最も共通に用いられている方法は,Σを横切ってGとその空間微分係数が連続であるように設定することです。 

これは体積Vに対して計算するには最も容易なグリーン関数です。

 

もしも,このに対して体積力がないときにはn(,t)=∫-∞dτ∫Σ([ui(ξ,τ)]Cijpq(ξj{∂Gnp(ξ,τ-t;,0)/∂ξq})dΣ となります。

 この表現式は断層上での変位があらゆる場所における変位を決めるのに十分であることを示していますが,これは「定量的地震学2」で述べた"一意性定理"から予測できることです。

 

(↑ 複素関数論のコーシーの積分定理にも似てますね。)

※注:"一意性定理"の内容を再掲します。

[一意性定理]:表面境界Sを持つ体積Vの弾性体内の到るところで,与えられた時刻t0における変位と粒子速度の値,およびt>t0における(ⅰ)体積力と供給される熱,(ⅱ)表面S=S1+S2の一方の部分S1上での応力Π,(ⅲ)残りの部分S2上での変位が既知なら,時刻t0より後の時刻tおける変位(,t)は一意的に決まる。

 

(注終わり※)

 

しかし,一見したところこの表現式においてソースの伝播を記述するグリーン関数についてΣにおける境界条件が与えられていないことは驚くべきことです。

断層上で生じた運動は,それ自身が断層面によってある形式で回折される波を作り上げると予測されますが,この相互作用はスリップ関数[(ξ,τ)]の決定を複雑にはするけれどグリーン関数の決定には入り込まないからですね。

こうして,多くの地震学者はある仮定されたスリップ関数によって作り上げられた運動を記述するために,この公式:un(,t)=∫-∞dτ∫Σ([ui(ξ,τ)]Cijpq(ξj{∂Gnp(ξ,τ-t;,0)/∂ξq})dΣを用いてきました。

以下では,たった今記述した直接には如何なる体積力も含んでない地震源の表現の体積力の等価物を求めることを試みます。

 

n(,t)=∫-∞dτ∫Σ([ui(ξ,τ)]Cijpq(ξj{∂Gnp(ξ,τ-t;,0)/∂ξq})dΣは,各々が体積力によって作り上げられたグリーン関数にわたって積分することにより,(,t)における変位を与えることを示していると見ることもできます。

 

こうした見方によれば,活断層も体積力の面的寄与と見なせるような意味があるに違いないと確信されます。

この方法で予測されるような体積力の等価物を決定するため,再び元の一般的な表示定理の式(1)の表現から始めます。

すなわち,再掲するとun(,t)=∫-∞dτ∫V[fp(η,τ)Gnp(,t-τ;η,0)]dVη+∫-∞dτ∫Σ[ui(ξ,τ)Cijpq(ξj{∂Gnp(ξ,τ-t;,0)/∂ξq}-{Gnp(ξ,t-τ;,0)Tp((ξ,τ),ν)}]dΣです。

ただし,GはΣに対してなお透明,つまりΣを横切ってGもその空間微分係数も連続とする仮定を採用します。

 

ここで,Σを横切る[],[(,ν)]については,先の境界条件[(,ν)]=0 のような条件を全く仮定せず,応力の不連続性も存在するようなより一般的な状況とします。

 

するとun(,t)=∫-∞dτ∫V[fp(η,τ)Gnp(,t-τ;η,0)]dVη+∫-∞dτ∫Σ([ui(ξ,τ)]Cijpq(ξj{∂Gnp(ξ,τ-t;,0)/∂ξq}-[Tp((ξ,τ),ν)]Gnp(ξ,t-τ;,0))dΣと書けます。

 

こで,デルタ関数δ3(ηξ)を用いればΣにおける不連続性はVの中に局在化できることを利用します。

 

例えばΣにおける(応力×面積=力):[]dΣは,体積力としての寄与∫VdV{[3(ηξ)dΣ}を与えます。

したがって,応力の不連続性[]がゼロでない場合,その寄与-∫-∞dτ∫Σ[Tp((ξ,τ),ν)]Gnp(ξ,t-τ;,0)}dΣは-∫-∞dτ∫V(∫Σ{[Tp((ξ,τ),ν)]δ3(ηξ)dΣ}Gnp(ξ,t-τ;,0))dVで表現できます。

それ故,Σ上の応力不連続性の体積力の等価物は[T](η,τ)≡-∫Σ{[Tp((ξ,τ),ν)]δ3(ηξ)dΣで与えられます。

変位の不連続性は応力よりも物理的解釈が困難なのですが,数学的に考えて∂Gnp(ξ,τ-t;,0)/∂ξq=-∫V{∂δ3(ηξ)/∂ηq}Gnp(ξ,τ-t;,0)dVなる恒等式を用いれば応力不連続性に等価式に同等な表現を得ます。

すなわち,変位の不連続性による寄与∫-∞dτ∫Σ([ui(ξ,τ)]Cijpq(ξj{∂Gnp(ξ,τ-t;,0)/∂ξq}は,-∫-∞dτ∫V(∫Σ{([ui(ξ,τ)]Cijpq(ξj{∂δ3(ηξ)/∂ηq}dΣ)Gnp(ξ,t-τ;,0)}dVと表現されます。

そこで,Σ上の変位不連続性の体積力の等価物はfp[](η,τ)≡-∫Σ{[([ui(ξ,τ)]Cijpq(ξj{∂δ3(ηξ)/∂ηq}dΣで与えられることがわかります。

今日はここで終わります。 

参考文献:K.Aki,& P.G.Richards「Quantitative Seismology(Theory and Method)」

 

PS:昨日はヘルパーの現場実習の初日として赤羽の訪問介護施設の現役のヘルパーさんに同行して自宅訪問の実習を受けてきました。

 

 実は,実際の利用者様に迷惑をかけてはいけない実習なので,ほとんど見学同然と思ってはいましたが,当初から最も危惧していたのがこの同行訪問実習でした。

 

 というのも,この訪問介護実習で不可避な自転車での移動というのは「糖尿病+心不全+足の動脈硬化」である私にとって可能かどうか?という体力的に最も不安な材料だったからです。

 

 実際,その通りでしたね。

 

 北区の赤羽付近は坂が多く"心臓破りの坂"もあって,心臓+足が動かず,平気を装ってはいましたが一瞬死ぬかとも思うこともありました。

 

 楽な下りがあるということは,帰りは上りですから,動かなくなってしまったら同行者の方に迷惑だと思いましたが,降りて歩いて押していくしかありません。(逆に1人だったら少しは楽だったかも。。。)

 

 利用者様の自宅までの往復と買い物支援,さらに終わって最後に事務所に帰る際,はぐれて道に迷ったのも含め,約2時間も自転車に乗ったのは過去の健康時代を考えても,めったにないことでした。

 

 ともあれ,雨も降らず鬼門の同行訪問実習は無事終了しました。移動以外の実習は体力的には休息時間でしたが,問題ないと思います。

(例によって,何事があっても常にニコニコ笑顔だった(ツモリな)のは,八方美人の性分ですね。)

 

 後の施設実習でも体力は必要でしょうが,まだ給金を頂いて責任を取ることを余儀なくされる本格的な仕事ではなく,学生としての実習だろうし私の最も得意な?他人とのコミュニケーション(=おしゃべり?)の方が主体ではないかとも思い,とにかく少なくとも自転車での移動がないであろうことを考えて楽観しています。

 

 本日は寝て曜日で完全休養日にしよう。アレ?岡山のおフクロの89歳の誕生日は今日17日か21日かどっちだったかな?

 

(↑ イキアタリバッタリで脳天気だなあ。。ウーン,キム・ヨナはいいなあ。。。忘れてた。今日(11/17)は「将棋の日」でもありました。)

  

 (実は,昨日朝は訪問介護事務所の付近まで最近懇意にして頂いている自宅の隣人のMさんに車で送って頂いたのです。Mさんどうもありがとうございました。

 

 しかし,住所だけを頼りにして到着したためか,頂いた矢印付きの地図と写真の通りに徒歩で行けば自然にわかったらしい事務所の入口を間違えてしまって指定時間通りには入ることができず最初から躓きました。

 

 実は私はかなりヒドイ方向音痴です。事務所の担当の方とスクールで対応して頂いた方,ゴメンナサイ。

 

 16年くらい昔は住所だけを頼りに当日自分が交通ガードマンとして立つべき場所を探していたのですから,まるっきりの方向音痴でもないでしょうが。。)

 

 オマケに昼休み近くのコンビニまで往復の際,靴の左側を間違えて履いていって帰ってから気付いて皆様に笑われましたが,靴を間違えてご迷惑かけてゴメンナサイ。

 

 こういうことで笑いを取ってはいけませんね。^^;)

 

PS2:市橋事件について,何日も食事を取っていないなどの報道はどこまで真実なのでしょうか? 取り調べ段階での官憲の守秘義務とマスコミリークとの関係等々についてはよく知りませんが。。。

  

 クリックして投票してください。(1クリック=1投票です。1人1日1投票しかできません。クリックすると「人気blogランキング」に跳びます。)

にほんブログ村 科学ブログへ にほんブログ村 科学ブログ 物理学へクリックして投票してください。(ブログ村科学ブログランキング投票です。1クリック=1投票です。1人1日1投票しかできません。クリックするとブログ村の人気ランキング一覧のホ-ムページームに跳びます。)

iconヤーマン プラチナゲルマローラー 1日3分コロコロエステ!ローラー型プラチナ配合美顔器  

ブックオフオンライン 

お売りください。ブックオフオンラインのインターネット買取 展開へ! ▼コミック 尾田栄一郎 「ONE PIECE(52)」 icon ▼コミック 「ONE PIECE」をオトナ買い icon 

三国志特集 ▼コミック 横山光輝 「三国志全巻セット」 icon 「三国志(文庫版)全巻セット」 icon  「三国志(ワイド版)全巻セット」 icon  ▼書籍 「三国志」/吉川英治 icon  「三国志」/北方謙三 icon  「三国志」/宮城谷昌光

iconオンライン書店 boople.com(ブープル)

| | コメント (0) | トラックバック (1)

その他のカテゴリー

001. 目次 | 002. 募金・ボランティア | 003. 日記・回想 | 004 訃報 | 005. 心身・思想・哲学 | 006. 社会・経済・政治 | 007. 病気(診察・薬) | 008. 恋愛・異性 | 009 宗教・神話 | 010 歴史(日本,世界) | 011. 将棋 | 012. TV(ニュース・ドラマ) | 013 スポーツ(ニュース・イベント) | 014 ノン・フィクション | 015 小説・詩・評論 | 016 漫画・劇画・アニメ | 017 演劇・映画・舞踊 | 018 音楽(日本・西洋・他) | 019 タレント(俳優・お笑い) | 020 ミュージシャン | 021 アイドル・ヒーロー | 022 創作 | 023 シャレ・ギャグ等 | 024 競馬・toto・賭け事 | 025 ファッション・風俗 | 100. 物理学一般 | 101 教育・学校(物理) | 102. 力学・解析力学 | 103. 電磁気学・光学 | 104. 熱力学・統計力学 | 105. 相対性理論 | 106. 星・ブラックホール・一般相対性 | 107. 重力・宇宙・一般相対性 | 108. 連続体・流体力学 | 109. 物性物理 | 110. 複雑系・確率過程・非線型・非平衡 | 111. 量子論 | 112. 原子・分子物理 | 113. 原子核物理 | 114 . 場理論・QED | 115. 素粒子論 | 116. 弦理論 | 118. 観測問題・量子もつれ | 119. 電気回路 | 200. 問題・解答 | 201. 自然科学一般 | 202. 気象・地学・環境 | 203. 生物学・生理学・生化学 | 204. 経済学(ミクロ・マクロ・マルクス) | 300 数学一般・算数 | 301. 集合・位相 | 302. 論理学・数学基礎論 | 303. 代数学・数論 | 304. 解析学 | 305. 複素数・複素関数論 | 306. 線型代数学 | 307. 幾何学(トポロジー・他) | 308. 微分方程式 | 309. 確率・統計 | 310. 関数解析・超関数 | 311 .数値計算・調和解析・離散数学 | 312. 公式・特殊関数 | 501. 商用宣伝・アフィリエイト