FMUのパラメータ推定ライブラリmodestpyを使う

概要

表題のとおりmodestpyはFMU(Functional Mockup Unit)のパラメータ推定をするpythonのライブラリで、ライセンスはBSD 2-clauseライセンスである。FMUの計算自体はJModelicaのpyFMIを使っている。ある時系列データに対してシミュレーション結果とのRMSEが小さくなるようなパラメータを探索するようになっている。資料はAPIのドキュメントとカンファレンスのプレプリントPDFがある。

使用バージョン
- modestpy 0.0.7
- JModelica2.4 →Modelica標準ライブラリ3.2.2

インストール

condaからインストールするのが推奨とサイトには書いてあるけれど、JModelicaのpython上で使えたほうがよさそうなのでコマンドプロンプトで以下のようにしてgithubの最新版のmedestpyをインストールした。

C:\JModelica.org-2.4\Python27\Python_64\python.exe -m pip install https://github.com/sdu-cfei/modest-py/archive/master.zip

そうするとC:\JModelica.org-2.4\Python27\Python_64\Lib\site-packages\modestpyにインストールされる。インストールしたらJModelicaのIPythonを起動してテストファイルを実行してみる。とりあえず何かしら計算をしているようだがよくわからない。

from modestpy.test import run
run.tests()

サンプルの実行

exampleフォルダがgithubにあるのでとりあえず作業フォルダにコピーしてくる。example\simpleの中にスクリプトがいくつかあり、resourcesフォルダの中にいくつか設定用のファイルがある。とりあえず計算に使うモデルSimple2R1C_ic.moをOpenModelicaで開いてみる。フォルダの中にはいくつかのプラットフォーム用のfmuファイルがあるものの、すべての元はこのファイルである。以下のようなモデルになっている。
f:id:kinonotofu:20181127221828p:plain

熱容量1000[J/W]の点に熱抵抗1[K/W]が二つつながっており、抵抗の反対側は片方がTi1[K]、もう片方がTi2[K]になっている。熱容量の点は初期温度Tstart=20で、温度Tはセンサーで取得している。
このモデルのパラメータを推定するらしい。モデルのなかでparameterとして定義されているのは二つの熱抵抗R1とR2、初期温度Tstart、熱容量Cの4つである。温度の時系列データを入力して、Tstart以外の3つのパラメータを推定するのだが、まずは広く探索するソルバーで計算して途中から局所的に探索していくらしい。よくわからないのでとりあえず実行してみる。
まずはexample\simple\simple.pyを作業フォルダにコピーしてくる。 simple.pyを開いてif __name__ == "__main__":コメントアウトし、そこから先を一括選択してShift+tabでインデントを元に戻す。(操作はエディタにもよるかもしれない)
あとはJModelicaのIPythonで以下のようにすれば実行される。(とくに編集せずにコマンドプロンプトからpython simple.pyでもよいと思う。)

cd '作業フォルダ'
import simple

なにやらエラーのような終わり方をしている。ただし計算はできていそうだった。
f:id:kinonotofu:20181127225250p:plain

計算結果は作業フォルダ/examples/simple/workdirに保存されるようになっている。
summary_1.csv(左)とsummary_2.csv(右)に以下のようにパラメタスタディのログが出力されている。GAが広域探索として遺伝的アルゴリズムを、PSが局所探索としてパターンサーチを使用している部分である。局所探索は他にもScipyにある探索が使えるらしい。
f:id:kinonotofu:20181127225846p:plain

best_per_run.csvにそれぞれのなかでもっともよいパラメータを書いている。
f:id:kinonotofu:20181127230035p:plain

そしてfinal.csvにもっともよい結果が出力される。
f:id:kinonotofu:20181127230128p:plain

それからいくつかの画像が出力される。
GAの結果(計算回数の枚数だけ出力)
f:id:kinonotofu:20181127230908p:plain

PSの結果(計算回数の枚数だけ出力)
f:id:kinonotofu:20181127230936p:plain

RMSEの推移 f:id:kinonotofu:20181127231041p:plain

推定したパラメータでの計算結果と推定に用いた実測値の比較
f:id:kinonotofu:20181127231132p:plain

とりあえずうまくパラメータを推定できたようである。

設定の確認

simple.pyは以下のようになっている。

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
import json
import os
import pandas as pd
from modestpy import Estimation
from modestpy.utilities.sysarch import get_sys_arch

#if __name__ == "__main__":
"""
This file is supposed to be run from the root directory.
Otherwise the paths have to be corrected.
"""
# DATA PREPARATION ==============================================
# Resources
platform = get_sys_arch()
assert platform, 'Unsupported platform type!'
fmu_file = 'Simple2R1C_ic_' + platform + '.fmu'
fmu_path = os.path.join('examples', 'simple', 'resources', fmu_file)
inp_path = os.path.join('examples', 'simple', 'resources', 'inputs.csv')
ideal_path = os.path.join('examples', 'simple', 'resources', 'result.csv')
est_path = os.path.join('examples', 'simple', 'resources', 'est.json')
known_path = os.path.join('examples', 'simple', 'resources', 'known.json')

# Working directory
workdir = os.path.join('examples', 'simple', 'workdir')
if not os.path.exists(workdir):
    os.mkdir(workdir)
    assert os.path.exists(workdir), "Work directory does not exist"
# Load inputs
inp = pd.read_csv(inp_path).set_index('time')
# Load measurements (ideal results)
ideal = pd.read_csv(ideal_path).set_index('time')
# Load definition of estimated parameters (name, initial value, bounds)
with open(est_path) as f:
    est = json.load(f)
# Load definition of known parameters (name, value)
with open(known_path) as f:
    known = json.load(f)

# MODEL IDENTIFICATION ==========================================
session = Estimation(workdir, fmu_path, inp, known, est, ideal,
                        lp_n=2, lp_len=50000, lp_frame=(0, 50000),
                        vp=(0, 50000), ic_param={'Tstart': 'T'},
                        methods=('GA', 'PS'),
                        ga_opts={'maxiter': 5, 'tol': 0.001, 'lhs': True},
                        ps_opts={'maxiter': 500, 'tol': 1e-6},
                        scipy_opts={},
                        ftype='RMSE', seed=1,
                        default_log=True, logfile='simple.log')
estimates = session.estimate()
err, res = session.validate()

はじめの方は読み込みファイルのパスを変数に入れたりしており、MODEL IDENTIFICATION 以降が計算の実行部となっている。
session = Estimation(~)で変数を一括で設定し、
estimates = session.estimate()でパラメータ推定、
err, res = session.validate()で推定した値を使用して実験値と計算値の比較をしている。

設定に使用するファイルは以下の通り。

  • fmu_file:FMUファイル、JModelicaなどで作成できる。
  • inputs.csv:Ti1とTi2の時系列データ
  • result.csv:inputs.csvの条件の時のTの実測データ
  • est.json:パラメータ探索の設定
  • known.json:Tの初期値Tstartの設定

パラメータ探索の設定est.jsonは以下のようになっている。だいたいの目安の値、最小値、最大値の順番に記載している。 "R1": [0.08, 0.001, 0.5],
"R2": [0.08, 0.001, 0.5],
"C": [1000.0, 100.0, 10000.0]

その他は以下のように設定している。

  • lp_n=2:探索の回数
  • lp_len=50000:1回の探索の長さ
  • lp_frame=(0, 50000):探索回数の最小値と最大値
  • vp=(0, 50000):探索後の検証の計算の長さ
  • ic_param={'Tstart': 'T'}:初期値の対応関係
  • methods=('GA', 'PS'):使用するソルバーの指定、現状GAは固定
  • ga_opts={'maxiter': 5, 'tol': 0.001, 'lhs': True}:GAのソルバーの設定、詳細はmodestpy/estim/ga/ga.pyのコメント参照
  • ps_opts={'maxiter': 500, 'tol': 1e-6}:PSのソルバーの設定、詳細はmodestpy/estim/ps/ps.pyのコメント参照
  • scipy_opts={}:Scipyのソルバーの設定、詳細はmodestpy/estim/scipy/scipy.pyのコメント参照
  • ftype='RMSE':目的関数現状はNRMSEかRMSEのみ。
  • seed=1:乱数使用時の種、指定しなければ現在時間などOS指定の値を使う。
  • default_log=True, logfile='simple.log':ログの設定

終わりに

今回は第6回Modelicaライブラリ勉強会でパラメタ最適化している話を聞いて、そういえば何かあったよなぁと思って調べてみた。まだ発展途上中のシンプルなライブラリだがしっかり計算はできるようなのでなかなかよさそうなライブラリだと思った。

ModelicaのBuildingsライブラリのHeatTransferを学ぶ_その1

概要

