輸送方程式のソース項の実装について

[初心者]

始めまして。初心者の者です。

次の講習会資料に興味があり勉強してみたいと思いました。

・2012年4月22日 OpenFOAM勉強会 for beginner@関東　今野 雅(東京大学)

　「輸送方程式のソース項の実装」

　　　http://www.slideshare.net/MasashiImano/20120421-of-beginner


この講習会での計算ケースファイルは公開されているのでしょうか。

「buoyantBoussinesqSimpleFoam」に湿度場を追加しているのは分かるのですが、

湿度の実装では、公式で、

　pw＝((101325+p_rgh)*w)/(.62198+.37802*w)

　pws = exp(?0.58002206*10^4*T^?1 + 0.13914993*10 ? 0.48640239*10^?1*T + 0.41764768*10^?4*T^2 ? 0.14452093*10^?7*T^3 + 0.65459673*10*ln(T))

　RH=pw/pws (%)

として相対湿度RHを計算していると思うのですが、どのファイルにどのように記載すれば
良いのかが経験不足では分かりませんでした。


また、「ＰＭＶ」算出も公式で、
 
PMV=(0.303e-0.36M+0.028)[(M-W)-3.05*10-3*{5733-6.99(M-W)-pa}-0.42*{(M-W)-58.15}-1.7*10-5M(5867-pa)-0.0014M(34-ta)-3.96×10-8fcl×{(tcl+273)4-(tr+273)4}-fcl･hc(tcl-ta)]
 
  ただし
  tcl=35.7-0.028(M-W)-Icl[3.96*10-8fcl*{(tcl+273)4-(tr+273)4}+fcl･hc(tcl-ta)]
  hc=2.38(tcl-ta)0.25　 　　　　2.38(tcl-ta)0.25 >12.1(var)0.5  のとき
  hc=12.1(var)0.5　　　  　　 　2.38(tcl-ta)0.25 <2.1(var)0.5  のとき
 
  fcl=1.00+1.290Icl　　　　Icl<0.078m2･℃/Wのとき　　つまり0.5cloより小さい
  fcl=1.05+0.645Icl 　　 　Icl>0.078m2･℃/Wのとき　　つまり0.5cloより大きい
 
  　M:体表面積の新陳代謝率(1met=58W/m2)
  　W:外部仕事（ほとんどの活動では0）[W/m2]
  　Icl:衣服の熱抵抗[m2･℃/W]　　←clo単位ではありません
  　fcl:裸体時皮膚表面積に対する着衣時の皮膚露出表面積の割合
  　ta:空気温度（℃）
  　tr:平均輻射温度（℃）
  　var:相対風速(m/s)
  　Pa:水蒸気分圧(Pa)
  　hc:対流熱伝導率[W/(m2･℃)]
  　tcl:衣服の表面温度（℃）
 
  PMV:+3（暑い）、+2（暖かい）、+1（やや暖かい）、0（どちらでもない）、-1（やや涼しい）、-2（涼しい）、-3（寒い）

なのは分かるのですが、どのファイルにどのように記載すれば良いのかが
経験不足では分かりませんでした。
trの平均輻射温度はMRTの事でしょうか？グローブ温度や放射温度の計算が
分かりませんでした。


さらに、漏気もどのように設定しているのかが分かりませんでした。

相当隙間面積αAを壁面の横や縦で５等分程度に分割して隙間に５箇所程穴を開けているのでしょうか。


熱・湿気・漏気・人体からの発熱湿などの複合的な計算事例であり、
勉強して理解してみたいと思いました。

ケースファイルが公開されていれば教えて頂きたく存じます。
無い場合は似たようなケースファイルが存在すれば教えた頂きたく存じます。
大変有意義な資料だと思いましたので質疑させて頂きました。

[初心者]

今野　様

大変失礼致しました。

前回は実際の部屋のstlファイルでの解析で、今回の資料は

簡略化したｓｔｌファイルかと勘違いしており、簡略化した講習会

資料と思っていました。

PMVは式を見ると収束計算が必要で難しいと感じました。

絶対湿度や相対湿度のように逐次計算できないようでした。

hotRoomチュートリアルに湿度場を実装した演習は

成功しましたが、相対湿度の関数をどのファイルに記載すれば

良いかが分かりませんでした。

[初心者]

