第一回電脳建築最適化世界選手権に参加しませんか

f:id:kinonotofu:20190305102942p:plain
このコンテストの関係者ではないが宣伝してみる。
電脳建築最適化世界選手権(WCCBO:World Championship in Cybanetic Building Optimization)という空気調和衛生工学会主催のコンテストがある。サイトをみるだけでは内容が分かりにくいが一応フライヤーに概要が書いてはある。
参加が決まっている人はマニュアルを読み込めば問題ないので、面白そうなら参加したいという人に向けて全体像を軽く把握できるようにしたいと思う。なお、解釈が間違っている可能性はあるし、今後内容が変更されたりするかもしれない。

概要

このコンテストは誤解を恐れずにざっくり言うと室や熱源の「設定温度の年間のスケジュール」を最適化するコンテストである。

電脳建築最適化というタイトルだが平面プランや構造などの最適化は全く関係がない。「空調制御を最適化して(省エネ性×快適性)の値の大きさを競う」というコンテストである。温度が変われば快適性や機器効率が変わるので上手にバランスを取る必要がある。実物件では条件が揃えられないため、エミュレータを使用して人間の位置や室内環境や熱源システムを模擬することで参加者全員同一条件としている。エミュレータでのバーチャルな世界を想定してるためか「電脳」とついている。デジタルツインのようなものに応用したいのかもしれない。

内容としてはエミュレータ上で設定温度などの年間のスケジュールを決めた「空調制御についての設定ファイル」を提出して「省エネ性×快適性」によるスコアを競うものになる。さらに任意参加ではあるが、リアルタイム(時間はエミュレータ上で圧縮されている)に設定温度を送信して空調制御を行う部門も存在し、自由度が飛躍的に高まることになる。

エミュレータは配布されるので自分で好きなだけ試すことができ、エミュレータ内の時間の進行速度はある程度(1~960倍)調整できる。ただし、管理者権限が必要なので会社で参加したいが権限が貰えない様な場合が仮にあるならば面倒かもしれない。これはBACnet(建築設備の自動制御で使うプロトコル)通信を行うためなのだが、BACnetについての事前知識は全く必要ない。また、任意参加の部門では運営のサーバーにVPN接続を行う必要がある。エミュレータWindowsの実行ファイルだが、C#ソースコードも配布される。MacLinuxで自力でビルドするのが現実的かについては現在不明である。

スケジュール

(~4/12金)参加申し込み

参加申し込みが必要で4/12までに申し込む必要がある。参加費なし賞金なし、誰でも匿名でも参加可能でチームでも個人でもよい。資料が郵送されるようなので代表者の住所は必要。
参加申し込み後に資料が送付されるように書いてあるが、すでに公式サイトに資料がアップロードされている。ただし、開始するまでに内容が更新されるかもしれない。

(6/7火)説明会、オフライン部門開始

空衛学会の会議室で説明会のようなものがあるが、参加は任意。遠隔地用に動画の配信があるかもしれない。この日にとりあえず一度提出をして翌日には順位が公開される。

(6/7金~7/7日)オフライン部門

エミュレータ上で年間の空調制御を決めて作成した「空調制御についての設定ファイル」を提出し、主催のサーバー上で計算した結果を競う。6/7に一度提出すればその後何もしなくても問題はない。サーバーでの計算は1日かかるが、期間内であれば何度でも再提出してスコアを更新できる。

(7/8月~8/7水)オンライン部門、任意参加

この期間は運営サーバーでエミュレータの速度を12倍にして1ヶ月で1年間を模擬している。参加者はリアルタイムに状況を把握しながらVPN接続で制御の変更を行う。基本的にプログラムを用意して一ヶ月間休みなくPCを動かすか、クラウドサービスを利用することになると思う。PCに張り付いて手動で操作してもよい(そのための実施期間?)。やり直しが効かない。

対象

マニュアルの付属資料に対象建物の平面図や機器表と系統図やゾーニング、制御のブロック図と機器効率の曲線などがある。対象はオンライン、オフライン共通と思われる。公平性から考えると参加申し込み後に送付される内容が現在と違ってもおかしくはないと思うが、そこまで配慮しないかもしれない。

  • 7階建で延べ床10,000m2くらいのテナントビル
  • 人員1000人近くが確率的に動き回る(乱数のシード値は参加者共通)
  • 中央熱源式で主な機器はガス焚き吸収式冷温水発生機と冷却塔と空冷モジュールHPで温度成層型の蓄熱槽がある
  • 電気室や中央監視室などの空冷PACは対象外

入力

  • 設定温度等の制御パラメータを指定する
  • 一つ一つ設定するのは大変だがエクセルで一括設定ができる(その他のインターフェースもある)
  • オフライン部門の制御パラメータは最大20の日程のグループに分けてそれぞれ設定する。設定ファイルを提出する
  • オンライン部門ではリアルタイムに設定値を送ることで独自制御アルゴリズムを実現可能。VPN接続で直接設定する

出力

  • 月1回の快適性アンケート
  • GUIで表示されるデータの1分間隔のログのcsvファイル
  • その他のGUIで表示されない情報は基本的に得られない
  • デフォルトの制御に対する省エネ性と快適性の値(目的変数)は毎分得られる

GUI

以下のようなBASを模擬したものが用意されている。
f:id:kinonotofu:20190305110315p:plain
f:id:kinonotofu:20190305110339p:plain

おわりに

電脳建築最適化世界選手権は以下のようなものである。

  • 設定温度等を入力して(省エネ性×快適性)の大きさを競う
  • 省エネ性や快適性を算出するエミュレータは与えられる
  • 誰でも無料で参加可能
  • 参加申し込み(4/12まで)
  • オフライン部門(6/7~7/7)、膨大な制御パラメータの最適化
  • オンライン部門(7/8~8/7)、最適な制御モデルの構築

空調に詳しくない人にも興味を持ってもらえればと思うのだけれど、年間計算最短9時間10分くらいかかるようなので力技の最適化は少しつらいかもしれない。第一回なので運営側も手探りなところが多いと思うけれど盛り上がればいいなと思う。

ModelicaのBuildingsライブラリのHeatTransferを学ぶ_その3 PCMの計算

概要

OpenModelicaでLBNLのBuildingsライブラリを使って熱伝導の計算を学ぶ。間違い等あればいろいろ教えて頂けるとありがたい。
前回は単層壁の1次元熱伝導モデルBuildings.HeatTransfer.Conduction.SingleLayerのPCM以外の部分を見たので今回はPCMの計算部分を確認する。

使用バージョン

  • Buildings 5.1.0

PCMの計算方法の概略

基本的に物質は相変化のために熱を使用している間は温度が一定となる。PCM(Phase Change Material)はその性質を利用して利用したい温度帯での見かけの熱容量を大きくして温度を安定させる材料となる。ただし、相変化の温度帯の調整や建材として使用するなどの理由により複数の材料を混ぜることで、相変化の間でもPCMの温度は変化していることが多い。BuildingsライブラリのPCMの計算で参考している文献にあるPCMを30%混ぜた板材の熱量と温度の関係のモデルはだいたい以下のようになる。
f:id:kinonotofu:20190130014315j:plain