OpenModelicaでLBNLのBuildingsライブラリを使う。初心者なのでいろいろ教えて頂けるとありがたい。
前回はThermalZones.Detailed.MixedAirの自然室温計算を見ていった。そのままMixedAirを使っていろいろやってもよかったのだけれど、MixedAirに組み込まれている伝熱現象を理解するためにも先にBuildings.HeatTransferのパッケージについて読んだほうがいい気がしたので少し脱線してこちらに取り組む。全体をざっと見た後に単層壁の非定常熱伝導の計算をする。

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

Buildings.HeatTransferの構成

基本的な伝熱はModelica標準ライブラリだとModelica.Thermal.HeatTransferを使うことになるが、建築で使う形としてモデルを整えているのがこのパッケージである。以下の様な構成になっており、伝導、対流、放射の伝熱の3要素やそれらが複合的に作用する窓のモデル、素材のデータのモデルなどがある。今回はこれらの全体像を見ていくことにする。ちなみに熱伝導のいくつかの部分や窓のモデルにOpenModelicaだと動かない物があるようだ。

名称 説明
UsersGuide ユーザーガイド
Conduction 熱伝導モデルのパッケージ
Convection 対流熱伝達のパッケージ
Radiosity 放射熱伝達のパッケージ
Sources 熱ソース
Windows 窓のモデルのためのパッケージ
Data 熱伝達モデルのためのデータ
Types 型定義のパッケージ
Examples モデルの使用方法とテストモデルを示すモデル群
Interfaces 熱伝達モデルのインターフェースのパッケージ

Buildings.HeatTransfer.Conductionの構成

nSupPCMという定数がConductionフォルダのpackage.moに書きこんである。この定数はPCM(phase change material)という使用温度帯で相変化をすることで見かけ上の熱容量を大きくした素材の熱伝導を計算するのに使用する。
SingleLayerは壁などの面の厚さ方向の一次元熱伝導計算を単一素材で計算するもので、MultiLayerはSingleLayerを複数並べたのもになっている。境界条件は表面のHeatPortで与えている。SingleLayerCylinderは配管の熱損失などに使うと思うのだが単一素材のモデルしかないので複数素材のときは自分で複数のモデルをつないでやる必要がある。
熱伝導のサンプルファイルはHeatTransfer全体のExampleのところにある。

名称 説明
MultiLayer 多層壁の熱伝導のためのモデル
SingleLayer 単層壁の熱伝導のためのモデル
SingleLayerCylinder 円筒の熱伝導
nSupPCM=6 u(T)の近似関係を設定する点数 PCMにのみ使用する
BaseClasses Buildings.HeatTransfer.Conductionのための部分クラス

Buildings.HeatTransfer.Convectionの構成

対流熱伝達は屋外側と室内側に分かれており、流体は空気のみなので液体の対流熱伝達率に使うとすれば固定値で入力することしかできない。室内側は面の向き毎(壁、床、天井)の温度差による関数で計算しており斜面の計算式はなく床でも天井でもなければすべて壁として同一の値になる。屋外側は室内側で使用した温度差の計算式による熱伝達に、風速と風向の関数による強制対流の熱伝達を加えたものになっている。境界条件は固体側と流体側それぞれ表面のHeatPortで与えている。

名称 説明
Exterior 屋外側の対流熱伝達モデル
Interior 室内側の対流熱伝達モデル
Functions 対流熱伝達のための関数
Examples モデルの使用方法とテストモデルを示すモデル群
BaseClasses Buildings.HeatTransfer.Convectionのための部分クラス

サンプルファイルは以下の二つである。

名称 説明
Exterior 屋外側対流熱伝達のテストモデル
Interior 室内側対流熱伝達のテストモデル

Buildings.HeatTransfer.Radiosityの構成

熱放射は放射する物質の温度によって波長が異なってくるが、建築の分野では主に短波(日射)と長波(赤外線、身の回りの物質の放射や夜間放射)に分けて計算する。単純に放射熱伝達と言うと長波のことを指すことが多いが、ここでのモデルも長波のためのモデルを取り扱っている。計算式は基本的に窓の計算のためのモデルTARCOG 2006によるが、OutdoorRadiosityは不透明壁体にも使えるように変更があるらしい。
放射は正の値のみをとるコネクタとして放出Buildings.HeatTransfer.Interfaces.RadiosityOutflowと吸収Buildings.HeatTransfer.Interfaces.RadiosityInflowがそれぞれあり、その差が放射熱伝達量となりHeatPortの熱流として得ることができる。線形近似による計算もできるようなのだが、Heatportが一つしかなく相手側の温度はParameterで与えているため使いどころがちょっとわからない計算になっている。
IndoorRadiosityは放射率1で反射なし、OpaqueSurfaceは反射あり。OutdoorRadiosityは天空温度と外気温度(地表面温度)と天空形態係数から室外側の放射計算をしている。RadiositySplitterは一つの放射のコネクタを割合を指定して2つに分割している。形態係数を使って放射を割り振るにはもう少し機能が欲しいところではないかと思う。

名称 説明
Constant 定常放射熱流束
IndoorRadiosity 室内側熱放射のモデル
OpaqueSurface 不透明表面のモデル
OutdoorRadiosity 窓に到達する屋外側熱放射のモデル
RadiositySplitter 入力信号により入射した放射を二つに分割するモデル
Examples モデルの使用方法とテストモデルを示すモデル群
BaseClasses Buildings.HeatTransfer.Radiosityのための部分クラス

Constantは定数ではなくConstantという名前のBlockであり放射熱流束のコネクタに固定値を与えている。計算上はModelica.Blocks.Sources.Constantとほぼ変わらないのだが、用途が明確になる。

block Constant
  parameter Real k(min=0, start=0)
  extends Modelica.Blocks.Icons.Block;
  Interfaces.RadiosityOutflow JOut;
equation 
  JOut = k;
end Constant;

例は以下の二つである。

名称 説明
OpaqueSurface 不透明表面の室内側のソースのテストモデル
OutdoorRadiosity 屋外側の放射のテストモデル

Buildings.HeatTransfer.Interfacesの構成

Interfacesは放射のものしかなくシンプルなのでここで説明しておく。最低値が0の熱流量のみのコネクタでISO 1の標準基準温度(293.15K、20℃)のときの放射量がデフォルト値となっている。
定義は以下の通り。

名称 説明
RadiosityInflow 放射熱流束の放出側のコネクタ
RadiosityOutflow 放射熱流束の吸収側のコネクタ
connector RadiosityInflow = input Real(min=0, final unit="W", nominal=419)
connector RadiosityOutflow = output Real(min=0, final unit="W", nominal=419)

Buildings.HeatTransfer.Sourcesの構成

Modelica.Thermal.HeatTransfer.Sources.FixedHeatFlowなどをシンプルにしたもの。Prescribedがつくと他のモデルから値をもらってくるのも同じである。FixedHeatFlow0による除算の可能性をなくすためにパラメータのalphaとT_refをなくしているらしいがFixedTemperatureは同じものである。

名称 説明
FixedHeatFlow 固定熱流束境界条件
FixedTemperature 絶対温度[K]での固定温度境界条件
PrescribedHeatFlow 熱流束を入力する境界条件
PrescribedTemperature 絶対温度[K]を入力する境界条件

Buildings.HeatTransfer.Windowsの構成

窓の計算はけっこう複雑である。今回は何があるかを見るだけにして後でもう少し細かく見ていきたい。
BeamDepthInRoomは任意の高さの平面上で日射が窓からどこまで差し込むかを計算する。窓には庇もつけられ、太陽位置は気象データを読み込むようである。ExteriorHeatTransferは外ブラインド等の日射遮蔽物付きの窓の対流と放射の熱伝達のモデルである。屋外側用のモデルなので天空放射や風による対流熱伝達の上昇も見込んでいる。InteriorHeatTransferConvectiveは内ブラインド等の日射遮蔽物付きの窓の対流と放射の熱伝達のモデルで室内側のモデルなので屋外側よりは少しシンプルになっているが少し複雑な印象をうける。FixedShadeは庇と袖壁があったときの窓面日射量と日射面積の割合を出力するモデルで窓のデータのレコードを入力する必要がある。OverhangとSideFinsはそれぞれ庇と袖壁の独立したモデルである。Windowがいわゆる窓の計算のモデルでWindowというソフトのバージョン5時点で同じ計算法(TARCOG 2006 )で実装したらしい。その他BaseClassesにPartialでない普通のmodelやらblockやらがたっぷりある。

名称 説明
BeamDepthInRoom 室内での直達日射の深さ
ExteriorHeatTransfer (日射遮蔽物のある)窓の屋外側表面の対流熱伝達と放射熱伝達のモデル
FixedShade 庇や袖壁による日影のモデル
InteriorHeatTransferConvective (日射遮蔽物のある)窓の室内側表面の対流熱伝達のモデル
Overhang 庇のある窓のモデル
SideFins 袖壁のある窓のモデル
Window 窓のモデル。
Functions 窓の放射モデルに使用する関数
Examples モデルの使用方法とテストモデルを示すモデル群
BaseClasses Buildings.HeatTransfer.Windowsのための部分クラス