PMVを算出するFORTRANコードを下記のとおり作成し計算してみました。

この内容をhotRoomのチュートリアルで実装してみたいと思っているのですが、

0,const,systemのどのフォルダーのどのファイルにどのように記載して実施

すれば良いのか分かりませんでした。

ご教授頂けると有難いです。

　　　　　　　　　　　　　　　　　　　記

      REAL FAI,V0,RM,RMET,CLO,TA,TMRT,XA,FCL,PMV,PPD
C
      open(2,file='result.csv')
C
      WRITE(*,'(A\)') '  室内気温(℃)？===＞'
      READ(*,*) TA
C
      WRITE(*,'(A\)') '  室内相対湿度(％)？===＞'
      READ(*,*) FAI
C
      WRITE(*,'(A\)') '  室内気流速度(m/s)？===＞'
      READ(*,*) V0
C
      WRITE(*,'(A\)') '  代謝量(met)？===＞'
      READ(*,*) RMET

C///////////////////////////////////////
C
C     代謝量　1met＝50(kcal/m2h)
C
C///////////////////////////////////////
C
      WRITE(*,'(A\)') '  着衣量(clo)？===＞'
      READ(*,*) CLO
C
      WRITE(*,'(A\)') '  平均放射温度(℃)？===＞'
      READ(*,*) TMRT
 
      RM=50.0*RMET     

C
      CALL XS(TA,FAI,XA)

      CALL CALPMV(TA,XA,TMRT,V0,RM,0.0,CLO,FCL,PMV)

      PPD=100.0-95.0*EXP(-0.0353*PMV**4-0.2179*PMV**2)

      WRITE(*,601) PMV,PPD
      WRITE(2,601) PMV,PPD

  601 FORMAT(' ＰＭＶ=',F6.3,' ＰＰＤ=',F5.1)
C
      close(2)

      STOP

      END

C-----------------------------------------------------
C
C     PREDICTED MEAN VOTE CALCULATION SUBROUTINE
C
C-----------------------------------------------------

C////////////////////////////////////////////////////////
      SUBROUTINE CALPMV(TA,XA,RMRT,V,RM,W,CLO,FCL,PMV)
C////////////////////////////////////////////////////////
C
      DATA E,FEFF/0.97,0.71/
C
      FCL=1.0+0.2*CLO
      IF( CLO.GT.0.5 ) FCL=1.05+0.1*CLO
C
      TSK=35.7-0.032*(RM-W)
      PSK=1.92*TSK-25.3
      PA=760.0*XA/(0.622+XA)
      IC=0
      TCL=0.5*(TSK+TA)
      HCO=10.4*SQRT(V)
      R=4.96E-8*E*FEFF*FCL
      IF ( CLO.EQ.0.0 ) THEN
        TCL=35.7-0.032*(RM-W)
        GOTO 2
      END IF
C
      CLO1=1.0/(0.18*CLO)
C
    1 IC=IC+1
      HC=2.05*(ABS(TCL-TA))**0.25
      HC=AMAX1(HC,HCO)
      F=(TSK-TCL)*CLO1-R*((TCL+273.0)**4-(RMRT+273.0)**4)
     +                -FCL*HC*(TCL-TA)
      DF=-CLO1-R*4.0*(TCL+273.0)**3-FCL*HC
      DTCL=-F/DF
      TCL=TCL+DTCL
      IF(ABS(DTCL).LT.0.001) GOTO 2
      IF(IC.LE.20) GOTO 1
      WRITE(6,600)
  600 FORMAT(1H ,'***** PMV(TCL) DOES NOT CONVERGE *****')
      STOP
C
    2 CONTINUE
      RL=RM-W-0.35*(PSK-PA)-0.42*(RM-W-50.0)-0.0023*RM*(44.0-PA)
     +    -0.0014*RM*(34.0-TA)-R*((TCL+273.0)**4-(RMRT+273.0)**4)
     +    -FCL*HC*(TCL-TA)
      PMV=(0.352*EXP(-0.042*RM)+0.032)*RL
      RETURN
      END

C-----------------------------------------------------
C
C     SUBROUTINE FOR CALCULATION OF ABSOLUTE HUMIDTY
C
C-----------------------------------------------------

C////////////////////////////////////////////////////////
      SUBROUTINE XS(TAIR,FAI,XA)