この図は計算モデルを表しているものであり、比熱は固相、液相、固液混相域で同じ値を用いたものであるが、縦軸の比エンタルピー(熱量)と横軸の温度の関係を見ると25℃付近での温度変化に対する必要熱量が大きくなっていることがわかると思う。本来は液相の方がPCMの比熱は少し大きくなるらしい。Buildingsライブラリの実装では比熱や密度や熱伝導率は固相と液相で同じ値を使用している。さらに熱量に関しては比エンタルピーではなく比内部エネルギーを定義して比内部エネルギーと温度の関係としている。内部エネルギーの変化量=エンタルピーの変化量+仕事量なので圧力や体積が変わらない(無視できる範囲内)ならば仕事をせず同じように扱えるということなのだと思う。ちなみに「比」というのは「単位質量あたりの」という意味である。

Cubic Hermite splineについて

Buildingsライブラリでは比内部エネルギーと温度の関係を表すために3次エルミートスプラインを使用している。基本的には固相や液相の時と固液混相のときでそれぞれ傾きの異なる直線としてモデル化したかったのだと思うのだが、おそらく傾きが不連続になると数値計算上収束性が悪くなるなどの問題がでるためにこのような方法をとっているものと考えられる。ちなみに参考文献(上のPCMの資料)では双曲線関数を使用していたようなのだが、どちらがよいのかはよくわからない。3次エルミートスプラインのほうが適当に近似曲線をつくるのは楽そうではある。
3次エルミートスプラインとは以下のようなものである。

  • エルミート補間は各点の値と微分を一致させる補間。だいたい3次のスプライン。
  • スプラインは区分多項式のこと。区分(点と点の間)毎に数式を定義する。
  • 3次式とするのは各区分の両端の値と微分の合計4つの値を条件とすると3次式が一意に定まるから。

数式の内容は英語のWikipediaCubic Hermite splineに書いてあるとおりである。Wikipedeiaの式を今回の実装に対応させて書き直すと以下のようになる。
点(x1,y1) の微分dy/dx1 と点(x2,y2) の微分dy/dx2 が与えられたとき x1 < x < x2 である x に対応する y を求める。仮にx1=0、x2=1であったときにはyの値は以下のようになる。
f:id:kinonotofu:20190131124123p:plain
f:id:kinonotofu:20190131124137p:plain
f:id:kinonotofu:20190131124153p:plain
f:id:kinonotofu:20190131124206p:plain
f:id:kinonotofu:20190131124217p:plain
f:id:kinonotofu:20190131124239p:plain

そしてこれにアフィン変換(今回は回転やせん断はなく平行移動と拡大縮小のみ)を施して一般的な形にすると以下のようになる。h00~h11の係数は上の式と同じである。
f:id:kinonotofu:20190131124303p:plain
f:id:kinonotofu:20190131124315p:plain
f:id:kinonotofu:20190131124327p:plain

この計算モデルがModelica標準ライブラリで用意されている。微分dy/dxについては自分で用意する必要があるのだが、今回の実装では⊿を各点を結んだ直線の傾きとすると、各点につながる直線の傾きの平均を各点の微分としている。
f:id:kinonotofu:20190131124805p:plain
f:id:kinonotofu:20190131124847p:plain

両端については片側にしか直線がないのでその値を使用する。
f:id:kinonotofu:20190131124923p:plain
f:id:kinonotofu:20190131124953p:plain

さらに補間に使用する微分の値に補正をかけることで単調性を確保することができる。つまり与えられた点が単調増加や単調減少なら補間した値も単調性を保つようにできるのだ。こちらはBuildingsライブラリで実装されており、以下の論文を参考にしているらしい。
F.N. Fritsch and R.E. Carlson, Monotone piecewise cubic interpolation. SIAM J. Numer. Anal., 17 (1980), pp. 238-246
論文は有料だが同著者の似たような記事でやや古いPDFがあるようだ。また、英語のWikipediaMonotone cubic interpolationもこの内容である。Wikipediaの図を見れば補正がなぜ必要になるのかがわかると思う。
各点の微分を修正していくのだが、各点をつなぐ直線を左から調べて、両端の点の微分を順次修正していくため、右側の点は次の直線を調べたときに再度微分の値が修正されることになる。
まず直線の傾きが0になる場合はその両端の点の微分は0であるとする。
f:id:kinonotofu:20190131125140p:plain
f:id:kinonotofu:20190131125210p:plain

それ以外の場合はまず、単調性を確保させる際に発散させないように以下の条件を調べる。
f:id:kinonotofu:20190131125344p:plain
f:id:kinonotofu:20190131125433p:plain
f:id:kinonotofu:20190131125444p:plain

今回の実装ではこれが満たされない区間は補正をしないだけになり警告がでないため、必ず単調性が確保されているわけではないかもしれない。あらかじめ出てきた曲線を調べておいたほうがよさそうである。条件を満たす場合には以下のような処理を行う。
f:id:kinonotofu:20190131125511p:plain
f:id:kinonotofu:20190131125529p:plain
f:id:kinonotofu:20190131125543p:plain

このような処理を繰り返していけば単調性が確保されるらしい。証明までは踏み込まなかったけれど、とりあえずこのようにして計算できるようだ。

初期条件の定義とPCMの計算に用いるパラメータの設定

まず、複数のクラスで共有する定数としてnSupPCMという定数がConductionフォルダのpackage.moに書きこんである。この数は6に設定してあるのだが、PCMの層を物理的に6分割して節点を求めているのではなく、補間に使用する点が6箇所であるという意味である。補間に使用する点数を増やしたい場合には他の部分のコードも修正する必要がでてくるので注意が必要である。

SingleLayerにおいてPCMの計算のみに使用する変数として、3次エルミートスプライン用の比内部エネルギーud(横軸)、温度Td(縦軸)、微分dT_duが要素数6つの配列として定義されており、これらの値を使用して温度を計算することになる。

  parameter Modelica.SIunits.SpecificInternalEnergy ud[Buildings.HeatTransfer.Conduction.nSupPCM](
    each fixed=false);
  parameter Modelica.SIunits.Temperature Td[Buildings.HeatTransfer.Conduction.nSupPCM](
    each fixed=false);
  parameter Real dT_du[Buildings.HeatTransfer.Conduction.nSupPCM](
    each fixed=false,
    each unit="kg.K2/J");

近似式は共通のものを使っているものの、温度や比内部エネルギーや熱抵抗や熱容量の定義点の位置は前回の通りであり、場所によってPCMの温度などの状態(固相か液相か固液混相か)は変わり得る。PCMの状態が全ての温度節点で同じ状態である必要はない。

初期条件の定義部分は以下の通り。前回も少し触れたのだが、初期化のときにsteadyStateInitial=trueの場合は微分0の条件が入るが、そうでない場合温度に初期値は与えていても内部エネルギuには値を与えていない。微分0にしても比内部エネルギーを定義するコネクタがないので実質比内部エネルギーの初期条件が定義できていないような気がする。equationの方に温度と比内部エネルギーの関係式があるのでそれを使用して初期温度から比内部エネルギーの初期値を定義してくれているのかもしれない。このあたりを自動でやってくれるのがModelicaのよいところなのだろう。
そして、先ほどの3つの配列に対してmaterialの値とBuildings.HeatTransfer.Conduction.BaseClasses.der_temperature_uという関数を用いて値を定義している。parameterの変数はequationでは値を変更できないが、initial equationの部分であれば数式で値を定義することができるらしい。ensureMonotonicityは3次エルミートスプラインの説明で触れた微分に補正をかけるかどうかの選択でtrueにすると単調性を確保してくれる。

    if not material.steadyState then
      if steadyStateInitial then
        if material.phasechange then
          der(u) = zeros(nSta);
        else
          der(T) = zeros(nSta);
        end if;
      else
        if stateAtSurface_a then
          T[1] = T_a_start;
          for i in 2:nSta loop
            T[i] =T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i-1);
          end for;
        else // stateAtSurface_a == false
          for i in 1:nSta loop
            T[i] = T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i);
          end for;
        end if;
    end if;
   if material.phasechange then
     (ud, Td, dT_du) = Buildings.HeatTransfer.Conduction.BaseClasses.der_temperature_u(
       c =  material.c,
       TSol=material.TSol,
       TLiq=material.TLiq,
       LHea=material.LHea,
       ensureMonotonicity=material.ensureMonotonicity);
   else
     ud    = zeros(Buildings.HeatTransfer.Conduction.nSupPCM);
     Td    = zeros(Buildings.HeatTransfer.Conduction.nSupPCM);
     dT_du = zeros(Buildings.HeatTransfer.Conduction.nSupPCM);
   end if;