窓のサンプルモデルは以下のとおり。BeamDepthInRoomやBoundaryHeatTransferは計算はできており、OverhangやSideFinsはエラーをだしながらも一応計算している。ただしその他の窓の物性値を使ったモデルの計算がうまくいかなさそうである。

名称 説明
BeamDepthInRoom 室内での直達日射の深さのテストモデル
BoundaryHeatTransfer 窓の境界条件の熱伝達のテストモデル
ElectrochromicWindow 電気で透明度が変わる窓
FixedShade 固定日影のテストモデル
Overhang 窓の庇のテストモデル
SideFins 窓の袖壁の使用法
Window 窓のテストモデル

Windowのエラー
f:id:kinonotofu:20181111203906p:plain

Buildings.HeatTransfer.Dataの構成

壁体などに使う素材の物性値をまとめたrecord(変数のセット)用のパッケージである。
ある程度の素材や素材の組み合わせはデフォルトで作られているが、それぞれGenericというrecordに値を入れて好きな素材の物性値を使うことができる。GasesのAirはMediaのAirとは違いGlazingSystemsの窓の中空層で使われている。中空層のAirがよくエラーを吐いているような気がするのだがたぶん直すのはWindowの方だろう。

名称 説明
BaseClasses 部分クラス
BoreholeFillings 掘削穴の充填材
Gases 窓の充填ガスの熱物性
Glasses 窓のガラスの熱物性
GlazingSystems 窓の熱物性
OpaqueConstructions 床や壁などの不透明材の構成
OpaqueSurfaces 不透明材の表面の熱物性
Resistances 熱抵抗
Shades 日射遮蔽物(ブラインドなど)の熱物性
Soil 熱コンダクタンス、密度、比熱容量により記述された土壌
Solids 熱コンダクタンス、密度、比熱容量により記述された固体
SolidsPCM 熱コンダクタンス、密度、比熱容量により記述された固体(PCM)

Buildings.HeatTransfer.Typesの構成

3つの列挙型が定義されている。

名称 説明
ExteriorConvection 屋外側表面の対流熱伝達モデルを定義する列挙型
InteriorConvection 室内側表面の対流熱伝達モデルを定義する列挙型
SurfaceRoughness 表面の粗度を定義する列挙型
  type ExteriorConvection = enumeration(
    Fixed   "Fixed coefficient (a user-specified parameter is used)",
    TemperatureWind   "Wind speed and temperature dependent")

  type InteriorConvection = enumeration(
    Fixed   "Fixed coefficient (a user-specified parameter is used)",
    Temperature   "Temperature dependent")

  type SurfaceRoughness = enumeration(
    VeryRough   "Very rough",
    Rough   "Rough",
    Medium   "Medium rough",
    MediumSmooth   "Medium smooth",
    Smooth   "Smooth",
    VerySmooth   "Very smooth") ;

Buildings.HeatTransfer.Examplesの構成

熱伝導の例ばかりで対流や放射や窓はそれぞれ独立に例を持っているので、Conductionに入れてしまえばいいような気がする。ConductorSingleLayerPCMは変換ができず、ConductorSteadyStateTransientとConductorStepResponseは計算開始直後にエラーを吐いて終了してしまう。

名称 説明
ConductorInitialization 熱伝導の初期化のテストモデル
ConductorMultiLayer 多層壁の熱伝導のテストモデル
ConductorMultiLayer2 表面の層aに条件を与えた多層壁熱伝導のテストモデル
ConductorMultiLayer3 表面の層aとbに条件を与えた多層壁熱伝導のテストモデル
ConductorSingleLayer 単層壁の熱伝導のテストモデル
ConductorSingleLayer2 表面の層に条件を与えない単層壁熱伝導のテストモデル
ConductorSingleLayerCylinder 円筒の熱伝導のテストモデル
ConductorSingleLayerPCM PCMの熱伝導のテストモデル
ConductorSteadyStateTransient 定常条件での過渡応答の熱伝導のテストモデル
ConductorStepResponse ステップ応答の熱伝導のテストモデル

ConductorSingleLayerPCMのエラー
f:id:kinonotofu:20181111201335p:plain

ConductorSteadyStateTransientやConductorStepResponseのエラー
f:id:kinonotofu:20181111201211p:plain

ここでDebug moreを押すと以下のような画面が表示される。使い方がよくわからないがうまく使えばもう少しデバッグがやりやすくなるかもしれない。ただし変換ができないものには使えなさそうである。
f:id:kinonotofu:20181111202914p:plain

PCMの単層壁の計算

なにがあるかをなんとなく見ていっただけになるのも何なので単層壁の熱伝導の計算を確認してみる。PCMの熱伝導の計算はサンプルファイルではエラーがでていたが、PCM自体は普通に計算できるようなので使ってみる。
PCM20というデフォルトで用意されている素材を厚さ0.1m、初期温度30℃として片側を30℃、もう片側を対流熱伝達(壁の温度差で計算)をはさんで0℃として計算した。モデルは以下のようになる。上のほうにあるのがPCMのrecordである。
f:id:kinonotofu:20181111195204p:plain

PCMの設定は以下のとおり。厚さ以外は特にいじっていない。
f:id:kinonotofu:20181111195111p:plain
f:id:kinonotofu:20181111195132p:plain

計算結果は以下の通り。壁体はSinglelayerで素材が一つでも計算は3層に分かれて計算されており、0℃に近い側の層が層変化の温度体(23℃~27℃)で緩やかな温度変化になっている。真ん中の層も27℃になると温度変化がゆるやかになっている。
f:id:kinonotofu:20181111195704p:plain

おわりに

今回はHeatTransferの全体像を見ていき、単層壁の伝熱計算を確認した。
このあたりのモデルは別に建築に限らず使ってもよい気もするけれど建築独特の計算法とか結構あるのだろうか。建築の中でも計算のモデル細かさが違ってきたりするのでなんともいえない部分はあるとは思うがどうだろう。デバッグに関してはMixedAirがOpenModelicaで動かないのは今回の熱伝導のいくつかのサンプルや窓のモデルが動いていないのが原因だと思うのでもう少し細かく見ながらあれこれいじれたらと思う。

今度のModelica勉強会でブログでやったことを少しまとめて発表しようと思います。ブログを書き始めて半年経つそうなのでそろそろ振り返って手を加えたりもしたいしよい機会になればよいかなぁと。ブログでとりあえずライブラリを探検していって勉強会でしっかりまとめるとかでもいいかもしれない。

次回は熱伝導のモデルをもう少し詳しく確認しようと思う。一つの素材の中に状態定義点が複数あるときに表面で定義するかどうかでどういう扱いになっているかとか。あとサンプルでなぜエラーを吐いているかも確認したい。エラーを吐いているのはモデルが多めなので一度減らして少しずつ動作の確認をすればなんとかなるはずだと思いたい。

ModelicaのBuildingsライブラリのThermalZones.Detailed.MixedAirで室温/熱負荷計算を学ぶ_その1

概要

OpenModelicaとJModelicaでLBNLのBuildingsライブラリを使う。初心者なのでいろいろ教えて頂けるとありがたい。
前回までは換気回路網計算を見ていたが、今回は熱回路網計算によって自然室温計算をしているモデルBuildings.ThermalZones.Detailed.Examples.MixedAirFreeResponseを見ていく。本当はBuildings.ThermalZones.Detailed.Validation.BESTEST.Cases6xx.Case600FFの方を見たかったのだけれどどうも計算が途中で止まるのでまずはこちらをみていくことにする。理解が深まったらなにか原因がわかるかもしれない。OpenModelicaでは実行できなかったのでコードを読むのはOpenModelica、実行するのはJModelicaとした。
使用バージョン
-OpenModelica1.12.0 →Modelica標準ライブラリ3.2.2
-JModelica3.4 →Modelica標準ライブラリ3.2.2
-Buildings 5.1.0

Buildings.ThermalZonesについて

ThermalZonesは建物の年間熱負荷計算とか室温変動の非定常計算などをするためのパッケージなのだと思う。室温の計算のためにBuildings.FluidやBuildings.HeatTransferなどのモデルを組み合わせてモデルを作ってくれているようだ。下の表のようにまずは大きく二つに分かれていてDetailedが詳細計算で、ReducedOrderの方はVerein Deutscher Ingenieureというドイツ技術者協会の規格のVDI 6007にある室と建物の過渡的な熱応答の計算方法を実装しているらしい。

名称 説明
Detailed 室モデルのパッケージ
ReducedOrder VDI 6007に基づく簡易化したモデル

どちらも気になるのだがDetailedの方をみることにする。Detailedに含まれるパッケージは下の表の通り。