C////////////////////////////////////////////////////////

      T=TAIR+273.15
      P1=11.344*(1.0-T/373.15)
      P2=-3.49149*(373.15/T-1.0)
      AA=ALOG(10.0)
      E1=EXP(P1*AA)
      E2=EXP(P2*AA)
      A=5.02808*ALOG10(373.15/T)-7.90298*(373.15/T-1.0)
     + -1.3816E-7*(E1-1.0)+0.0081328*(E2-1.0)
      PT=760.0*EXP(A*AA)
      XA=0.622*PT*FAI/(76000.0-PT*FAI)
      RETURN
      END

【計算事例】

①入力

  室内気温(℃)？===＞ 20.
  室内相対湿度(％)？===＞ 50.
  室内気流速度(m/s)？===＞ 0.1
  代謝量(met)？===＞ 1.1
  着衣量(clo)？===＞ 0.5
  平均放射温度(℃)？===＞ 23.49

②結果

 ＰＭＶ=-1.323 ＰＰＤ= 41.8

[初心者]

補足です。

室温はT,気流はＵ、圧力はp_rgh、平均放射温度は便宜的にT,

相対湿度は絶対湿度を実装して、RHを算出した結果を使用すればと

考えています。

2015年7月18日土曜日 23時20分32秒 UTC+9 初心者:

[jimi3776]

初心者様

初心者様の希望されている「(OpenFOAMへの)実装」は現実的に見て得策ではないと感じます。
なぜなら、OpenFOAMは名「Field Operation And Manipulation」のごとく、空間分布場を計算するので
初心者様の求めたいPMVの「空間分布場」を求めたければ、代謝量や着衣料などの空間分布場を0ディレクトリに
定義してあげることになります。
ここで、恐らく代謝量や着衣料などは空間分布ではなく、点分布でもよいのではないかと想像します。
この場合には、気温や気流速度などをいくつかのサンプリング点で抽出して、その後ポスト処理として、
その点毎にPVMの収束計算を行う方が得策に思います。
（少なくとも冒頭資料中のRH分布などはサンプリング点14点についての分布と思われます。）
そのサンプリング点の抽出方法は以下の資料が役に立つと思われます。
http://www.opencae.jp/data/OpenFOAMLectureCourses/200905/20090513_imano_sampleAndPlot_v2.pdf

さて、PMV計算に必要な「平均放射温度」とあるのですが、これは温度ソースにあたるわけですよね？
この部分が冒頭資料のソース項にあたるのではないでしょうか？
※輻射問題のように貴殿の問題のソース項が室内気温Tの関数ならば全く別問題です。

このソース項の実装および設定について、冒頭資料は２つの方法を提示していることに注意してください。
→その区別（＝資料p12と19の区別）ができるかが重要だと思います。
　個人的には、前者の方が初心者向けだが、細かい設定は煩雑？？(＝setFieldsの使い方次第)
　後者は、cellSet及びcellZoneの使い方次第でかなり自由度の高い設定が可能。（＝topoSetの使い方次第）
　因に、後者の実装は2.2.x以降でfvOptionsとして実装済みだと思います。

以上、PMVがどういう指標なのかも理解してないCFD寄りの者の意見でした。

[初心者]

jimi3776  様

ご回答有難うございます。

PMVの代謝量は事務作業で1met、着衣量はクールビズより0.5cloとして

取り敢えず、定数として与えたいと思うのですが、どこにどのように設定して

良いか分かりませんでした。

また、平均放射温度は本来ならグローブ温度を使うのですが、最初は

簡便のために温度Ｔと同じ値としたいと思っています。

本当は壁面からの対流熱伝達率から計算するのが良いのでしょうが・・・。

よって、ＰＭＶの６つの指標のうち、気温・気流はＴ，Ｕファイルがあります。

また、相対湿度は絶対湿度ソルバーをカスタマイズして算出可能かと思います。

残りの３つの指標は上述のとおり２つは定終、１つはＴと同じです。

ＰＭＶファイルも同じ計算ソルバーの中で出来るのではと思っているのですが、

失敗ばかり続きます。

ポスト処理として初心者の浅知恵で、U,T,RHファイルから計算して

最終ステップのフォルダーにＰＭＶを格納してparaviewで読み込ませようと

しましたがうまく行きませんでした。

ソース項の実装は、hotRoomのtutorialでsetFieldsの使い方は