関数der_temperature_uの中身

Buildings.HeatTransfer.Conduction.BaseClasses.der_temperature_uという関数によってud、Td、dT_duが定義されているので中身を確認する。入力と出力は以下の通り。

入力

名称 デフォルト値 説明
SpecificHeatCapacity c - 比熱[J/(kg.K)]
Temperature TSol - 固相線[K]
Temperature TLiq - 液相線[K]
SpecificInternalEnergy LHea - 相変化時の潜熱[J/kg]
Boolean ensureMonotonicity false Trueの時に強制的に3次エルミートスプラインを単調にする

出力

名称 説明
SpecificInternalEnergy ud[nSupPCM] 定義点でのud[J/kg]
Temperature Td[nSupPCM] 定義点でのTd[K]
Real dT_du[nSupPCM] 定義点でのdT/du[kg.K2/J]

内部で使用するパラメータとしてscale、Tm1、Tm2がある。

protected
  parameter Real scale=0.999;
  parameter Modelica.SIunits.Temperature Tm1=TSol + (1 - scale)*(TLiq - TSol);
  parameter Modelica.SIunits.Temperature Tm2=TSol + scale*(TLiq - TSol);

関数のアルゴリズムは以下の通り。assert文でBuildings.HeatTransfer.Conduction.nSupPCM==6とTLiq>TSolを要求している。

algorithm 
  assert(Buildings.HeatTransfer.Conduction.nSupPCM == 6,
    "The material must have exactly 6 support points for the u(T) relation.");
  assert(TLiq > TSol, "TLiq has to be larger than TSol.");
  // Get the derivative values at the support points
  ud:={c*scale*TSol,
       c*TSol,
       c*Tm1 + LHea*(Tm1 - TSol)/(TLiq - TSol),
       c*Tm2 + LHea*(Tm2 - TSol)/(TLiq - TSol),
       c*TLiq + LHea,
       c*(TLiq + TSol*(1 - scale)) + LHea};
  Td:={scale*TSol,
       TSol,
       Tm1,
       Tm2,
       TLiq,
       TLiq + TSol*(1 - scale)};
  dT_du := Buildings.Utilities.Math.Functions.splineDerivatives(
      x=ud,
      y=Td,
      ensureMonotonicity=ensureMonotonicity);

近似式に使用する6つの点が何なのかがここで分かった。⊿T=TSol*(1 - scale)と表すと6つの温度が表す点は以下の通り。そしてそれぞれの温度の時の比内部エネルギーudと微分dT_duが定義されている。

要素番号 Td
1 固相線-⊿T
2 固相線
3 固相線+⊿T
4 液相線-⊿T
5 液相線
6 液相線+⊿T

関数splineDerivativesの中身

dT_duを計算しているBuildings.Utilities.Math.Functions.splineDerivativesという関数について見ていく。この関数は三次エルミートスプライン補間に用いる微分を求める関数であり、任意の長さの配列に対して使用することができる。

入力

名称 デフォルト値 説明
Real x[:] - 横軸の値
Real y[size(x, 1)] - 横軸の値
Boolean ensureMonotonicity isMonotonic(y, strict=false) Trueの時に強制的に3次エルミートスプラインを単調にする

出力

名称 説明
Real d[:] 微分の値

内部で使用する変数として配列の長さn=size(x, 1)、傾きdelta[n - 1]、係数alpha、beta、tauがある。
関数のアルゴリズムは以下の通り。assert文にでてくるBuildings.Utilities.Math.Functions.isMonotonicという関数は与えられた配列a[:]が単調増加または単調減少の場合にtrueを返す関数で、strict=falseの場合にa[n]=a[n+1]となる部分があっても(途中で平らな部分があっても)trueを返すようになっている。Buildings.Utilities.Math.Functions.を省略しているのでModelica標準の関数なのか少しややこしい。ここではxの配列は常に増加しなくてはならず、yはensureMonotonicity=falseならば値に制限がなく、ensureMonotonicity=trueの時は単調増加でも単調減少でもよく途中で同じ値が続くことも許容している。
微分の与え方については、まずn=1のときに微分は0、n=2のときにどちらの点の微分も二つの点を結んだ直線の傾きになる。それ以外ははじめに説明したとおりに実装されておりensureMonotonicity=trueのときに補正を加えるようになっている。