名称 説明
UsersGuide ユーザーガイド
CFD CFDにより室内空気を計算するモデル
MixedAir 完全混合空気による室のモデル
Constructions 室モデルで使われる壁体構成のためのモデル
FLEXLAB FLEXLAB(多分LBNLの実験施設)をモデル化するために使うモデル
Types 型定義のパッケージ
Examples モデルの使用方法とテストモデルを示すモデル群
Validation 検証モデル群
BaseClasses Buildings.ThermalZones.Detailedのための基底クラスのパッケージ

メインとなるモデルはCFDとMixedAirの二つで、CFDの方は室内の温度などの分布をFast Fluid Dynamics (FFD) のプログラムで計算してModelicaと連成しているらしくBuildings.ThermalZones.Detailed.Examples.FFD.UsersGuideあたりを読みながら使うのだと思う。面白そうだが 計算時間とか精度とかどんなものなのだろう。いずれ使ってみたい。
今回はMixedAirの方を使った検証モデルについて見ていくことにする。このモデルは室内空気が完全混合しているものとして一つの代表点で室温をあらわしている。よくある熱回路網計算である。ちなみにConstructionsはMixedAirの実装に使われているモデルで、FLEXLABはMixedAirを使ってLBNLの実験施設のFLEXLABのセルを再現しているものだと思われる。

MixedAirの方にはコンピュータプログラムによる建物のエネルギー性能評価の試験方法のアメリカの規格(ANSI/ASHRAE Standard 140-2007)に含まれるBESTESTによる方法に従って検証をしたファイルがある。この規格の最新版は2017なので最新の規格に対応できているかまではわからない。ASHRAEはいくつかの規格を公開しているものの140は公開していない。買うと110ドルするのでどこかで見れないものだろうか。いくつかファイルがあるのだが、Case600FFが単室の自然室温の計算をしているもので、これに空調を加えたものがCase600、そしてCase600を継承して多くのファイルが派生している。Cases6xx系統は軽量な(熱容量の小さい)壁体のモデルでCases9xx系統が重い(熱容量の大きい)壁体のモデルになっている。これにも取り組みたいのだが、まずはExamplesにあるもう少しシンプルなモデルMixedAirFreeResponseをみることにする。

MixedAirFreeResponseの構成

最初にモデル全体の構成を確認する。
- 室のモデル(MixedAir)roo
- 気象データweaDat
- 土壌温度TSoi
- 間仕切conOutを介した隣室温度TBou
- 室内の発熱multiplex3_1とそれにつないだ固定値の3種の発熱
- 窓のブラインド等の開閉制御信号replicatorとそれにつないだ固定値
これらを用いて二日間分の単室の室温を計算している。
f:id:kinonotofu:20181030224508p:plain
図において線でつながっていないアイコンがあるがこれはrecord(変数のまとまり)である。

今回パラメータとして定義しているのは室を構成する素材の数と、recordによる素材の構成の定義である。
- 窓の数 parameter Integer nConExtWin = 1;
- 土壌に接する壁の数 parameter Integer nConBou = 1;
- 隣室に接する壁の数 parameter Integer nSurBou = 1;
の3つをパラメータとしており外壁の数もparameter Integer nConExt = 2;としてよいはずだがこれはなぜかrooの方で定義されている。土壌や隣室は勝手にそう呼んでいるだけで、nConBouは室rooの中でConBouで指定した壁面に接する温度、nSurBouは室rooの空気に直接接する温度をつなぐ面の数のことである。

素材の構成はライブラリに登録されたものをそのまま使うことも、素材を指定して厚みを与えて自分で構成を作ることもできる。素材の物性値を変更したい場合は自分で素材のrecordを作ってやる必要がある。
外壁parameter Buildings.HeatTransfer.Data.OpaqueConstructions.Insulation100Concrete200 matLayExt;
間仕切壁parameter Buildings.HeatTransfer.Data.OpaqueConstructions.Brick120 matLayPar;
屋根(上が外側として指定のはず)

  parameter Buildings.HeatTransfer.Data.OpaqueConstructions.Generic matLayRoo(
        material={
          HeatTransfer.Data.Solids.InsulationBoard(x=0.2),
          HeatTransfer.Data.Solids.Concrete(x=0.2)},
        final nLay=2);

床(上が外側として指定のはず)

  parameter Buildings.HeatTransfer.Data.OpaqueConstructions.Generic matLayFlo(
        material={
          HeatTransfer.Data.Solids.Concrete(x=0.2),
          HeatTransfer.Data.Solids.InsulationBoard(x=0.15),
          HeatTransfer.Data.Solids.Concrete(x=0.05)},
        final nLay=3);

窓の素材

  parameter Buildings.HeatTransfer.Data.GlazingSystems.DoubleClearAir13Clear glaSys(
    UFra=2,
    shade=Buildings.HeatTransfer.Data.Shades.Gray(),
    haveInteriorShade=false,
    haveExteriorShade=false);

以上が素材のrecordで、当然これらもアイコンをダブルクリックして変数に値を入力できる。
屋根のmatLayRooの場合以下のような感じになる。両面の長波放射率と日射吸収率はデフォルト値がはいっている。外気側のroughnessは何に使う値かわからない。対流熱伝達率でもかわるのだろうか。
f:id:kinonotofu:20181030225114p:plain

モデルの頭でpackage MediumA = Buildings.Media.Airとしてrooの中で再宣言するのはいつもと同じである。

気象データweaDat

気象データはconnect(weaDat.weaBus, roo.weaBus)としてweaBusをつなぎこんでいる。ファイルを適当に指定して、湿球温度の計算は切っておく。

  Buildings.BoundaryConditions.WeatherData.ReaderTMY3 weaDat(
    filNam=Modelica.Utilities.Files.loadResource("modelica://Buildings/Resources/weatherdata/USA_IL_Chicago-OHare.Intl.AP.725300_TMY3.mos"),
      computeWetBulbTemperature=false)

間仕切conOutと隣室温度TBou

surf_surBouという室rooの空気に直接接続するheatportの配列に多層壁のモデルの配列をconnect(roo.surf_surBou, conOut.port_a)としてつないでいる。多層壁の反対側は15℃固定とした隣室側表面温度TBouにconnect(TBou.port,conOut.port_b)としてつないでいる。多層壁モデルはstateAtSurfaceをtrueにすると表温度を計算するので、今回のように片側の温度を固定している場合はFalseにしてやる必要がある。Aに壁の面積、layersに壁体構成のrecordを入力する。port_aがrecordで上側に定義したものとつながるようになっている。

  Buildings.HeatTransfer.Sources.FixedTemperature TBou[nSurBou](each T=288.15);
  Buildings.HeatTransfer.Conduction.MultiLayer conOut[nSurBou](
    each A=6*4,
    each layers=matLayPar,
    each steadyStateInitial=true,
    each stateAtSurface_a=true,
    each stateAtSurface_b=false);

土壌温度TSoi

室rooの壁体conBouに接続するsurf_conBouというheatportの配列に10℃の固定温度TSoiをconnect(TSoi.port, roo.surf_conBou)として与えている。

  Buildings.HeatTransfer.Sources.FixedTemperature TSoi[nConBou](each T = 283.15)

室内の発熱multiplex3_1

室内の発熱として放射、対流、潜熱流の順番の3つのheatportの配列qGai_flowにmultiplex3_1から発熱をconnect(multiplex3_1.y, roo.qGai_flow)として与えている。放射だけならheaPorRad、対流だけならheaPorAirというheatportがそれぞれ単独で存在するようだ。潜熱流としてMixedAirに内部発熱を定義したいときは要素3つの配列にしてやる必要があるのだと思う。
Modelica.Blocks.Routing.Multiplex3は3つのinputを順にならべた配列にして1つのoutputにするモデルである。inputが配列でも問題ないが今回は単独の数値3つをまとめているのでoutputは要素3つの配列である。値はすべて0なので実質発熱はない。

  Modelica.Blocks.Routing.Multiplex3 multiplex3_1;
  Modelica.Blocks.Sources.Constant qRadGai_flow(k=0);
  Modelica.Blocks.Sources.Constant qConGai_flow(k=0);
  Modelica.Blocks.Sources.Constant qLatGai_flow(k=0);
  connect(qRadGai_flow.y, multiplex3_1.u1[1]);
  connect(qConGai_flow.y, multiplex3_1.u2[1]);
  connect(qLatGai_flow.y, multiplex3_1.u3[1]);

窓のブラインド等の開閉制御信号replicator

uShaという窓のブラインド等の開閉制御信号を入力する配列にReplicatorで作成した配列をconnect(roo.uSha, replicator.y)としてつないでいる。Replicatorはinputの値がnout個並ぶ配列を出力するモデルであり、複数の窓に一括で同じ制御をしたいときなどに便利だと思う。今回は定数で0をconnect(uSha.y, replicator.u) としてつなげて複製しているので日射遮蔽はないものとしていると思われる。

  Modelica.Blocks.Sources.Constant uSha(k=0)
  Modelica.Blocks.Routing.Replicator replicator(nout=max(1,nConExtWin))

