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勉強会のネタにできればいいのだけれどどうなることやら。次回は来年になりそうだけれどマイペースに頑張ります。