algorithm 
  if (n>1) then
    assert(x[1] < x[n], "x must be strictly increasing.
  Received x[1] = " + String(x[1]) + "
           x[" + String(n) + "] = " + String(x[n]));
  // Check data
    assert(isMonotonic(x, strict=true),
      "x-values must be strictly monontone increasing or decreasing.");
    if ensureMonotonicity then
      assert(isMonotonic(y, strict=false),
        "If ensureMonotonicity=true, y-values must be monontone increasing or decreasing.");
    end if;
  end if;

  // Compute derivatives at the support points
  if n == 1 then
    // only one data point
    d[1] :=0;
  elseif n == 2 then
    // linear function
    d[1] := (y[2] - y[1])/(x[2] - x[1]);
    d[2] := d[1];
  else
    // Slopes of the secant lines between i and i+1
    for i in 1:n - 1 loop
      delta[i] := (y[i + 1] - y[i])/(x[i + 1] - x[i]);
    end for;
    // Initial values for tangents at the support points.
    // End points use one-sided derivatives
    d[1] := delta[1];
    d[n] := delta[n - 1];

    for i in 2:n - 1 loop
      d[i] := (delta[i - 1] + delta[i])/2;
    end for;

  end if;
  // Ensure monotonicity
  if n > 2 and ensureMonotonicity then
    for i in 1:n - 1 loop
      if (abs(delta[i]) < Modelica.Constants.small) then
        d[i] := 0;
        d[i + 1] := 0;
      else
        alpha := d[i]/delta[i];
        beta := d[i + 1]/delta[i];
        // Constrain derivative to ensure monotonicity in this interval
        if (alpha^2 + beta^2) > 9 then
          tau := 3/(alpha^2 + beta^2)^(1/2);
          d[i] := delta[i]*alpha*tau;
          d[i + 1] := delta[i]*beta*tau;
        end if;
      end if;
    end for;
  end if;
end splineDerivatives;

Modelica.Constants.small=1.e-60であり、PCが表示する限界の桁数の小さな値という意味らしい。 前回出てきたのはModelica.Constants.eps=1.e-15でこちらは1.0+eps=1.0となる最大の数値だそうなのでそちらでも十分なような気がするというか使い分けがよく分からない。

PCMの比内部エネルギーと温度の計算

初期条件の与え方の確認が終わったので実際の計算部分を確認する。コネクタとのやり取りや分割した各点での熱流の計算は前回確認した通りである。
定常計算のときには比内部エネルギーの変化がなく両端の温度と熱抵抗から各点の温度を計算している。おそらく比内部エネルギーは温度から逆算されている。
非定常計算の場合にはder(u[i]) = (Q_flow[i]-Q_flow[i+1])*mInv[i]とあるように温度ではなく比内部エネルギーの変化量を熱流により定義している。そしてBuildings.HeatTransfer.Conduction.BaseClasses.temperature_uという関数に先ほどの3つの配列と比内部エネルギーを与えて温度を計算している。

equation 
    port_a.Q_flow = +Q_flow[1];
    port_b.Q_flow = -Q_flow[end];
    port_a.T-T[1]    = if stateAtSurface_a then 0 else Q_flow[1]*RNod[1];
    T[nSta]-port_b.T = if stateAtSurface_b then 0 else Q_flow[end]*RNod[end];
    for i in 1:nSta-1 loop
       T[i]-T[i+1] = Q_flow[i+1]*RNod[i+1];
    end for;

    // Steady-state heat balance
    if material.steadyState then
      for i in 2:nSta+1 loop
        Q_flow[i] = port_a.Q_flow;
      end for;

      for i in 1:nSta loop
        if material.phasechange then
          T[i]=Buildings.HeatTransfer.Conduction.BaseClasses.temperature_u(
                    ud=ud,
                    Td=Td,
                    dT_du=dT_du,
                    u=u[i]);
        else
          u[i]=0; 
        end if;
      end for;
    else
      if material.phasechange then
        for i in 1:nSta loop
          der(u[i]) = (Q_flow[i]-Q_flow[i+1])*mInv[i];
          T[i]=Buildings.HeatTransfer.Conduction.BaseClasses.temperature_u(
                    ud=ud,
                    Td=Td,
                    dT_du=dT_du,
                    u=u[i]);
        end for;
      else
        for i in 1:nSta loop
          der(T[i]) = (Q_flow[i]-Q_flow[i+1])*CInv[i];
        end for;
        for i in 1:nSta loop
          u[i]=0; // u is not required in this case
        end for;
      end if;
    end if;

関数temperature_uの中身

Buildings.HeatTransfer.Conduction.BaseClasses.temperature_uという関数は線形補間や3次エルミートスプライン補間により比内部エネルギーから温度を求める関数である。この関数では与えられた内部エネルギーがどの区間にいるのかを判定しているだけで、実際に温度を計算してるのはBuildings.Utilities.Math.Functions.cubicHermiteLinearExtrapolationという関数である。dT_duの要素数をnSupPCMに指定していない理由がわからないが、動いているプログラムをいじりたくはないのでもやもやする。

入力

名称 デフォルト値 説明
SpecificInternalEnergy ud[nSupPCM] - 定義点でのud[J/kg]
Temperature Td[nSupPCM] - 定義点でのTd[K]
Real dT_du[:] - 定義点でのdT/du[kg.K2/J]
SpecificInternalEnergy u - 比内部エネルギー[J/kg]

出力

名称 説明
Temperature T 温度[K]

関数のアルゴリズムは以下の通り。

algorithm 
  i := 1;
  for j in 1:size(ud,1) - 1 loop
    if u > ud[j] then
      i := j;
    end if;
  end for;
  T :=  Buildings.Utilities.Math.Functions.cubicHermiteLinearExtrapolation(
     x=u,
     x1=ud[i],
     x2=ud[i + 1],
     y1=Td[i],
     y2=Td[i + 1],
     y1d=dT_du[i],
     y2d=dT_du[i + 1]);

iは6つの点から隣り合う2つの点を選ぶインデックスであり、1から5のどれかになる。そして比内部エネルギーが6つの点の外側にあった場合には線形外挿されることになる。図で表すと以下のようになる。分かりやすさのため⊿Tを大きく書いたが実際は0.01K程度である。傾きの異なる二つの直線を滑らかにつなぐためにこのようにしているのだと思う。
f:id:kinonotofu:20190131012127p:plain

関数cubicHermiteLinearExtrapolationの中身

3次エルミートスプラインは与えられた点の内側の補間のみをするようになっているため、Buildings.Utilities.Math.Functions.cubicHermiteLinearExtrapolation関数では補間に用いる6つの点の外側で線形外挿する処理を行っている。3次エルミートスプラインはModelica.Fluid.Utilities.cubicHermiteという関数を使用して行っている。

入力

名称 デフォルト値 説明
Real x - 補間値を得たい横軸の値
Real x1 - 補間に用いる横軸の値
Real x2 - 補間に用いる横軸の値
Real y1 - 補間に用いる縦軸の値(x1に対応する値)
Real y2 - 補間に用いる縦軸の値(x2に対応する値)
Real y1d - 補間に用いる傾き(x1に対応する値)
Real y2d - 補間に用いる傾き(x2に対応する値)

出力

名称 説明
Real y 補間して得た縦軸の値
algorithm 
  if (x > x1 and x < x2) then
    y:=Modelica.Fluid.Utilities.cubicHermite(
      x=x,
      x1=x1,
      x2=x2,
      y1=y1,
      y2=y2,
      y1d=y1d,
      y2d=y2d);
  elseif x <= x1 then
    y:=y1 + (x - x1)*y1d;
  else
    y:=y2 + (x - x2)*y2d;
  end if;

関数cubicHermiteの中身

Modelica.Fluid.Utilities.cubicHermiteが3次エルミートスプラインに使用する関数である。これはModelica標準ライブラリの関数なのだがなぜがFluidパッケージの中にある。

入力

名称 デフォルト値 説明
Real x - 補間値を得たい横軸の値
Real x1 - 補間に用いる横軸の値
Real x2 - 補間に用いる横軸の値
Real y1 - 補間に用いる縦軸の値(x1に対応する値)
Real y2 - 補間に用いる縦軸の値(x2に対応する値)
Real y1d - 補間に用いる傾き(x1に対応する値)
Real y2d - 補間に用いる傾き(x2に対応する値)

出力

名称 説明
Real y 補間して得た縦軸の値

他に内部で使用する変数としてh、t、h00、h01、h10、h11などはじめに確認した式そのままの変数とaux3(tの三乗)やaux2(tの二乗)などの変数がある。アルゴリズムは以下の通り。

algorithm
  h := x2 - x1;
  if abs(h)>0 then
    // Regular case
    t := (x - x1)/h;
    aux3 :=t^3;
    aux2 :=t^2;
    h00 := 2*aux3 - 3*aux2 + 1;
    h10 := aux3 - 2*aux2 + t;
    h01 := -2*aux3 + 3*aux2;
    h11 := aux3 - aux2;
    y := y1*h00 + h*y1d*h10 + y2*h01 + h*y2d*h11;
  else
    y := (y1 + y2)/2;
  end if;
end cubicHermite;

これで一応すべて確認したことになる。

おわりに

PCMについては実用上はこれで使えるのだと思う。固相と液相の物性値を変えたい場合も比内部エネルギーのように温度との関係式を定義した関数を用意するとよいのかもしれない。ただし、いろいろと計算を細かくしていくと計算量がどんどん膨大になっていきそうな気もする。比内部エネルギーから温度を算出する関数を定義すれば、勝手に温度から比内部エネルギーの算出もしてくれるのはModelicaの便利なところだと思った。
OpenModelicaのバージョンが1.13.1がリリースされていた。replaceableへの対応に関しては1.14以降になったらしい。さらにModelica標準ライブラリも3.2.3が正式リリースになったようなので使ってみてもよいかもしれない。
次回からは窓の計算モデル(Buildings.HeatTransfer.Windows.Window)を見ながら構成を確認したり日射周りの計算などを見たりしていきたいと思う。

PyGemのインストールメモ

概要

PyGemはメッシュモーフィングのライブラリで1400万格子まで対応するらしい。以下の記事をみて面白そうだと思ったのでインストール方法だけメモ。
qiita.com

Linux環境の構築をサボっているので今回はColabにインストールしたけれど、普通の人は自分のLinux環境に入れればいいと思う。

はじめに

Colabに行き「ファイル-Python3の新しいノートブック」を選択。あとは以下のコードを入力してShift+Enterで実行していけばいい。前の処理が終わるのを待つ必要はない。

依存ライブラリのインストール

!pip install numpy scipy matplotlib vtk numpy-stl sphinx nose

Anacondaのインストール

%%bash
wget -c https://repo.continuum.io/archive/Anaconda3-5.1.0-Linux-x86_64.sh
bash ./Anaconda3-5.1.0-Linux-x86_64.sh -b -f -p /usr/local

Anacondaのアップデート

!conda update -y -n base conda -c defaults

OCCのインストール

!conda install -y -c conda-forge -c dlr-sc -c pythonocc -c oce pythonocc-core==0.18.1 python=3.6

PyGemのインストール

%%bash
wget https://github.com/mathLab/PyGeM/archive/v1.1.tar.gz
tar -zxvf v1.1.tar.gz
cd PyGeM-1.1
2to3 -w ./
python setup.py install

Tutorial1の実行

PyGemがインスト-ルできたので、STLファイルをFFDで変形するサンプルを実行する。

%matplotlib inline
import sys
sys.path.append('/usr/local/lib/python3.6/site-packages/')
sys.path.append('/usr/local/lib/python3.6/site-packages/pygem-1.1-py3.6.egg/')
from pygem import FFDParameters, FFD, StlHandler

params = FFDParameters()
params.read_parameters(filename='./PyGeM-1.1/tests/test_datasets/parameters_test_ffd_sphere.prm')

stl_handler = StlHandler()
mesh_points = stl_handler.parse('./PyGeM-1.1/tests/test_datasets/test_sphere.stl')
fig = stl_handler.plot(plot_file='./PyGeM-1.1/tests/test_datasets/test_sphere.stl')

free_form = FFD(params, mesh_points)
free_form.perform()
new_mesh_points = free_form.modified_mesh_points
stl_handler.write(new_mesh_points, 'test_sphere_mod.stl')
fig = stl_handler.plot(plot_file='test_sphere_mod.stl')

以下のように球体の変形ができている(上が元のSTL、下が変形後のSTL
f:id:kinonotofu:20190108164634p:plain

おわりに

PyGeM-1.1/pygem/params/の中に変形手法毎のファイルがありそれぞれのread_parametersクラスにどんなパラメータがあるか書いてある。OpenFOAMでの使い方はこの記事が分かりやすかった。ただし、昔は指定していたけれど現在は内部で計算するため@propertyをつけて読取専用としたプロパティがあるので注意。

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

概要

OpenModelicaでLBNLのBuildingsライブラリを使って熱伝導の計算を学ぶ。初心者なのでいろいろ教えて頂けるとありがたい。
すこし間があいてしまったが前回はBuildings.HeatTransferの全体も見てPCM(Phase Change Material、潜熱蓄熱材)を使った簡単な伝熱計算をしてみた。今回は1次元熱伝導計算の基本となるBuildings.HeatTransfer.Conduction.SingleLayerのソースを詳しく確認する。ただし、PCMの計算の細かいところは次回に行う。

使用バージョン
-Buildings 5.1.0

材料の物性値

まずは熱伝導に使用する材料の物性値について確認する。SingleLayerの材料の物性値はreplaceable parameter Data.BaseClasses.Material materialとして定義してあるのだがこう書くものなのだろうか。replaceable record ConMaterial = Data.BaseClasses.MaterialとしてConMaterial materialと書いたほうがよい気がするのだがどうなのだろう。recordにparameterの指定をするのもなくはないのかもしれないけれど、どうなのだろう。問題なく動いているなら下手に手を加えるわけにも行かないのかもしれない。Buildings.HeatTransfer.Data.BaseClasses.Materialに定義されている変数は以下の通り。

種類 名称 デフォルト値 説明
パラメータ Modelica.SIunits.Length x - 材料の厚み
パラメータ Modelica.SIunits.ThermalConductivity k - 熱伝導率
パラメータ Modelica.SIunits.SpecificHeatCapacity c - 比熱
パラメータ Modelica.SIunits.Density d - 密度
パラメータ Real R(unit="m2.K/W") - 熱抵抗
パラメータ Integer nStaRef(min=0) 3 参照部材(0.2mのコンクリート)の状態変数の数
パラメータ Integer nSta(min=1) max(1, integer(ceil(nStaReal))) 材料の実際の状態変数の数
パラメータ Boolean steadyState (c < Modelica.Constants.eps or d < Modelica.Constants.eps) Trueの時に定常熱伝導で計算するフラグ
パラメータ Real piRef 331.4 参照部材(0.2mのコンクリート)のx/sqrt(alpha)
パラメータ Real piMat if steadyState then piRef else xsqrt(cd/k) x/sqrt(alpha)
パラメータ Real nStaReal(min=0) nStaRef*piMat/piRef 実数としての状態の数
パラメータ Modelica.SIunits.Temperature TSol - 固相温度、PCMでのみ使用
パラメータ Modelica.SIunits.Temperature TLiq - 液相温度、PCMでのみ使用
パラメータ Modelica.SIunits.SpecificInternalEnergy LHea - 状態変化の潜熱
定数 Boolean ensureMonotonicity false Trueの時に微分dT/duを強制的に単調にする
定数 Boolean phasechange false 材料がPCMの時にTrueにする

steadyStateのところのModelica.Constants.epsは極小値をあらわす定数で1.e-15である。比熱か密度が0のときに定常としている。また、TSol以降はPCM用の設定で最後の二つの変数の修飾子はparameterではなくconstantになっている。

熱伝導計算に実際に使用するBuildings.HeatTransfer.Data.Solidsの一般的なものは以下のようになっていて、熱抵抗の式とPCM関係の設定をfinalで変更できないようにしている。

record Generic "Thermal properties of solids with heat storage"
    extends Buildings.HeatTransfer.Data.BaseClasses.Material(final R=x/k,
                                                             final TSol=293.15,
                                                             final TLiq=293.15,
                                                             final LHea=0,
                                                             final phasechange=false);
end Generic;

これに具体的に値を与えるときは以下のようにする。実際に使用するときはこの書き方だけ分かっていればよいと思う。

record Concrete = Buildings.HeatTransfer.Data.Solids.Generic (
    k=1.4,
    d=2240,
    c=840) "Concrete (k=1.4)";

材料の分割数の決定

Buildings.HeatTransfer.Data.SolidsにnStaRef、nSta、piRef、piMat、nStaRealあたりの説明があり、非定常熱伝導の現象の無次元数であるフーリエ数Fo(伝熱量/畜熱量)を使用して材料の物性値と厚みにより分割数を決定している。
フーリエ数は以下のように定義される。
f:id:kinonotofu:20181226102546p:plain

αは熱拡散係数、tは特徴時間、Lは特徴長さである。任意の時間に対して分割数を決定するので、次のようなフーリエ数を用いた定数Πを導入する。
f:id:kinonotofu:20181226102628p:plain

ここで以下のように定数が一致するように各素材の分割数N'(=nStaReal)を決定する。Πx(=piMat)が分割数を求める部材の定数、Nref(=nStaRef=3)とΠref(=piRef=331.4)が参照部材の値。上の式でL=x/Nとして等式を立てていると考えると分かりやすいと思う。
f:id:kinonotofu:20181226102645p:plain

参照部材は0.2mのコンクリート(α=3.64E-7m/s)としており、デフォルトの分割数は3になっている。伝熱現象の解像度がこれと等しくなるようにその他の材料の分割数を決定している。実際にはN'の小数点以下切り上げとしたNx(=nSta=max(1, integer(ceil(nStaReal))))を使用しているのでもう少し細かくなる。なお、実際には等間隔に分割されるとは限らない。また、ここでの分割数が1であり熱伝導の計算で表面に節点を定義する場合は分割数を2として計算する。詳細は熱伝導のところで触れる。定常状態では分割数は関係なくなると思うが、デフォルトでnStaRefと同じ値3になる。

PartialConductorの継承

SingleLayerはPartialConductorを継承している。PartialConductorの内容は以下のシンプルなものでコネクタはHeatPort二つ、面積Aと熱抵抗Rはパラメータで後から何かしら値を与えなければならず、熱コンダクタンスUAや熱貫流率Uはfinalで定義されている。また、equationにてdTが二つのポート間の温度差として定義されている。定常でも非定常でもよいようにするためか、熱容量に関する値は定義されていない。

partial model PartialConductor
  extends Buildings.BaseClasses.BaseIcon;
  parameter Modelica.SIunits.Area A;
  parameter Modelica.SIunits.ThermalResistance R;
  final parameter Modelica.SIunits.ThermalConductance UA = 1/R;
  final parameter Modelica.SIunits.CoefficientOfHeatTransfer U = UA/A;
  Modelica.SIunits.TemperatureDifference dT "port_a.T - port_b.T";
  Modelica.Thermal.HeatTransfer.Interfaces.HeatPort_a port_a;
  Modelica.Thermal.HeatTransfer.Interfaces.HeatPort_b port_b;
equation 
  dT = port_a.T - port_b.T;
end PartialConductor;

そしてPartialConductorはSingleLayerで継承される際にRを材料に定義された値として指定されている。

  extends Buildings.HeatTransfer.Conduction.BaseClasses.PartialConductor(
   final R=if (material.R < Modelica.Constants.eps) then material.x/material.k/A else material.R/A);

Modelica.Constants.epsは1.e-15なのでこのif文はRの値が0のときには 厚み/熱伝導率/面積、それ以外のときは熱抵抗/面積となるようにしている。Aは先ほど見たようにPartialConductorの持つ変数であり、PartialConductorに与える値の定義にそのまま記述できるらしい。材料の方のRは[m2.K/W]でPartialConductorのRは[K/W]]なので面積で除算をしている。

SingleLayerの計算定義点の設定

計算上の層の分割がどのようになっているのかを確認する。

まず表面に状態を定義するかどうかを指定する変数がある。デフォルトでは両面とも定義するようになっている。stateAtSurface_a=trueとすればT[1]からA面の表面温度を取得できるが、falseの場合はport_a.Tから表面温度を取得するようになる。この変数によって材料の分割位置がかわってくる。

  parameter Boolean stateAtSurface_a=true;
  parameter Boolean stateAtSurface_b=true;

それから分割数に関わる変数がいくつかある。nStaが温度節点の数、nR=nSta+1は熱抵抗の数になる。nSta2が先ほど見た部材で定義した分割数である。定義しなおすとDymolaでエラーが出るらしいがfinalはつけられないのだろうか。表面に温度節点を定義した場合には最低でも2分割になる。

  parameter Integer nSta2=material.nSta;
protected 
  final parameter Integer nSta=
    max(nSta2,
        if stateAtSurface_a or stateAtSurface_b then 2 else 1);
  final parameter Integer nR=nSta+1;

さらに、SingleLayerで定義されている各節点の状態値は以下のとおり。温度(PCMは内部エネルギーも)と熱流が変数としてあり、それを計算するためのパラメータがprotectedの部分に書かれている。protectedとすると例えば他のモデルでSingleLayer sLayを定義したときにsLay.RNod[1]などとして値を参照することができなくなる。Q_flow[nSta+1]はQ_flow[nR]にすればよいのではないかと思う。

種類 名称 説明
変数 Modelica.SIunits.Temperature T[nSta] 温度
変数 Modelica.SIunits.HeatFlowRate Q_flow[nSta+1] 熱流
変数 Modelica.SIunits.SpecificInternalEnergy u[nSta] 比内部エネルギー(PCM用)
パラメータ Modelica.SIunits.ThermalResistance RNod[nR] 熱抵抗、nR=nSta+1
パラメータ Modelica.SIunits.Mass m[nSta] 質量
パラメータ Real mInv[nSta] 質量の逆数
パラメータ Modelica.SIunits.HeatCapacity C[nSta] 熱容量
パラメータ Real CInv[nSta] 熱要路湯の逆数

これらの変数が層のどの範囲を受け持っているのかを確認する。下の方に図にまとめたのでそちらを先に参照したほうがよいかもしれない。
まずは熱抵抗の設定、それぞれの表面に温度節点を定義しているのかどうかで変わり、表面に節点があるときは材料と反対側の抵抗は使わないので0にしている。それ以外は表面が半分の抵抗になることで変わらない。

  parameter Modelica.SIunits.ThermalResistance RNod[nR]=
    if (stateAtSurface_a and stateAtSurface_b) then 
      if (nSta==2) then 
        {(if i==1 or i==nR then 0 else R/(nSta-1)) for i in 1:nR}
      else 
        {(if i==1 or i==nR then 0 elseif i==2 or i==nR-1 then R/(2*(nSta-2)) else R/(nSta-2)) for i in 1:nR}
    elseif (stateAtSurface_a and (not stateAtSurface_b)) then 
        {(if i==1 then 0 elseif i==2 or i==nR then R/(2*(nSta-1)) else R/(nSta-1)) for i in 1:nR}
    elseif (stateAtSurface_b and (not stateAtSurface_a)) then 
       {(if i==nR then 0 elseif i==1 or i==nR-1 then R/(2*(nSta-1)) else R/(nSta-1)) for i in 1:nR}
    else 
      {R/(if i==1 or i==nR then (2*nSta) else nSta) for i in 1:nR};

次に温度の設定。指定した初期温度で定常条件での温度を指定している。先ほど定義した抵抗の値を使っているので上の図の位置に温度節点があると考えられる。T_a_startとT_b_startは温度の初期値を決めるパラメータで293.15「K](20℃)がデフォルト値である。cateは配列をつなげて第一引数で指定した次元の配列を作る関数なのだが、ここの数式は場合分けしているものの、別にelse以降の式だけで問題ないように思うのだけれど違うのだろうか。(stateAtSurface_a=trueならばRNod[1]=0)

  Modelica.SIunits.Temperature T[nSta](start=
   if stateAtSurface_a then 
     cat(1,
       {T_a_start},
       {(T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i-1)) for i in 2:nSta})
   else 
    {(T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i)) for i in 1:nSta},
   each nominal=300);