これとは別にuWinという窓の状態(電気のオンオフで透過率が変わるなど?)を表すReal型の入力ポート(たぶん日射をカットする割合を入力する)がある。uShaもuWinも0~1の値を入力することになる。窓周りは詳しく調べておきたいところだがとりあえずこんなものがあるのか程度で留めておく。

壁や窓を含む室モデルroo

このモデルが今回の本題である。設定は以下のようになっている。

  Buildings.ThermalZones.Detailed.MixedAir roo(
    redeclare package Medium = MediumA,
    AFlo=6*4,
    hRoo=2.7,
    nConExt=2,
    datConExt(layers={matLayRoo, matLayExt},
           A={6*4, 6*3},
           til={Buildings.Types.Tilt.Ceiling, Buildings.Types.Tilt.Wall},
           azi={Buildings.Types.Azimuth.S, Buildings.Types.Azimuth.W}),
    nConExtWin=nConExtWin,
    datConExtWin(
              layers={matLayExt},
              each A=4*3,
              glaSys={glaSys},
              each hWin=2,
              each wWin=4,
              ove(wR={0},wL={0}, gap={0.1}, dep={1}),
              each fFra=0.1,
              each til=Buildings.Types.Tilt.Wall,
              azi={Buildings.Types.Azimuth.S}),
    nConPar=1,
    datConPar(layers={matLayPar}, each A=10,
           each til=Buildings.Types.Tilt.Wall),
    nConBou=1,
    datConBou(layers={matLayFlo}, each A=6*4,
           each til=Buildings.Types.Tilt.Floor,
           each stateAtSurface_a = false),
    nSurBou=1,
    surBou(each A=6*3,
           each absIR=0.9,
           each absSol=0.9,
           each til=Buildings.Types.Tilt.Wall),
    energyDynamics=Modelica.Fluid.Types.Dynamics.FixedInitial,
    T_start=273.15+22,
    lat=0.73268921998722) "Room model"

床面積、天井高さ、各構成要素の数とその設定などが入力されている。各構成要素の数の入力にパラメーターの変数を使っていないところがあるがミスだろうか。ConExtが外壁で気象データとつながる壁体、ConExtwinが窓、ConParが室内壁(壁の両面がこの室に面している)、ConBouが土壌、SurBouが今回の間仕切壁となっている。それぞれの要素を表にすると以下のような使い分けになっている。この名称の頭にnをつけて要素数、datをつけてそれぞれの設定を記述している。

名称 用途 境界条件を与えるポート
ConExt 外壁 weabus
ConExtwin weabus
ConPar 室内の壁の蓄熱とか放射の効果を見たいとき なし(両面とも室内空気に面する)
ConBou 境界の温度を気象データ以外から与えたいとき surf_conBou
SurBou 壁のモデルをMixedAir以外でつくりたいとき surf_surBou

surBouは壁体がモデルの外でモデル化されるので素材の指定はないのだが、これの指定の変数名はdatSurBouでなくsurBouであることには注意が必要かも知れない。長波放射率(室内の壁同士の放射熱交換用)absIR=0.9と日射吸収率absSol=0.9(窓から入ってきた日射の計算用)を設定している。
Buildingsライブラリで使っている定数として自然換気で風向を指定するときにAzimuthについては確認した。今回は壁体の角度tilに使っている角度の定数Buildings.Types.Tiltを確認する。こちらも単位はラジアンで角度の定数が入っていて、壁面の方向が0[rad]が上向き(水平屋根面など)となっている。

  constant Modelica.SIunits.Angle Ceiling = 0 ;
  constant Modelica.SIunits.Angle Floor =   Modelica.Constants.pi;
  constant Modelica.SIunits.Angle Wall =    Modelica.Constants.pi/2;

窓の設定は袖壁が設定できるようになっている。hWinが窓高さ、wWinが窓の幅、ove(窓の右端から庇の右端の距離,窓の左端から庇の左端の距離, 窓の上面から庇まで距離,庇の深さ)が庇の設定で、このあたりのパラメータは窓のガラス部分のサイズを基準に入力するはずである。右というのが外側から窓に向かって右側であるらしい。またfFraは窓面全体面積にたいするサッシ面積の割合となっている。ちなみに袖壁はsidFin(窓の上面より上に出る長さ,窓の側面からの距離, 袖壁の深さ)袖壁は窓の最下端以下の部分からはじまり、窓の両側に対象的に配置される。

初期温度は室の空気だけでなくここで定義した壁体すべてにも適用される。緯度は気象データのものと一致している必要があるようなので太陽位置などを計算するのだろう。太陽位置は日射量を直達成分と拡散成分に分解したときと壁面の角度に応じて計算するときで一致していなければならない。
今回は使用していないが、MixedAirにはportsとC_flowはfluidportとCO2などの物質の濃度のポートがある。このポートとAirflowのモデルをつかって換気計算も同時に行うことができる。(計算負荷がどの程度になるかは不明)

OMEdit上での設定も確認してみる。
f:id:kinonotofu:20181030230621p:plain
dummyConとdummyGlaSysというダミーのパラメータをインスタンスにあたえるための設定があるが使いどころがよく分からない。
linearizeRadiationは放射計算を温度の4乗で行うと計算負荷が高くなるので放射伝達率などで近似的に線形で計算する設定だと思われる。
hIntFixedは室内側対流熱伝達率でhExtFixedは外気側対流熱伝達率である。ただしデフォルトではこの値は使われておらず、下のintConModとextConModでFixedを選択することでこの値を使用することになると思われる。

f:id:kinonotofu:20181030230838p:plain
steadyStateWindowは窓を定常計算モデルとするかの設定である。窓は熱容量が小さいので定常計算(熱容量なし)とするのはよくあるが、デフォルトでは熱容量を考慮しておりしかもこれがシミュレーションの高速化に寄与するらしい。
mSenFacの室の家具などの熱容量の設定がここにあるが、今回は特別見込まず空気の熱容量のみになっている。

f:id:kinonotofu:20181030231218p:plain
ホモトピー法を使用するかどうかの設定。

f:id:kinonotofu:20181030231339p:plain
ここは実験的に時間のサンプリングの処理により計算の高速化を図ろうとする設定らしいが将来的になくなるかもしれないらしい。

f:id:kinonotofu:20181030231716p:plain
圧力と温度の初期値の設定。

計算の実行

OpenModelicaで実行しようとすると以下のエラーがでる。
f:id:kinonotofu:20181030232145p:plain [3]~[5]はこのファイルを開いたときの警告、[1]と[2]は実行時のエラーである。壁体の配列のあたりでなにかよくないことが起こっているのかもしれない。窓のペアガラスの中空層の空気の部分もファイルを開いたときにエラーとなっている。

とりあえずJModelicaで実行してみる。今回は前回JModelicaを使ったときよりバージョンがあがって2.4を使っている。ファイルが見つからないという警告がでなくなった。少し初期値がみつからないと警告してくるがとりあえず計算はできた。

室温roo.air.vol.Tと外気温度roo.weaBus.DryBul、屋根roo.conEXT[1]と間仕切りconOutの表面温度は以下の通り。計算開始時にがくっと温度が変化している。屋根表面の温度が外気より下がるのは夜間放射の影響と思われる。昼間は日射の影響により屋根や間仕切りの室内側表面温度に比べて室温が上昇している。家具の熱容量は考慮していないがRCの壁を使っているので室温変動はこんなものなのだろうか。
f:id:kinonotofu:20181030235331p:plain

まとめ

とりあえずMixedAirの使い方が分かった。あいかわらず計算結果から知りたい値を探し出すのがたいへんである。伝熱の各要素などの理解を深めたい気もするが、次回はとりあえずBESTESTの計算をどうにかうまくできないかと思う。計算が止まる原因の解決が難しそうなら多数室の計算とか空調時の熱負荷の算出に挑戦してもよいかもしれない。

先日LBNLがModelicaやFMI関係のスタッフを募集していた。

BuildingsライブラリやOpen Building Control, Spawn of EnergyPlus それからIBPSA librariesにも関わる事をするようだ。Modelica関係の開発がさらに活発化するかもしれない。

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にした。
f:id:kinonotofu:20181007185513p:plain

風圧係数の関数について

風圧係数については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にする。それから、テキストビューで媒体の指定を追加してやる。
f:id:kinonotofu:20181007192941p:plain

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上の設定は次のとおり。
f:id:kinonotofu:20181007195525p:plain
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にする。こうすることでスケジュールデータの時刻の単位を秒に変換できる。これで設定ができた。
f:id:kinonotofu:20181007212920p:plain

