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%混ぜた板材の熱量と温度の関係のモデルはだいたい以下のようになる。
この図は計算モデルを表しているものであり、比熱は固相、液相、固液混相域で同じ値を用いたものであるが、縦軸の比エンタルピー（熱量）と横軸の温度の関係を見ると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の値は以下のようになる。
そしてこれにアフィン変換（今回は回転やせん断はなく平行移動と拡大縮小のみ）を施して一般的な形にすると以下のようになる。h00～h11の係数は上の式と同じである。
この計算モデルがModelica標準ライブラリで用意されている。微分dy/dxについては自分で用意する必要があるのだが、今回の実装では⊿を各点を結んだ直線の傾きとすると、各点につながる直線の傾きの平均を各点の微分としている。
両端については片側にしか直線がないのでその値を使用する。
さらに補間に使用する微分の値に補正をかけることで単調性を確保することができる。つまり与えられた点が単調増加や単調減少なら補間した値も単調性を保つようにできるのだ。こちらはBuildingsライブラリで実装されており、以下の論文を参考にしているらしい。
F.N. Fritsch and R.E. Carlson, Monotone piecewise cubic interpolation. SIAM J. Numer. Anal., 17 (1980), pp. 238-246
論文は有料だが同著者の似たような記事でやや古いPDFがあるようだ。また、英語のWikipediaのMonotone cubic interpolationもこの内容である。Wikipediaの図を見れば補正がなぜ必要になるのかがわかると思う。
各点の微分を修正していくのだが、各点をつなぐ直線を左から調べて、両端の点の微分を順次修正していくため、右側の点は次の直線を調べたときに再度微分の値が修正されることになる。
まず直線の傾きが0になる場合はその両端の点の微分は0であるとする。
それ以外の場合はまず、単調性を確保させる際に発散させないように以下の条件を調べる。
今回の実装ではこれが満たされない区間は補正をしないだけになり警告がでないため、必ず単調性が確保されているわけではないかもしれない。あらかじめ出てきた曲線を調べておいたほうがよさそうである。条件を満たす場合には以下のような処理を行う。
このような処理を繰り返していけば単調性が確保されるらしい。証明までは踏み込まなかったけれど、とりあえずこのようにして計算できるようだ。
初期条件の定義と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程度である。傾きの異なる二つの直線を滑らかにつなぐためにこのようにしているのだと思う。
関数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)を見ながら構成を確認したり日射周りの計算などを見たりしていきたいと思う。