次に熱容量関係の設定。ここは質量mを設定して残りのmInv、C、CInvはmを演算しているだけである。面積×厚み×密度を分割しており表面に温度節点を定義していれば表面と表面から2番目の層が他の半分の厚さになっている。

  parameter Modelica.SIunits.Mass m[nSta]=
   (A*material.x*material.d) *
   (if (stateAtSurface_a and stateAtSurface_b) then 
      if (nSta==2) then 
        {1/(2*(nSta-1)) for i in 1:nSta}
      elseif (nSta==3) then 
        {1/(if i==1 or i==nSta then (2*(nSta-1)) else (nSta-1)) for i in 1:nSta}
      else 
        {1/(if i==1 or i==nSta or i==2 or i==nSta-1 then (2*(nSta-2)) else (nSta-2)) for i in 1:nSta}
    elseif (stateAtSurface_a and (not stateAtSurface_b)) then 
      {1/(if i==1 or i==2 then (2*(nSta-1)) else (nSta-1)) for i in 1:nSta}
    elseif (stateAtSurface_b and (not stateAtSurface_a)) then 
      {1/(if i==nSta or i==nSta-1 then (2*(nSta-1)) else (nSta-1)) for i in 1:nSta}
    else 
      {1/(nSta) for i in 1:nSta});