時系列の換気量の計算

設定したファイルをつなげてやる。weaBusなどはインスタンスの名前がついていないことがあるので注意。
f:id:kinonotofu:20181007213524p:plain
今回追加した部分の接続は以下の通り。

  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)
f:id:kinonotofu:20181007223237p:plain

計算開始後の換気量(上)と風速(下)の関係は以下の通りで風速の大きいところでしっかり換気量が増えているのがわかる。風向毎に風圧係数がかわったり浮力(温度差)の影響があったりするので必ずしもこんなにきれいに相関は出ないのが普通かもしれないけど今回はきれいにでた。
f:id:kinonotofu:20181007225954p:plain
csvで出力すると発熱開始と終了のたびにデータが3つでていた。データの集計の際に面倒なことになるので何かよい方法は無いものかなぁと思った。

おわりに

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

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

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

概要

OpenModelicaでLBNLのBuildingsライブラリを使う。初心者なのでいろいろ教えて頂けるとありがたい。
その1ではBuildings.Airflow.Multizoneパッケージを見てExamples.OneRoomにある温度差換気の計算を確認した。風力換気のモデルが見つからなかったので今回は自分で作って使ってみた(実は風力換気用のモデルがあったのでそれも使ってみた)。

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

作成するモデルの概要

その1で扱ったOneRoomにおいて風が風速2.0[m/s]で吹くものとし、高いほうの窓の風圧係数が-0.2であるとする。低いほうの窓は特に変更しないので風圧係数0.0ということになる。さらに、手計算と比較して検算をしやすくするため、室温と外気温を同じ値にして浮力駆動の換気が起こらないようにする。その他の条件はOneRoomと同じである。

既存モデルの改変

Buildings.Airflow.Multizone.Examples.OneRoomのモデルをベースにするため、モデルを複製する。ファイル→Modelicaクラス新規作成からOneRoomTestというPackageを作り、そこにBuildings.Airflow.Multizone.Examples.OneRoomを複製する。
また、風圧力を考慮したモデルをつくるためにBuildings.Airflow.Multizone.MediumColumnをMediumColumnWithWindという名前で複製しておく。場所はOneRoomTestにする。 このあたりの操作はOpenModelicaでSpaceCoolingを使うときにしたものと同じである。

ファイルを複製したらMediumColumnWithWindに風圧係数と風速を受け取るコネクタを追加し、風圧力を計算する項を追加する。
まず、OMEditのアイコンビューかダイアグラムビューでModelica.Blocks.Interfaces.RealInputをドラッグ&ドロップして名前をvとCpにする。ライブラリブラウザでRealInputで検索をかけてやると見つけやすい。さらにアイコンビュー上でテキストを追加してvとCpとする。既存のテキストを右クリックして複製すると色やサイズなどの設定が楽である。
f:id:kinonotofu:20180908165137p:plain

そして、風圧力を考慮するように以下のように記述する。ここではport_a側(位置が高い方、Top側)に風圧力がかかるようにした。hの最小値は0のままにしているので、このモデルでは高さが低い側に風圧力をかけることはできないことに注意する。

model MediumColumnWithWind "Vertical shaft with no friction and no storage of heat and mass"
(略)
equation 
(略)
  dp = port_a.p - port_b.p;//ここはそのまま
//  dp = -h*rho*Modelica.Constants.g_n
  dp = -h*rho*Modelica.Constants.g_n+Cp*0.5*rho*v^2;//ここは書き換え
(省略)
end MediumColumnWithWind;

このモデルで風圧力を求めるために使っている風圧係数とは風の動圧(速度圧)が建物壁面で静圧に変わる割合のことである。風が吹いているときの壁面の圧力を屋根面平均高さとか軒の高さの動圧で割って算出している。そして風速が大きくなればその動圧に比例して壁面圧力も大きくなるということにしている。風速が変わっても流れ場の性状は風圧係数を求めたときとだいたい同じということを仮定しているのだと思うが、この仮定が有効な風速の条件とかは勉強不足でわからない。レイノルズ数が大きく違ってくるとよくないような気はするがどうなんだろう。

モデルのテスト

風圧力を考慮したモデルができたのでOneRoomで使用するように書き換える。colOutTopのBuildings.Airflow.Multizone.MediumColumnをMediumColumnWithWindに書き換えてやる。ダイヤグラムビューでみると風速と風圧係数をつなぎこむところが使いにくいので右クリックで横反転をしてやる。風圧係数と風速は固定値で与えるのでModelica.Blocks.Sources.Constantをドラッグ&ドロップしてそれぞれCp、とvという名前にする。kにそれぞれ-0.2、2.0を与え、コネクタをつないでやる。
f:id:kinonotofu:20180908183251p:plain

さらに、今は風圧力だけの換気を計算するようにしたいので、浮力が発生しないようにvolOutのInitializationタブのT_startを273.15+20にしてvolEasと同じ値にする。これでモデルを全角文字を含まないパスに保存して、実行する。
f:id:kinonotofu:20180908183316p:plain

手計算だと流量係数0.65、面積0.01m2の開口部を直列合成したときの有効開口面積は0.00459619[m]、密度を1.19684[kg/m3]として圧力差0.478736[Pa]で風量は0.00410554[m3/s]になる。計算結果が0.0041059[m3/s]でだいたい計算はできていそうである。

動圧込みの外気モデル

実はBuildings.Fluid.Sources.Outside_CpBuildings.Fluid.Sources.Outside_CpLowRiseでは気象データで得られる大気圧に風速と風圧係数から算出した動圧を加えた圧力を使っている。Outside_Cpは風圧係数を定数かコネクタから入力でき、Outside_CpLowRiseはwindPressureLowRiseで風圧係数を計算する。ただし、大気圧と風速の入力はweaBusからなのでReaderTMY3からデータを入力する必要がある。そういうことでこれを使って計算してみることにした。
上の開口部はOutside_CpでCpを-0.2にして下の開口部はOutsideを使うようにした。どちらもredeclare package Medium = Mediumをテキストビューで追加し、接続はnPortsを1にしてport[1]にcolOutTopやcolOutBotにつないだ。
気象データの方はReaderTMY3でpAtmを1.01325[bar]、TDryBulを20[degC]、winSpeを2.0[m/s]としてそれぞれParameterから入力とした。filNamは指定しないと計算できなかったので適当に"C:/OpenModelica1.12.0-64bit/lib/omlibrary/Buildings 5.1.0/Resources/weatherdata/USA_IL_Chicago-OHare.Intl.AP.725300_TMY3.mos"を指定した。それから接続するときダイアグラムビューで接続するとconnect(weaBus, weaBus)などとなってしまったのでconnect(weaDat.weaBus, out.weaBus)やconnect(weaDat.weaBus, outcp.weaBus) に修正する必要があった。モデルと計算結果は以下のとおり。
f:id:kinonotofu:20180908200730p:plain
f:id:kinonotofu:20180908201013p:plain
風量は0.00406035[m3/s]でだいたい同じ値である。こちらをはじめから使えばよかったのだけれど気づかなかった。サンプルファイルもBuildings.Fluid.Sources.Examples.Outside_Cpにありだいたい同じようなことをやっている。

おわりに

とりあえず風力換気の計算ができた。次はもう少しだけ実践的に気象データを使って時刻変動する換気量を計算したい。

OpenModelicaでSpaceCooingを動かす

概要

OpenModelicaでLBNLのBuildingsライブラリを使う。チュートリアルSpaceCoolingはStstem2、System3がOpenModelicaではエラーが発生して動かなかったのだが、動くように修正する方法をfinbackさんに教えていただいたので紹介する。finbackさんにはお世話になりっぱなしで本当にありがとうございます。

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

SpaceCoolingと関連モデルの複製とパスの修正

まずは修正したモデル一式を格納するパッケージを作成する。ファイル→Modelicaクラス新規作成を選択しSpaceCoolingTestという名前のPackageを作成する。
f:id:kinonotofu:20180831072244p:plain

SpaceCoolingをパッケージまるごとSpaceCoolingTestパッケージの中に複製する。Buildings.Examples.Tutorial.SpaceCoolingを右クリックして複製を選択。パスをSpaceCoolingTestにしてOKを押す。
f:id:kinonotofu:20180831070526p:plain
f:id:kinonotofu:20180831072430p:plain

さらにBuildings.Fluid.HeatExchangers.WetCoilCounterFlowとBuildings.Fluid.HeatExchangers.DryCoilCounterFlowを同様にして複製する。WetCoilCounterFlowはDryCoilCounterFlowを継承しており、DryCoilCounterFlowがOpenModelicaで問題となっているモデルということになる。複製が終わると以下のようになる。
f:id:kinonotofu:20180831072606p:plain

そしてパスを修正していく。
SpaceCoolingTestの冒頭でimport文によりBuildings.Boundaryconditionsを読み込むようにする。