簡単な事例は理解出来ました。

教えて頂いたファイルは勉強させて頂きます。

出来ればポスト処理ではなく、計算で他のファイルと同様に出力させたいと

考えています。

2015年8月1日土曜日 12時57分48秒 UTC+9 jimi3776:

[jimi3776]

paraviewなどで空間分布としてPMVを見たいわけですね。(快適空間はどこか？とか…)
ただ、初心者様の「計算中に同時に求める」ためにはかなりのC++及びOpenFOAMコードに精通している必要があると思います。
具体的にいえば以下のディレクトリにあるソースをほぼ理解できる必要があります。→同等なものを実装することになります。
https://github.com/OpenFOAM/OpenFOAM-2.3.x/tree/master/src/postProcessing

個人的に一番の近道と思われる方法は、やはり初心者様が試みられたT,U,RHファイル(ASCII形式)をポスト処理する方法だと思います。
大事に思われるのは、keyword internalField、boundaryFieldで、後者は
type zeroGradient;
にして、前者を
ununiform 〜〜;
と記述していくことだと思います。
計算結果のTファイルを眺めていれば、ルールがわかると思います。

[jimi3776]

PS
PMVファイルを最終ステップだけ置くだけではparaviewが変数として認識しないと思います。他のディレクトリ、特に0/にもダミーのPMVファイルを置いてください。

[初心者]

jimi3776  様

色々とアドバイス有難うございます。

ご紹介頂いたDIRソースはやはり難しいと感じました。

少しポスト処理を勉強して再挑戦してみたいと思います。

別添は簡単な箱モデルによるT,U,RH（humidity）を計算した

結果です。

このファイルを使ってPMVファイルを作成したいと思っているのですが

上手く行きませんでした。

1met、0.5clo、MRT=Tと仮定して、収束計算はニュートンラプソン法で

ＰＭＶを別途計算してから、PMVファイルとして計算フォルダーに

戻してからparaviewで見ようとしたのですが失敗ばかりでした。

[jimi3776]

読み込みに失敗した、24000要素の正しいPMV分布のASCIIファイルも添付してもらえれば、差し当たって達成したいことは解決すると思います。

[初心者]

jimi3776　様

お世話になっております。

説明不足でしたが、PMV失敗は相対湿度を算出するコードを追加して、それに

ＰＭＶもコードから一緒に計算しようとするところで失敗していました。

wmakeする方法でした。

ポスト処理する場合は、途中のboundary条件がU,T等で存在したりなかったりと

意味が分からなく途中で止めていました。

最初の数値の羅列は、ポイントの値だと分かり、全てこれだけであれば計算可能と

思ったのですが。

[jimi3776]

boundaryFieldは偏微分方程式を解く場合には非常に重要ですが、ポスト処理で求めるスカラー変数については
あまり重要ではない場合が多いと思います。
→PMVも空間分布が重要であって壁などの値はそれほど重要ではないと感じます。

T,U,humidityのboudaryFieldは無視して、internalFieldだけ読み込んで
以下のようなPMVファイルを作ればparaviewに読み込めると思います。
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "16000";
object PMV;
}
dimensions [0 0 0 0 0 0 0];
internalField nonuniform List<scalar>
24000
(
〜
〜
)
;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
symmetry
{
type symmetryPlane;
}
top
{
type zeroGradient;
}
ground
{
type zeroGradient;
}
wall
{
type zeroGradient;
}
}

[初心者]

jimi3776　様

有難うございます。

挑戦してみます。

16000stepで計算が終了しているので、そのフォルダーにPMVファイルを作成してみます。

0フォルダーにも何か必要でしょうか。

T,Uなどはファイルがあるので必要のような感じがしますが。

また、Uファイルだけベクトルのため(x,y,z)が出力されています。

とりあえずこれを平均してスカラーの１つの値だけにしてからPMVを別途計算して

ファイルを作成してみます。

[jimi3776]

FoamFile 
{ 
    version     2.0; 
    format      ascii; 
    class       volScalarField; 

    location    "0"; 

    object      PMV; 
} 
dimensions      [0 0 0 0 0 0 0]; 

internalField   uniform 0;

[初心者]

jimi3776　様


別添のとおり実施してみました。