残りの定義は以下の通り。熱容量は質量×比熱である。

  final parameter Real mInv[nSta]=
    if material.steadyState then zeros(nSta) else {1/m[i] for i in 1:nSta};
  final parameter Modelica.SIunits.HeatCapacity C[nSta] = m*material.c;
  final parameter Real CInv[nSta]=
    if material.steadyState then zeros(nSta) else {1/C[i] for i in 1:nSta};

分かりにくいので図にすると以下のようになる。表面温度を設定した部分では熱容量の配分と熱抵抗の配分が一致していないような気がする。説明にもそうしていると書いてあるのでバグではないようなのだけれど少しもやもやする。熱抵抗の中心で熱容量を分割したほうがモデルとしてわかりやすいのではないだろうか。
f:id:kinonotofu:20181226102717p:plain

initial equation

assert文があり、分割時の熱抵抗と熱容量の合計値が間違っていないかをチェックしている。
温度は変数を定義する際にstartで値を指定しても結局initial equationを定義してあげないと警告がでるようなので、再度定義をしていると思われる。steadyStateInitialは初期化条件を両端の温度から定常条件で与えるかを選ぶ変数でtrueの時に全ての節点の温度の時間微分が0という条件となる。そうでない場合は先ほどと同じである。この違いはsteadyStateInitial=trueのときには両端の値はコネクタから受け取ったものになるのだろうか。
内部エネルギーuの定義はないのだけれど、Modelica.SIunits.SpecificInternalEnergy u[nSta](each nominal=270000)の270000を使用するのだろうか。今回は使用しないが、PCMであっても特に定義していないように思う。何を定義するのかは独立変数の数だけにしないといけないのかもしれない。ud、Td、dT_duはPCM以外では使用しないので0が入っている。これらはparameterなのだけれどまとめて関数で値を設定するのでこちらで定義していると思われる。

  assert(abs(sum(RNod) - R) < 1E-10, "Error in computing resistances.");
  assert(abs(sum(m) - A*material.x*material.d) < 1E-10, "Error in computing mass.");

  // The initialization is only done for materials that store energy.
    if not material.steadyState then
      if steadyStateInitial then
        if material.phasechange then
          der(u) = zeros(nSta);
        else
          der(T) = zeros(nSta);
        end if;
      else
        if stateAtSurface_a then
          T[1] = T_a_start;
          for i in 2:nSta loop
            T[i] =T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i-1);
          end for;
        else // stateAtSurface_a == false
          for i in 1:nSta loop
            T[i] = T_a_start + (T_b_start - T_a_start)*UA*sum(RNod[k] for k in 1:i);
          end for;
        end if;
      end if;
    end if;

   if material.phasechange then
     (ud, Td, dT_du) = Buildings.HeatTransfer.Conduction.BaseClasses.der_temperature_u(
       c =  material.c,
       TSol=material.TSol,
       TLiq=material.TLiq,
       LHea=material.LHea,
       ensureMonotonicity=material.ensureMonotonicity);
   else
     ud    = zeros(Buildings.HeatTransfer.Conduction.nSupPCM);
     Td    = zeros(Buildings.HeatTransfer.Conduction.nSupPCM);
     dT_du = zeros(Buildings.HeatTransfer.Conduction.nSupPCM);
   end if;