package SpaceCoolingTest
  import Buildings.BoundryConditions;
  package SpaceCooling "Package with example for how to build a model for space cooling"
(略)

SpaceCoolingのSystem2とSystem3でそれぞれ複製した方のWetCoilCounterFlowを読み込むようにcooCoiの定義部分でBuildings.Fluid.HeatExchangers.WetCoilCounterFlowとなっている部分のBuildings.Fluid.HeatExchangers.を削除してWetCoilCounterFlowに修正する。コメントアウトした部分が元の記述。

(略)
//元はBuildings.Fluid.HeatExchangers.WetCoilCounterFlow cooCoi(略);
  WetCoilCounterFlow cooCoi(redeclare package Medium1 = MediumW, redeclare package Medium2 = MediumA, m1_flow_nominal = mW_flow_nominal, m2_flow_nominal = mA_flow_nominal, dp1_nominal = 6000, dp2_nominal = 200, UA_nominal = -QCoiC_flow_nominal / Buildings.Fluid.HeatExchangers.BaseClasses.lmtd(T_a1 = THeaRecLvg, T_b1 = TASup_nominal, T_a2 = TWSup_nominal, T_b2 = TWRet_nominal), show_T = true, energyDynamics = Modelica.Fluid.Types.Dynamics.FixedInitial) "Cooling coil" annotation(
    Placement(transformation(extent = {{-10, -10}, {10, 10}}, rotation = 180, origin = {-30, -26})));
(略)

DryCoilCounterFlowの冒頭にimport文を追加してBuildings.Fluid.HeatExchangers.BaseClassesを読み込む。

model DryCoilCounterFlow "Counterflow coil with discretization along the flow paths and without humidity condensation"
  extends Buildings.Fluid.Interfaces.PartialFourPortInterface(show_T = false);
  extends Buildings.Fluid.Interfaces.FourPortFlowResistanceParameters(final computeFlowResistance1 = false, final computeFlowResistance2 = false, from_dp1 = false, from_dp2 = false);
  import Buildings.Fluid.HeatExchangers.BaseClasses;
(略)

モデルの修正

DryCoilCounterFlowにおいて、replaceable BaseClasses.HexElementSensible ele[nEle] constrainedby BaseClasses.PartialHexElement(略)となっている部分を修正する。replaceable文ではクラスを交換することはあってもインスタンスを交換することはあまりないとのこと。また、constrainedbyでクラスを限定するような部分でカッコ内の引数をredeclareするのも一般的ではないのではないかとのこと。MElementというreplaceableなモデルをつくり、これでele[nEle]インスタンスを作成する。コメントアウトした部分が元の記述。

(略)
//  replaceable BaseClasses.HexElementSensible ele[nEle] constrainedby BaseClasses.PartialHexElement();
replaceable model MElement = BaseClasses.HexElementSensible constrainedby BaseClasses.PartialHexElement;
MElement ele[nEle](redeclare each package Medium1 = Medium1, redeclare each package Medium2 = Medium2, each allowFlowReversal1 = allowFlowReversal1, each allowFlowReversal2 = allowFlowReversal2, each tau1 = tau1 / nEle, each m1_flow_nominal = m1_flow_nominal, each tau2 = tau2, each m2_flow_nominal = m2_flow_nominal, each tau_m = tau_m / nEle, each UA_nominal = UA_nominal / nEle, each energyDynamics = energyDynamics, initialize_p1 = {i == 1 and not Medium1.singleState for i in 1:nEle}, initialize_p2 = {i == 1 and not Medium2.singleState for i in 1:nEle}, each deltaM1 = deltaM1, each deltaM2 = deltaM2, each from_dp1 = from_dp1, each from_dp2 = from_dp2, dp1_nominal = {if i == 1 then dp1_nominal else 0 for i in 1:nEle}, dp2_nominal = {if i == nEle then dp2_nominal else 0 for i in 1:nEle}) "Heat exchanger element" annotation(
     Placement(transformation(extent = {{0, 0}, {20, 20}})));
(略)

WetCoilCounterFlowでele[nEle]はredeclareされているので対応するように修正。こちらもインスタンスを再宣言するのではなくクラス(model)を再宣言するようにする。コメントアウトした部分が元の記述。

(略)
//extends Buildings.Fluid.HeatExchangers.DryCoilCounterFlow(redeclare replaceable package Medium2 = Modelica.Media.Interfaces.PartialCondensingGases, redeclare final Buildings.Fluid.HeatExchangers.BaseClasses.HexElementLatent ele[nEle]);
extends Buildings.Fluid.HeatExchangers.DryCoilCounterFlow(redeclare replaceable package Medium2 = Modelica.Media.Interfaces.PartialCondensingGases, redeclare final model MElement=Buildings.Fluid.HeatExchangers.BaseClasses.HexElementLatent);
(略)

計算実行

修正が終わったファイル一式を保存して実行する。たぶんファイルパスに日本語とか全角文字が含まれていると対応できない。
計算ができなくておかしいと思っていたが、WetCoilCounterFlowを右クリックしての順番をDryCoilCounterFlowの下に変えてやると計算できた。当然といえば当然だが継承元を先に記述しないといけないのだろう。
f:id:kinonotofu:20180831083102p:plain

System2の室温は以下のとおり。
f:id:kinonotofu:20180831081118p:plain

System3の室温は以下のとおり。
f:id:kinonotofu:20180831081530p:plain

おわりに

無事にOpenMoelicaで動作させることができた。多分同じような修正(インスタンスをreplaceable、redeclareしているような部分を書き直す)で他にも現在OpenModelicaで動いていない部分が動くようになったりするのではないかと思う。この修正でDymolaで困ることがなければプルリクやイシューで頼めば直してもらえるはず。
修正方法を教えて頂いたfinbackさんには本当に感謝します。

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

概要

OpenModelicaでLBNLのBuildingsライブラリを使う。初心者なのでいろいろ教えて頂けるとありがたい。
今回はBuildings.Airflow.Multizoneパッケージを見ていく。Buildings.Airflowパッケージでは異なる部屋間や部屋と外部環境との間の空気の流れを計算するためのモデルがあるのだが、直下にはMultizoneパッケージのみがぶら下がっていてそこに各要素のモデルがある。ダクトのネットワークはBuildings.Fluidを見ろと書いてあるのだが多分FixedResistancesとかMoversあたりを使うのだと思う。
このパッケージでできる換気回路網計算はCFDとは異なり、空間をそれぞれ1つの点で表現し、空間同士を開口による流れの抵抗のネットワークとしてつないでいくものだ。年間の自然換気量の変化の影響を見たいときなどCFDで非定常計算をするのはあまり実用的ではないのでこちらを使ったりする。定常計算でもこちらの方がざっくり設定できるし1ケースの計算時間もたぶん早いので都合がよいことも多い。

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

Buildings.Airflow.Multizoneの構成

名称 説明
UsersGuide ユーザーガイド
DoorDiscretizedOpen 高さに沿って離散化したドアモデル
DoorDiscretizedOperable 高さに沿って離散化したドアモデル
EffectiveAirLeakageArea 有効漏気面のモデル
MediumColumn 摩擦がなく熱と質量の蓄積がない垂直軸のモデル
MediumColumnDynamic 摩擦はないが熱と質量の蓄積がある垂直軸のモデル
Orifice オリフィス(単純開口)のモデル
ZonalFlow_ACS 1秒あたりの換気量の入力によるゾーン流のモデル
ZonalFlow_m_flow 1秒あたりの換気量の入力によるゾーン流のモデル
Types 型の定義のパッケージ
Examples モデルの使用法やテストモデルなどのパッケージ
Validation バリデーション用のモデルのパッケージ
BaseClasses Buildings.Airflow.Multizoneの基底クラスのパッケージ

ユーザーガイドと型定義やサンプルや検証や基底クラスなどのパッケージがあるが実際に使用するモデルは8個である。正直ユーザーガイドには大したことは書いていないのでExamplesを見ながらなんとなく組み上げていくのだと思う。
Typesパッケージで定義されているものは1つだけでdensitySelectionというMediumColumnモデルの密度の設定パラメータとなる列挙型の変数があるのみである。

type densitySelection = enumeration(
  fromTop   "Density from top port",
  fromBottom   "Density from bottom port",
  actual   "Actual density based on flow direction")
  "Enumeration to select density in medium column";

Examplesパッケージは結構たくさん例がある。