最終ステップの16000フォルダーにＰＭＶを格納しました。
期は初期ファイルの0フォルダー「internalField」
初ＰＭＶ＝０の空間として解析してみました。   

internalField   uniform 0;

別添pngファイルが中央断面です。
左上天井から１０℃の冷たい冷房をしているので、ＰＭＶも－２程度と
少し寒くなっており、右側は少し暖かくなっており、それなりに合っている
ようにも感じます。この方法で良かったのでしょうか。

また、速度ベクトルをＸ，Ｙ，Ｚ全て足して３で割った平均値を速度スカラーとして
計算しましたが良かったでしょうか。

便宜的にＭＲＴ＝Ｔとしていること、また、ニュートンラプソン法では収束
しないので２分法で計算しました。

ポスト処理してから戻すという方法しか初心者レベルでは出来ませんでしたが、
大変勉強になりました。ありがとうございました。

今後は今野様やjimi3776様のように、一貫してソースカスタマイズから
ＰＭＶが計算できるといいのですが、まだ無理そうです。
簡単な箱モデルで試してみましたが、定数項や収束計算するなど、
似たようなケースファイルなどがあると勉強になります。

同様にＰＰＤも計算してみたところ上手く行きました。
ＰＭＶ、ＰＰＤも作成出来て勉強になりました。

[jimi3776]

0ディレクトリに何が必要か分かれば以下のページで何をやってるか想像つくと思います。
http://www.geocities.jp/penguinitis2002/study/OpenFOAM/tankentai/07-sfield.html
頑張って下さい。
PS
U=sqrt(Ux**2 +Uy**2+Uz**2) by fortran
でしょうね。

[初心者]

jimi3776  様

アドバイスありがとうございました。
速度ベクトルのスカラー変換は訂正してみます。

今回はポスト処理の方法で解決する方法を知れて勉強になりました。
今後、似たような解析に応用できそうです。
ありがとうございました。

[初心者]

jimi3776　様

速度Ｕファイルを、

U=sqrt(Ux**2 +Uy**2+Uz**2)

として再計算しました。別添のとおりです。

ＰＭＶが－３程度も存在し、結構、不快な環境でした。
ＰＰＤも計算しましたが、99％も不快という場所もありました。
冷房時なので本来なら＋値が適当だと思いますが冷え過ぎています。

定性的・定量的にも吹き出し口近くは寒く不快であり、
離れると改善されてきますが、それでも良い値ではありません。

今回はＰＭＶを可視化する目的でしたので、部屋の大きさに
対する設備容量検討はしていません。

ＰＭＶが0.5程度になるように最適な設備設計の検討が本来は必要です。

また、平均放射温度ＭＲＴも温度Ｔと同じとして計算しています。
本来だと、この部分の計算が必要ですが、壁面表面温度や形態係数など、
結構大変そうでした。

ＭＲＴ計算には形態係数の計算が必要で、今回の箱モデルでも

24000個もデーターがあり、しかも何処が壁面か分かりません。

形態係数の計算も大変そうです。

ＯＦで計算出来るのかも分かりません。

速度スカラー変換を間違えており正しく実施したところ、前回より
悪い結果になり、またＰＰＤも追加しましたので参考までに添付します。

[Masashi Imano]

今野です．

> また、平均放射温度ＭＲＴも温度Ｔと同じとして計算しています。
> 本来だと、この部分の計算が必要ですが、壁面表面温度や形態係数など、
> 結構大変そうでした。

MRTを気温Tで代用するには，周囲からの放射熱伝達を全く考慮しないことになるので，

本来不適当です．

ただし，MRTを形態係数を用いて算出するのは，確かに難しいですし，計算時間もかかります．

格子の中心に人体をモデル化した微小立方体もしくは微小球体を配置し，

それらとパッチ間の形態係数を，全格子にわたって算出する必要があるからです．

参考文献

http://repo.flib.u-fukui.ac.jp/dspace/bitstream/10098/3403/1/AN00215401-046-02-001.pdf

そこで，室内のように閉空間の場合には，全パッチの面積重み付け平均温度を

MRTとする簡便な方法が用いられることもあります．

全パッチの面積重み付け平均温度を求めるプログラムを作成するのはそれほど難しくないですが，

全パッチにおける温度TについてpatchIntegrateを実行し，パッチの面積と，温度との積算値を表示

させる方法でも求めることができます．

patchIntegrate T パッチ名