equation

コネクタから熱流を受け取る(熱流の向きはi→i+1の方向)

    port_a.Q_flow = +Q_flow[1];
    port_b.Q_flow = -Q_flow[end];

コネクタから温度を受け取る、表面に節点があればそのまま、なければ定常計算で端の節点の温度を計算する。

    port_a.T-T[1]    = if stateAtSurface_a then 0 else Q_flow[1]*RNod[1];
    T[nSta]-port_b.T = if stateAtSurface_b then 0 else Q_flow[end]*RNod[end];

定常のときはT[:]の計算、非定常のときはQ_flow[:]の計算になる。両端はコネクタから受け取った値である。

    for i in 1:nSta-1 loop
       T[i]-T[i+1] = Q_flow[i+1]*RNod[i+1];
    end for;

定常計算のときは全ての熱流が等しいという条件で上のQ_flow[:]の計算Q_flow[i] = port_a.Q_flow、非定常のときはT[:]の計算der(T[i]) = (Q_flow[i]-Q_flow[i+1])*CInv[i]。比内部エネルギーは使わないので0を与えている。PCMの計算については次回確認してい。
表面に節点を定義としたところにFixedTemperatureを直接つなぎ込むと、その固定温度とここで計算する温度が一致させられずエラーになるので注意が必要である。

    // Steady-state heat balance
    if material.steadyState then
      for i in 2:nSta+1 loop
        Q_flow[i] = port_a.Q_flow;
      end for;

      for i in 1:nSta loop
        if material.phasechange then
          // Phase change material
          T[i]=Buildings.HeatTransfer.Conduction.BaseClasses.temperature_u(
                    ud=ud,
                    Td=Td,
                    dT_du=dT_du,
                    u=u[i]);
        else
          // Regular material
          u[i]=0; // u is not required in this case
        end if;
      end for;
    else
      // Transient heat conduction
      if material.phasechange then
        // Phase change material
        for i in 1:nSta loop
          der(u[i]) = (Q_flow[i]-Q_flow[i+1])*mInv[i];
          // Recalculation of temperature based on specific internal energy
          T[i]=Buildings.HeatTransfer.Conduction.BaseClasses.temperature_u(
                    ud=ud,
                    Td=Td,
                    dT_du=dT_du,
                    u=u[i]);
        end for;
      else
        // Regular material
        for i in 1:nSta loop
          der(T[i]) = (Q_flow[i]-Q_flow[i+1])*CInv[i];
        end for;
        for i in 1:nSta loop
          u[i]=0; // u is not required in this case
        end for;
      end if;
    end if;

おわりに