名称 説明
CO2TransportStep 浮力駆動によるCO2輸送モデル
ChimneyShaftNoVolume シャフトの定常状態モデルでの煙突効果を示すモデル
himneyShaftWithVolume シャフトの動的モデルによる煙突効果を示すモデル
ClosedDoors 3つの閉じたドアのモデル
NaturalVentilation 密度差による流れの逆転のテストモデル
OneEffectiveAirLeakageArea 有効漏気面を持つモデル
OneOpenDoor 1つの開いたドアと1つの閉じたドアのモデル
OneRoom マルチゾーン空気交換モデルの検証のための1つの室のモデル
Orifice オリフィスのモデル
ReverseBuoyancy 4つの室と逆流する浮力駆動型空気循環モデル
ReverseBuoyancy3Zones 3つの室と逆流する浮力駆動の空気循環モデル
ZonalFlow 2つの室で所定の空気交換を行うモデル

ValidationパッケージはThreeRoomsContamモデル一つだけである。Michael WetterのMultizone Airflow Model in Modelicaに詳細な内容が書いてあるようだ。

Buildings.Airflow.Multizone.Examples.OneRoom

とりあえずサンプルとして一室モデルを見ていく。
f:id:kinonotofu:20180807233745p:plain
室volEasと外気volOutを、上下の開口部oriOutTopとoriOutBotでつないで温度差換気を行っているモデルである。室と開口部の間はcolEasInTop、colEasInBot、colOutTop、colOutBotという高さ方向の距離を示すモデルでつないでいる。

  connect(colEasInTop.port_a, oriOutTop.port_a);
  connect(colEasInTop.port_b, volEas.ports[1]);
  connect(colEasInBot.port_a, volEas.ports[2]);
  connect(colEasInBot.port_b, oriOutBot.port_a);
  connect(oriOutBot.port_b, colOutBot.port_b);
  connect(colOutBot.port_a, volOut.ports[1]);
  connect(colOutTop.port_b, volOut.ports[2]);
  connect(colOutTop.port_a, oriOutTop.port_b);

空気の流れの計算なので媒体としてBuildings.Media.Airを使用している。

  package Medium = Buildings.Media.Air;

空間のモデルの設定

室volEasの設定はSpaceCoolingで見たときと比べて、熱容量の係数を使っていない(家具などの熱容量を考慮しない)、massDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitialを指定する、m_flow_nominal=0.001が0に近い小さな値になっているなどの違いがある。
以下に設定を示す。

  Buildings.Fluid.MixingVolumes.MixingVolume volEas(
    redeclare package Medium = Medium,
    energyDynamics=Modelica.Fluid.Types.Dynamics.FixedInitial,
    T_start=273.15+20,
    V=2.5*5*5,
    nPorts=2,
    m_flow_nominal=0.001,
    massDynamics=Modelica.Fluid.Types.Dynamics.SteadyStateInitial);

外気volOutはモデルとしては室と同じものを使用している。
体積を非常に大きくとっているのは室との換気で温度が変動しないようにするためだろうか。また、外気側では圧力を指定している。外気を圧力の基準とするような意味合いだと思うがMedium.p_defaultで大気圧を与えている。こちらのmassDynamicsはFixedInitialになっている。これで圧力固定になるのだろうか。

  Buildings.Fluid.MixingVolumes.MixingVolume volOut(
    redeclare package Medium = Medium,
    energyDynamics=Modelica.Fluid.Types.Dynamics.FixedInitial,
    T_start=273.15 + 10,
    V=1E12,
    p_start=Medium.p_default,
    nPorts=2,
    m_flow_nominal=0.001,
    massDynamics=Modelica.Fluid.Types.Dynamics.FixedInitial);

開口部のモデルの設定

上下の開口部oriOutTopとoriOutBotをまとめて示す。

  Buildings.Airflow.Multizone.Orifice oriOutTop(
    redeclare package Medium = Medium,
    A=0.01,
    m=0.5);
  Buildings.Airflow.Multizone.Orifice oriOutBot(
    redeclare package Medium = Medium,
    A=0.01,
    m=0.5);

OMEdit上の設定画面は以下のとおり。
f:id:kinonotofu:20180807234101p:plain
オリフィスの流量係数は0.65がデフォルトになっており今回はデフォルト値をつかっている。ForceErrorControlOnFlowAはよくわからないが流量を見たいときはtrueにしておくらしい。Aが開口面積、mは0.5が乱流時で1が層流時を示すらしい。useDefaultPropertiesはfalseにすると実際の密度と粘性係数を使用し、trueだとデフォルト値を使用するようになっている。dp_turbulentは層流と乱流が一致するときの差圧?だが0.1を推奨しているらしい。lWetはレイノルズ数の計算に使用するパラメーターなので水力直径だろう。密度の扱いは計算の安定性に関わるので設定があるのはわかるけれど、粘性係数は気にしたことがなかった。乱流層流の切替も気にしたことがなかったけど確かに流量が変わりそうな気はする。 
以下の残りの2つのタブは流れに関係する設定でよくあるやつだった。
f:id:kinonotofu:20180807234116p:plain
f:id:kinonotofu:20180807234126p:plain

高さ方向の距離を示すモデルの設定

Buildings.Airflow.Multizone.MediumColumnは開口部の高さ位置を指定するためのモデルなのだと思う。浮力換気のための高低差による圧力の変動を計算しているのだと思う。設定を以下に示す。

  Buildings.Airflow.Multizone.MediumColumn colOutTop(
    redeclare package Medium = Medium,
    h=1.5,
    densitySelection=Buildings.Airflow.Multizone.Types.densitySelection.fromBottom);
  Buildings.Airflow.Multizone.MediumColumn colOutBot(
    redeclare package Medium = Medium,
    h=1.5,
    densitySelection=Buildings.Airflow.Multizone.Types.densitySelection.fromTop);
  Buildings.Airflow.Multizone.MediumColumn colEasInTop(
    redeclare package Medium = Medium,
    h=1.5,
    densitySelection=Buildings.Airflow.Multizone.Types.densitySelection.fromBottom);
  Buildings.Airflow.Multizone.MediumColumn colEasInBot(
    redeclare package Medium = Medium,
    h=1.5,
    densitySelection=Buildings.Airflow.Multizone.Types.densitySelection.fromTop);

OMEdit上の設定画面は以下のとおり。
f:id:kinonotofu:20180808000457p:plain
f:id:kinonotofu:20180808000649p:plain
hが高さ方向の距離を示し、Topにつないだほうが上、Bottomにつないだほうが下という位置関係になるようだ。ここでの室の圧力の基準位置がどこだかわからないが、基準位置を床にして上に0.5と3.5などのように床面と開口面高さ方向中央位置で指定したほうがよさそうな気がする。このあたりは自分で整合性をとっておけばたぶん問題はないとは思う。
densitySelectionではモデルのどちら側の密度を使用するかを指定している。TopとBottomの間にはρghだけ圧力差があるのだが、このときのρをどちら側からとってくるかを選んでいるのだろう。actualという選択肢もあるのでそのときは風上側の密度を使用するのだと思う。

計算実行

温度差換気(浮力駆動の換気)で気にすることはだいたい以下のとおり。
* 外気はどの高さでも同じ温度、室内もどの高さでも同じ温度(完全混合の仮定)
* 温度(=室)によって密度ρが変わる。
* 高さにより「密度ρ × 重力加速度g × 高さh」だけ圧力が変わる。
* 温度(=密度)が異なる二つの空間の間の圧力差は高さによって変わる。
* 圧力差によって流れが発生。「開口面の高さでの圧力差により流れを計算」するので同じ室の間でも開口部ごとに流れの方向がかわることができる。

計算結果は以下の通り。とりあえず1秒だけ計算している。
f:id:kinonotofu:20180808002843p:plain
自分で計算をしてみると0.0065[m3/s]だったのでだいたいよさそうである。
1時間分計算してみると室温が下がっている様子がわかる。
f:id:kinonotofu:20180808003437p:plain
圧力は以下のとおり。外気の圧力が一定なのが確認できた。室の圧力が外気を上回ったり下回ったり変動しているように見えるがレンジの設定が細か過ぎるだけでほぼ外気と同じ圧力となっている。ちなみにcsvでデータ出力すると1.01325[bar]で一定値となってでてくる。桁落ちされるのは困るのでどうにかしたいが大気圧との差分をどこかで計算すれば場当たり的には対処できそうな気はする。しかしなんでグラフ上で見える値をそのまま書き出せないのだろう。ちなみに単位が[bar]なのは[Pa]に直すのを忘れていただけでグラフ上もcsvも[Pa]にできる。

f:id:kinonotofu:20180808011533p:plain

おわりに

なんとなく雰囲気は分かったような気がするけれど、まだあまり仕組みはよく分かっていない。Fluidportをつないであげれば開口部では圧力から流量と計算し、室のほうで流量の収支をとるように圧力を補正してくれているのだろうか。開口への風圧の与え方もわからないのだけれど、別のFluidポートをつなぎこんで圧力だけあたえるとかはできなさそうな気がするのでMediumColumnに手を加えて風圧力を組み込むしかないのかもしれない。あと一応ちゃんと論文を読んだほうがよさそうである。