Time = 時刻

    Area magnitude of patch パッチ名[0] = 31.844 ← パッチの面積 [m2]

    Integral of T over area magnitude of patch パッチ名[0] = 9730.28 ← パッチの界面の面積と温度の積算値l [m2 K]

2015年8月4日火曜日 20時39分46秒 UTC+9 初心者:

[初心者]

今野　様

アドバイス有難うございます。

その後、私もMRT算出方法をネットで調べていたところ、

今野様が指摘したように、形態係数算出は次の通り

レイトレーシング法など、かなり高度で大変な事が分かりました。

http://ir.library.osaka-u.ac.jp/dspace/bitstream/11094/27633/1/24598_%E8%AB%96%E6%96%87.pdf

他の簡易な方法は、今野様のアドバイス通りに面積加重平均の考え方もあることを知りました。

こちらは位置に寄らずに壁表面温度と面積だけで簡易に求まる事も知りました。

全パッチにおける温度TについてpatchIntegrateを実行し，パッチの面積と，温度との積算値を表示

させる方法でも求めることができるという理屈は、上述の面積加重計算で概念は理解できるのですが、

OpenFOAMで実施しようとすると、どのファイルから面積と温度を読み取って、コマンドでどのように

命令すれば良いか分かりませんでした。

patch関連のファイルはconstant/polyMeshファイルだと思われ、別添のとおり添付しますが、

その後の実施方法が分かりませんでした。

具体的な手順等を教えて頂くと有難いです。

コマンド命令はblockMesh、snappyHexMesh、appplication名しか実施したことがありません。

[Masashi Imano]

今野です．

constant/polyMesh/boundary の中を直接見るか，patchSummary コマンドで

パッチの情報を出力してください．

$ patchSummary -expand -latestTime

なお，パッチ界面もしくは，グルーピングされた界面間の形態係数を求めるだけなら，

OpenFOAM のviewFactorsGen で行えますが，一般に精度が悪い二重面積分法で算出されています．

形態係数の算出法の比較については，拙著も参考にしてください．

今野 雅,倉渕 隆,鎌田 元康 放射熱伝達解析における形態係数の計算法についての研究
http://ci.nii.ac.jp/naid/110004088161

2015年8月5日水曜日 20時17分16秒 UTC+9 初心者:

[初心者]

今野　様

有難うございます。

初めてのコマンドなので一つずつ試してみます。

少しお時間をください。

今野様の論文を拝見しようとアクセスしたところ、

日本建築学会の論文で会員になってないと閲覧出来ない内容となっていました。

まずは、patchSummary コマンドを試してみます。

初めて聞くコマンドばかりで少し勉強させてください。

[初心者]

今野　様

早速、「patchIntegrate T パッチ名」を実施してみました。

①top

Create time

Create mesh for time = 0

Time = 0
    Area vector of patch top[3] = (0 7.92 0)
    Area magnitude of patch top[3] = 7.92
    Reading volScalarField T
    Integral of T over vector area of patch top[3] = (0 2376 0)
    Integral of T over area magnitude of patch top[3] = 2376

Time = 16000
    Area vector of patch top[3] = (0 7.92 0)
    Area magnitude of patch top[3] = 7.92
    Reading volScalarField T
    Integral of T over vector area of patch top[3] = (0 2337.04 0)
    Integral of T over area magnitude of patch top[3] = 2337.04

End

②ground

Create time

Create mesh for time = 0

Time = 0
    Area vector of patch ground[4] = (0 -8 0)
    Area magnitude of patch ground[4] = 8
    Reading volScalarField T
    Integral of T over vector area of patch ground[4] = (0 -2480 0)
    Integral of T over area magnitude of patch ground[4] = 2480

Time = 16000
    Area vector of patch ground[4] = (0 -8 0)
    Area magnitude of patch ground[4] = 8
    Reading volScalarField T
    Integral of T over vector area of patch ground[4] = (0 -2442.83 0)
    Integral of T over area magnitude of patch ground[4] = 2442.83

End

③wall

Create time

Create mesh for time = 0

Time = 0
    Area vector of patch wall[5] = (0.08 0 -12)
    Area magnitude of patch wall[5] = 23.92
    Reading volScalarField T
    Integral of T over vector area of patch wall[5] = (24 0 -3600)
    Integral of T over area magnitude of patch wall[5] = 7176