とりあえず単層壁の非定常熱伝導の計算をみていった。計算モデルの内容はしっかり分かったのだが、1次元の熱伝導計算もこだわるといろいろでてくるのだなぁと思った。
層の分割数の決定方法については勉強になったが、結局分割幅が表面の節点の定義状況によって変わってしまうので本当にこれでいいのだろうかという気になった。それから現状では場合分けを理解するのに少し苦労するのではないかと思う。まず熱容量を等分割して熱容量のセル中心に温度節点を配置、熱抵抗分割位置に対応するように配置し、表面に温度節点があるときは端のセルの半分の熱容量を表面に与えるという仕組みだった。このとき熱抵抗の定義はかわらないため、熱容量の扱いのみが変わることになる。
次回は同じモデルの中のPCMの計算部分について見ていく。PCMでは相変化する際にはスプライン補間を使いながら内部エネルギーから温度を読み取るようになっている。このやり方は熱水分同時移動の物性値などにも転用できるのではないかなぁと思うけれどどうなんだろう。次回のModelica勉強会のネタにできればいいのだけれどどうなることやら。次回は来年になりそうだけれどマイペースに頑張ります。

Google ColaboratoryでJModelica

概要

Google Colaboratoryは機械学習の研究や教育用途であれば無償で使用できるツールで、Googleのアカウントがあればブラウザを使用してJupyter Notebookをクラウド上で動かすことができる。JModelicaが使用できたからといって何かいいことがあるかはよく分からないけれど、以下の記事でOpenFOAMを動かしているのを見て面白いと思ったのでJModelicaでもやってみた。時間帯にもよるかもしれないが、けっこう時間がかかった(30~40分くらい)ので決して手軽ではない。

qiita.com

注意点

基本的にはユーザーガイドのインストールの方法でよいが以下の部分だけ注意。

  • Python2のノートブックを使う。
  • ホームディレクトリ(~)が/rootになっているが初期位置の/contentでビルドした。
  • apt-getのときにfileとjccもインストールする。
  • JModelica.org/external/Assimuloをダウンロードしてくる。

JModelicaをビルドしてサンプルを実行

まずはGoogle ColaboratoryGoogleのアカウントにログインした状態で「ファイル→Python2の新しいノートブックを開く」で新規ファイルを作成する。あとは以下のコードのまとまりをそれぞれ入力してShift+Enterを押していけばよいだけである。前の処理が終わって無くてもShift+Enterを押しておけば順番に逐次実行していってくれる。
JModelicaをビルドしたらサンプルとしてBuildingsライブラリのBuildings.Examples.Tutorial.SpaceCooling.System3(室温をオンオフ制御して空調するもの)を実行してみる。

ビルドの環境を整える

%%bash
apt-get update
apt-get -y install g++ subversion gfortran ipython cmake swig ant openjdk-8-jdk python-dev python-numpy python-scipy python-matplotlib cython python-lxml python-nose python-jpype zlib1g-dev libboost-dev jcc file

Ipoptのビルド(JModelicaのビルドの前に必要)

%%bash
mkdir /content/Ipopt
cd /content/Ipopt
wget http://www.coin-or.org/download/source/Ipopt/Ipopt-3.12.12.tgz
tar xzf Ipopt-3.12.12.tgz
cd Ipopt-3.12.12/ThirdParty/Blas
./get.Blas
cd ../Lapack
./get.Lapack
cd ../Mumps
./get.Mumps
cd ../Metis
./get.Metis
cd ../../
mkdir build
cd build
../configure --prefix=/content/Ipopt/
make install

JModelicaのビルド

%%bash
export JAVA_HOME=/usr/lib/jvm/java-8-openjdk-amd64
export LD_LIBRARY_PATH=$JAVA_HOME/jre/lib/amd64:$JAVA_HOME/jre/lib/amd64/server
mkdir -p /content/JModelica
cd /content/JModelica
svn co https://svn.jmodelica.org/trunk JModelica.org
svn co https://svn.jmodelica.org/assimulo/trunk JModelica.org/external/Assimulo
cd JModelica.org
mkdir build
cd build
../configure --prefix=/content/JModelica --with-ipopt=/content/Ipopt/
make install
make casadi_interface

Buildingsライブラリの用意

%%bash
mkdir -p /content/JModelica/work
cd /content/JModelica/work
wget https://github.com/lbl-srg/modelica-buildings/releases/download/v5.1.0/Buildings-v5.1.0.zip
unzip Buildings-v5.1.0.zip
mv "Buildings 5.1.0" Buildings

JModelicaの計算設定ファイルの書き出し

with open('/content/JModelica/work/settings.py', 'w') as f:
    f.write('import os\n')
    f.write('from pymodelica import compile_fmu\n')
    f.write('from pyfmi import load_fmu\n')
    f.write('os.chdir(\'/content/JModelica/work\')\n')
    f.write('fmufile = compile_fmu(\'Buildings.Examples.Tutorial.SpaceCooling.System3\',\'Buildings\')\n')
    f.write('model = load_fmu(\'Buildings_Examples_Tutorial_SpaceCooling_System3.fmu\')\n')
    f.write('opts = model.simulate_options()\n')
    f.write('opts[\'ncp\'] = 144\n')
    f.write('opts[\'solver\'] = \'Radau5ODE\'\n')
    f.write('opts[\'Radau5ODE_options\'][\'rtol\'] = 1.0e-6 \n')
    f.write('model.simulate(start_time=15552000.,final_time=15638400., options=opts)\n')

計算の実行

%%bash
/content/JModelica/bin/jm_ipython.sh /content/JModelica/work/settings.py

計算結果のダウンロード

from google.colab import files
files.download('/content/JModelica/work/Buildings_Examples_Tutorial_SpaceCooling_System3_result.mat')

以上で結果ファイルをダウンロードできるので自分のPCで結果を見る。

f:id:kinonotofu:20181212004921p:plain

問題なく計算ができていそうである。

ファイルの入力方法

今回はpythonで設定ファイルを書き出したが、実際はファイルを用意してアップロードする方が楽だと思われる。以下を実行するとファイル選択ダイアログを開くことができる。

from google.colab import files
uploaded = files.upload()

または以下を実行して出たリンクからアカウント認証をしてコードを取得し、入力することでGoogleドライブのマウントができる。"/content/drive/My Drive"にマウントされるので自分のGoogleドライブに置いたファイルを使用することができるようになる。

from google.colab import drive
drive.mount('/content/drive')

時間の把握

時間は以下のようにして把握する。これが0.5daysを超えるとインスタンスが落とされるらしい。

!cat /proc/uptime | awk '{print $1 /60 /60 /24 "days (" $1 "sec)"}'

おわりに

Google Colaboratoryはセッションをはじめに立ち上げてから12時間、またはPCのスリープやブラウザを閉じるなどでセッションが切れて90分のどちらかでインスタンスが落ちるので注意が必要である。
Buildingsライブラリのユーザーガイドに熱流体を扱う場合のソルバーはDymolaだとデフォルトの残差1e-4のdasslでなく残差1e-6のradauがよいと書いてあったので、JModelicaでRadau5ODEにしてみた。別にデフォルトでも計算は問題なくできる。ちなみにJModelicaのデフォルトのソルバーは相対残差1e-4のCVodeというものになっている。アウトプットの時間間隔は10分になるように回数を144に設定した。
ポスト処理は自分のPC上で行ったが、matファイルはたぶんpythonで処理できると思うので、Google Colaboratory上でできるのではないかと思う。とりあえずJModelicaは動いたのでFMIを使った計算もできると思われる。Google Colaboratoryの本来の目的の機械学習との連携についてはまだよくわからない。

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

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