OpenModelicaのBuildingsライブラリで換気回路網計算を学ぶ_その3

概要

OpenModelicaでLBNLのBuildingsライブラリを使う。初心者なのでいろいろ教えて頂けるとありがたい。
前回Buildings.Airflow.Multizoneパッケージで風力換気の計算をしてみた。今回はその続きで気象データを読み込んで風向と風速が時間変化するときの換気量をみていきたい。

使用バージョン
-OpenModelica1.12.0 →Modelica標準ライブラリ3.2.2
-Buildings 5.1.0

作成するモデルの概要

前回作った風力換気の単室モデルに適当なスケジュールで発熱を与えて気象データを読み取って自然室温のようなものを計算する。外壁の熱貫流は無視しているので自然室温ではない。

べき乗則モデルの作成

前回と同様にBuildings.Airflow.Multizone.Examples.OneRoomのモデルをベースにするのでコピーしてくる。今回は外気側のボリュームを外部風向に合わせて風圧係数を計算できるBuildings.Fluid.Sources.Outside_CpLowRiseにしたいのだが、そのまま使うのも味気ないのでべき乗則で風速を補正して使うようにしてみる。気象データで得られる風速の計測高さと風圧係数を求めたときの風速の高さ（軒や屋根の高さ）は必ずしも一致しないので、気象データの風速をべき乗則などで補正してやる必要がある。気圧も高さが変わると補正が必要だと思うけれど、そもそも今は定数として標準大気圧を使っているので特に補正はしない。
今回使う気象データ（USA_IL_Chicago-OHare.Intl.AP.725300_TMY3.mos）は風速計高さを調べても不明だったので30mを6mに換算してみる。energyplusの気象データの風速計の高さをどこで見るのかご存知の方がいれば教えて頂けるとありがたい。TMY3のマニュアルをみたりしたけれどよくわからず観測所の情報も緯度経度と高度はわかっても風速計高さにたどりつけなかった。アメダスの場合は気象庁の地域気象観測所一覧から確認できる。拡張アメダス標準気象データは風速計の高さを補正している場合があるので補正時と異なるべき指数を使いたい場合は注意が必要である。それから、元の風速計の高さが低くい場合は周辺の建物の影響を受けていたり、十分な高さの気象官署からデータを取るにしても検討地点から距離が離れすぎててたりと悩ましいことが少なくないのだけれどそこは割り切るしかないかもしれない。（よい検討方法が実はすでにあるかもしれないけれど。）

まずは前回のように修正したいモデルOutside_CpLowRiseをコピーしてくる。
変数として気象データの風速計高さhWea、風圧算出用の風速高さhCor、補正用のべき指数nCor、風圧算出用の風速vCorを加えてべき乗則で風速を補正して使用するように書き換える。元のモデルでは気象データから読み込む風速はvWinという変数に入れられている。

model Outside_CpLowRise "Boundary that takes weather data as an input and computes wind pressure for low-rise buildings"
  extends Buildings.Fluid.Sources.BaseClasses.Outside;
  parameter Modelica.SIunits.Height hWea "height in weathedata";  //追加部分
  parameter Modelica.SIunits.Height hCor "height for Wind pressure coefficient";  //追加部分
  parameter Real nCor "coefficient for power low";  //追加部分
  Modelica.SIunits.Velocity vCor "Wind velocity at hCor";  //追加部分
(略)
equation
  alpha = winDir - surOut;
  CpAct = Buildings.Airflow.Multizone.BaseClasses.windPressureLowRise(Cp0 = Cp0, incAng = alpha, G = G);
  vCor = vWin * ( ( hCor / hWea ) ^ nCor ); //追加部分
//  pWin = 0.5 * CpAct * medium.d * vWin * vWin;  //修正前コメントアウト
  pWin = 0.5 * CpAct * medium.d * vCor * vCor;  //修正後
  pTot = pWea + pWin; 
(略)
end Outside_CpLowRise;

今回はこのような変更にしたけれど、connect(weaBus.winSpe, vWin);の部分を消して、別途べき乗則用のモデルを用意して補正した風速をvWinにつないでもよいかもしれない。

できたら左のライブラリブラウザからドラッグ＆ドロップしてきてoutNとoutSという名前をつけてやる。右クリックで横反転してやると見やすい。ここでは北壁を上の窓につなぐようにした。 ダブルクリックするとちゃんと設定にhWea、hCor、nCorが追加されているので30、6、0.28としてやる。べき指数は地表面の状態で変わる指数で、今回は平屋建物が立ち並ぶ郊外を想定して0.28にした。

風圧係数の関数について