Time = 16000
    Area vector of patch wall[5] = (0.08 0 -12)
    Area magnitude of patch wall[5] = 23.92
    Reading volScalarField T
    Integral of T over vector area of patch wall[5] = (17.5961 0 -3539.53)
    Integral of T over area magnitude of patch wall[5] = 7054.53

End

④inlet

Create time

Create mesh for time = 0

Time = 0
    Area vector of patch inlet[0] = (0 0.08 0)
    Area magnitude of patch inlet[0] = 0.08
    Reading volScalarField T
    Integral of T over vector area of patch inlet[0] = (0 23.44 0)
    Integral of T over area magnitude of patch inlet[0] = 23.44

Time = 16000
    Area vector of patch inlet[0] = (0 0.08 0)
    Area magnitude of patch inlet[0] = 0.08
    Reading volScalarField T
    Integral of T over vector area of patch inlet[0] = (0 23.44 0)
    Integral of T over area magnitude of patch inlet[0] = 23.44

End

⑤outlet

Create time

Create mesh for time = 0

Time = 0
    Area vector of patch outlet[1] = (-0.08 0 0)
    Area magnitude of patch outlet[1] = 0.08
    Reading volScalarField T
    Integral of T over vector area of patch outlet[1] = (-24 0 0)
    Integral of T over area magnitude of patch outlet[1] = 24

Time = 16000
    Area vector of patch outlet[1] = (-0.08 0 0)
    Area magnitude of patch outlet[1] = 0.08
    Reading volScalarField T
    Integral of T over vector area of patch outlet[1] = (-23.6706 0 0)
    Integral of T over area magnitude of patch outlet[1] = 23.6706

End

        S(m2)   S*A(m2K） 

inlet      0.08    23.676 
outlet     0.08   23.44 
top        7.92    2337.04 
ground    8       2442.83 
wall       23.92   7054.53 

合計        40      11881.516 
   
   
 MRT     297.0379 (K)
   
              23.8879 (℃）

このように、荷重平均によるMRTは23.8879℃として良いという事でしょうか。

また、次のコマンドは以下のようにエラーとなり原因が分かりませんでした。

 patchSummary -expand -latestTime

Usage: patchSummary [OPTIONS]
options:
  -case <dir>       specify alternate case directory, default is the cwd
  -constant         include the 'constant/' dir in the times list
  -latestTime       select the latest time
  -noFunctionObjects
                    do not execute functionObjects
  -noZero           exclude the '0/' dir from the times list, has precedence
                    over the -zeroTime option
  -parallel         run in parallel
  -region <name>    specify alternative mesh region
  -roots <(dir1 .. dirN)>
                    slave root directories for distributed running
  -time <ranges>    comma-separated time ranges - eg, ':10,20,40-70,1000:'
  -srcDoc           display source code in browser
  -doc              display application documentation in browser
  -help             print the usage

Using: OpenFOAM-2.0.x (see www.OpenFOAM.com)
Build: 2.0.x-68952d880745

--> FOAM FATAL ERROR:
Invalid option: -expand

FOAM exiting

boundaryファイルの中身は次の通りです。

FoamFile
{
    version     2.0;
    format      ascii;

    class       polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

6
(
    inlet
    {
        type            patch;
        nFaces          8;
        startFace       69400;
    }
    outlet
    {
        type            patch;
        nFaces          8;
        startFace       69408;
    }
    symmetry
    {
        type            symmetryPlane;
        nFaces          1200;
        startFace       69416;
    }
    top
    {
        type            wall;
        nFaces          792;
        startFace       70616;
    }
    ground
    {
        type            wall;
        nFaces          800;
        startFace       71408;
    }
    wall
    {
        type            wall;
        nFaces          2392;
        startFace       72208;
    }
)

[初心者]

今野　様

ＭＲＴを加重平均で算出し、met,clo,MRTを定数としてPMV,PPDを計算してみました。
別添のとおりです。

前回のMRT＝Tという結果と違い、快適環境は改善された結果となりました。
かなり数値が変わると感じました。

今野様の論文やレイトレーシング法などの高度な解析は少し難しそうですが、
この方法であれば位置に寄らずに簡易にＭＲＴを算出できて便利そうでした。

取り急ぎ検討した結果を報告致します。

ありがとうございました。