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

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

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関係の開発がさらに活発化するかもしれない。