風圧係数についてはOutside_CpLowRiseに組み込んである関数を使って算出するのだが。この関数はBuildings.Airflow.Multizone.BaseClasses.windPressureLowRiseにある関数で、BaseClassesと同じ場所にあるがfunctionである。 これを使わない場合の風圧係数は何らかの設計資料から似た形のものを引用したり、CFDから風圧係数を計算したり、場合によっては風洞実験をするなどの方法で風向別のデータを得ることができる。その場合は別途風向から風圧係数を算出するモデルを作るか、風向から風圧係数を事前に求めておいて時系列データとして読み込んでBuildings.Fluid.Sources.Outside_Cpに入力することになる。
今回使う気象データはチュートリアルのSpaceCooolingと同じものでは北風が0°で時計回りに東風が90°で南風が180°西風が270°となっており10°刻みで風向が定義されている。日本だと16風向で22.5°刻みのことが多いので不思議な感じがする。無風時はデータが空になっているらしいが今回はそのようなデータはなかった。ちなみに一般に風向の定義として南風は南から北へ向かって吹く風である。逆と間違えないようにしたい。

WindPressureLowRiseの関数はSwami and Chandra (1987) の論文によるものだ。 Cp0は風向が壁に正対するとき（南向きの壁に南風が吹いたとき）の風圧係数で一般的な低層の建物では0.6が推奨らしい。 sは（対象の壁の幅/対象の壁に隣接する壁の幅）である。隣接する壁の幅とは対象の壁に正対したときの建物の奥行き方向の長さになる。ここでは幅18m奥行き22mの建物と考えて0.81と入力してみる。 aziは壁の向きで南が0°、西が90°になっている。角度の単位はdegreeかradianか選べる他、東西南北についてはBuildings.Types.Azimuth.Nなどで指定してやることができる。これは次のように定義されているので単位はradianである。

within Buildings.Types;
package Azimuth "List of possible constant values for surface azimuth"
   constant Modelica.SIunits.Angle E = -Modelica.Constants.pi/2
   constant Modelica.SIunits.Angle N = Modelica.Constants.pi
   constant Modelica.SIunits.Angle S = 0
   constant Modelica.SIunits.Angle W = +Modelica.Constants.pi/2
end Azimuth;

これでパラメータの設定ができた。nportは流体の入出力ポートでこれを1にしておく。CO2とか粉塵とかの計算はしないのでuse_C_inはfalseにする。それから、テキストビューで媒体の指定を追加してやる。

OneRoomTest.Outside_CpLowRise outN(redeclare package Medium = Medium,azi(displayUnit = "rad") = Buildings.Types.Azimuth.N,hCor = 6, hWea = 30, nCor = 0.28, nPorts = 1, s = 0.81)

気象データと発熱スケジュールの読み込み

気象データはBuildings.BoundaryConditions.WeatherData.ReaderTMY3を使う。湿球温度は特にいらないのでcomputeWetBulbTemperatureをfalseにし、filNamはライブラリにはじめからついているUSA_IL_Chicago-OHare.Intl.AP.725300_TMY3.mosを使う。他のパラメータはデフォルトのままで大気圧はパラメーター固定値で他はファイルから、天空温度も使わないのでそのままにしている。

さらに換気による排熱の効果も見たいと思ったので室の内部発熱スケジュールを組み込んでみたいと思う。これはModelica標準ライブラリを使う。発熱スケジュールデータをModelica.Blocks.Sources.CombiTimeTable から読み込み、Modelica.Thermal.HeatTransfer.Sources.PrescribedHeatFlowを経由してheatportとしてvolEasにつなぎこんでやる。PrescribedHeatFlowはデフォルトの設定のままでよいがインプットをそのまま出力するだけなのでalphaが0であることを確認する。
CombiTimeTableによる時系列データの入力は3種類あり、直接記述、ファイル読み込み、cのソースに書き込みができるようだがcのソースに書き込みの使い方がいまいちよくわからなかった。ファイルの読み込みはASCIIファイルのフォーマット以外にMATLABのフォーマットを使うことができるらしい。OMEdit上の設定は次のとおり。

tableOnFileをfalseにしたときはtableのところに直接テーブルを記述する。[0,0;1,2;2,0;3,2]のように各時刻のデータをカンマで区切り、時刻が変わるときにはセミコロンで区切る。毎時刻の一番初めのデータは時刻（秒）になるようにして時刻は単調増加（減少しない事）にしないといけない。基本的にテーブルに無い時刻のデータは線形補間で計算される。不連続に値を変動させたい場合は[0,0;1,0;1,2;2,0;3,2]などのようにすると1秒のところで不連続になる。この時刻では計算が2回行われるのだと思う。
tableOnFileをtrueにした場合はfileNameが”NoName”だとcのソースに書き込みになるらしい。"C:\OpenModelica1.12.0-64bit\lib\omlibrary\Modelica 3.2.2\Resources\Data\Tables"とかにあるusertab.cやらusertab.hやらにサンプルらしきものがある。
そしてファイル名を指定してやるとファイルの読み込みになる。同じところにtest.txtというサンプルファイルがある。読み込みファイルの形式としては
* 1行目に#1と記載する。（データ形式のバージョンを示す） * double tab0(2,2)のようにデータ型（floatかdouble）、テーブル名（ここではtab0）、データの数（行数、列数）を示す。 * テーブルのデータをスペース、タブ、カンマ、セミコロンのいずれかで区切って記載していく。一列目には単位を秒として時刻を記述する。改行コードは特にその行のデータの終わりを意味しないので複数行に同じ時刻のデータを記載できる。 * 一つのファイルに複数のテーブルを記載できる。tableNameで読み込むテーブルを選択する。

今回はファイルを読み込んで使用するのでtableOnFileをtrueにし、tableはそのまま、tableNameは"heat"にして（ダブルクオーテーションを忘れないように）、fileNameは"C:/Users/Simulation/OModelica/heat.txt"とheat.txtを作成して指定する。verboseReadは読み込みのメッセージの出力設定なのでfalseでよいと思う。作成したファイルは次のとおり。9時から18時まで100Wの発熱を与えるようにする。時間が秒の単位になっていないがそれは別に設定する。テーブルが大きい場合はExcelでつくってcsvなどで書き出してやるとよいと思う。

#1
double heat(4,2)
0,0
9,100
18,0
24,0 

設定の続きにもどる。columnsは読み込む列の番号を指定するのだと思うけれど使い方があまりよくわからなかった。今回はたぶん2列のデータなのでデフォルトで特に問題なかった。
smoothnessはテーブルに記載の無い時刻でのデータの補間の方法を指定する。デフォルトはLinearSegmentsで線形補間だがConstantSegmentsだとその前の時刻の値を保持することになるのでこれを指定する。
extrapolationはテーブルの時刻の外側のデータの取り扱いでデフォルトはLastTwoPointsで線形外挿だがPeriodicとして周期的にテーブルを使用する。これで1日分のスケジュールを毎日繰り返すようになる。 offsetは時刻以外のすべてのデータに加える値で0のままに、startTimeはすべての時刻に加える値で0のままにする。
timeScaleはすべての時刻に掛ける値でこれを3600にする。こうすることでスケジュールデータの時刻の単位を秒に変換できる。これで設定ができた。

時系列の換気量の計算

設定したファイルをつなげてやる。weaBusなどはインスタンスの名前がついていないことがあるので注意。

今回追加した部分の接続は以下の通り。

  connect(colOutTop.port_b, outN.ports[1]);
  connect(colOutBot.port_a, outS.ports[1]);
  connect(weaDat.weaBus, outN.weaBus);
  connect(weaDat.weaBus, outS.weaBus);
  connect(combiTimeTable1.y[1], prescribedHeatFlow1.Q_flow);
  connect(prescribedHeatFlow1.port, volEas.heatPort);

計算を開始する。年間計算として終了時刻を31532400s、計算は間隔3600sにしてみる。年間計算したけれどそこまで計算時間が気にはならなかった。そして年間計算したけれど結局集計や詳細にデータを見るのがたいへんなので部分だけ取り出して計算ができてそうだというのを確認するだけになった。
室温の変化は以下の通りで昼間発熱して温度があがって夜間に換気で外気近くまで温度が下がる様子が計算できた。ただ熱容量や壁の伝熱がないと100Wでこんなに室温があがるものなのかと驚いた。(容積は2.5×5×5m3)

計算開始後の換気量（上）と風速（下）の関係は以下の通りで風速の大きいところでしっかり換気量が増えているのがわかる。風向毎に風圧係数がかわったり浮力（温度差）の影響があったりするので必ずしもこんなにきれいに相関は出ないのが普通かもしれないけど今回はきれいにでた。

csvで出力すると発熱開始と終了のたびにデータが3つでていた。データの集計の際に面倒なことになるので何かよい方法は無いものかなぁと思った。

おわりに

とりあえず時系列の換気計算ができた。理解が間違っている部分があれば指摘していただけるとありがたい。
風圧係数については8月に建築研究所が建築物の自然換気設計のための風圧係数データベースを公開した。まだ読めていないのだけれど、単なるデータベースではなく風圧係数についての知識や考察がたっぷり書いてあるようなのでとても勉強になると思う。このデータベースは32風向分用意されており、日本で得られる気象データは16風向なのでそのまま使える。ただし、今回使った気象データは10°刻みの36風向の気象データになっていたので、このデータベースを使うなら風圧係数を補間してやる必要がでてくる。（線形補間してよいのかは分からないのだけれど、分からなければとりあえず線形補間するのだろう）

次回は自然室温や熱負荷を計算するためのパッケージと思われるBuildings.ThermalZonesについて少し見ていきたいと思う。チュートリアルのSpaceCoolingでは定常計算の壁体熱貫流と換気負荷を見込んだ単純な室だったが、壁体を非定常熱伝導にしたり、窓を取り付けたり、室内の放射伝熱などの計算をしたりなど、より詳細な室の取り扱いができるといいなぁと思う（Buildings.HeatTransferに関連モデルがあるようだけれど室のサンプルはThermalZonesのはず）。次回見るサンプル(Buildings.ThermalZones.Detailed.Validation.BESTEST.Cases6xx.Case600FF)もOpenModelicaではエラーが出るようなのでJModelicaと併用していくと思う。