International Essays |
[このノートブックは以下のwolfram blogをLLMツールにより日本語に翻訳したものです:Wolfram blog,How Long to Boil an Egg? FEM Modeling with Wolfram Language - 2025年9月15日
https://blog.wolfram.com/2025/09/15/how-long-to-boil-an-egg-fem-modeling-with-wolfram-language/]
https://blog.wolfram.com/2025/09/15/how-long-to-boil-an-egg-fem-modeling-with-wolfram-language/]
卵をどれくらい茹でればよいですか?Wolfram言語によるFEMモデリング
卵をどれくらい茹でればよいですか?Wolfram言語によるFEMモデリング
2025年9月15日
Ricardo López, Physics Applications Developer Intern, Algorithms R&D; Gay Wilson, Food Data Curator, Wolfram|Alpha Scientific Content; Oliver Ruebenkoenig, Senior Kernel Developer, Algorithms R&D
朝食好きの方が気になる質問,「卵をどれくらい茹でればよいですか?」にお答えする時が来ました.とても簡単に思えるかもしれません——卵を沸騰したお湯に入れて待つだけ——ですが,完全に固ゆで卵だけが美味しくタンパク質を摂取する唯一の方法と言うのは正確ではありません.有限要素法(FEM)を利用して,水中の卵の状態をシミュレーションし,卵内部の温度変化を評価することで,理想的な温度と時間を見つけて完璧な卵を作ることができます.そして,半熟の黄身や,しっかり固まったパサパサの黄身など,さまざまな仕上がりまでにかかる時間を予測できます.
複雑な物理システムのモデリングは,しばしばそれらを管理しやすい部分に分解することから始まります.これはまさにFEMが行うことです.FEMは,複雑な構造をより小さく単純な要素に分割し,支配方程式を局所的に解き,その結果を組み合わせて全体のシステムの挙動を捉えます.
このアプローチにより,物理的なプロトタイプのコストや労力をかけずに,実際のシナリオをシミュレーションおよび解析することが可能になります.ここでは,Wolfram言語 における有限要素解析のワークフローを紹介し,基本を網羅しつつ,よく発生する問題点を取り上げる簡単な例を順を追って説明します.
なぜ有限要素法を使いますか?
なぜ有限要素法を使いますか?
多くのPDE(たとえば古典的な方程式であるPoisson方程式やSchrödinger equationsなど)は,特に方程式を解こうとする領域が複雑な場合,解析解が常に存在するとは限りません.FEMはPDEを数値的に解くための有用な手法であり,その理由は少なくとも2つあります.第一に,任意の形状の領域でPDEを解くことができます.第二に,Laplace equationからNavier–Stokes equationsまで,さまざまな種類の微分方程式を解くことができます.
有限要素解析に必要なツールです
有限要素解析に必要なツールです
◼
領域はメッシュに離散化されます.このメッシュは,要素と呼ばれる多数の小さくて単純な部分領域から構成されます.FEM(有限要素法)の重要な考え方は,それぞれの要素内でPDE(偏微分方程式)のより単純なバージョンを解き,これらの局所解を組み合わせて,元のPDEの挙動を領域全体で近似する大域的な解を得ることです.
◼
解きたい特定のPDEです.
◼
境界条件は,PDEを解が計算されている領域の外部の世界と結び付けます.
Wolfram言語では,FEMは関数,およびを通じて実装されております.ゆで卵のシナリオではNDSolveValueを使用いたします.
バージョン1:ワークフローを開始します
バージョン1:ワークフローを開始します
モデルを一歩ずつ構築していき,いくつかのバージョンを経て徐々にモデルの複雑さを高めていきます.まず,物体の大きさや単位に注意しながら,シンプルなモデルから始めます.
最初のステップは,構造物のジオメトリまたは領域を定義することです.卵は軸対称性を持っているため,卵の断面を解析することでその挙動を十分に説明できると仮定します.FEMを行うためには,この領域をメッシュに離散化する必要があります.また,二次元メッシュは三次元メッシュほど多くの要素を必要としないため,ジオメトリの近似に必要な計算時間も短くなります.これが軸対称モデルを用いる重要な理由の一つです.
平均的な卵は直径約5 cm,半径約2.5 cmである場合が多いです.モデルのパラメータにはSI単位を使用するのが最善の方法ですので,卵の半径を0.025 mとして定義します.この時点では,意図的に黄身と卵白の素材の違いを無視しています.私たちの目的は,この初期バージョンのモデルをできるだけ単純にすることです.
有限要素パッケージをロードします:
In[]:=
Needs["NDSolve`FEM`"]
ここでは, を0に設定します.これにより,セッション内に保存される過去の結果の数がゼロに制限され,計算用メモリを最小限のコストで節約できます.
$HistoryLength を0に設定します:
In[]:=
$HistoryLength=0;
卵のシンプルな幾何学形状を作成します:
In[]:=
radius=0.025;Ω=RegionUnion[Circle[{0,0},radius,{-π/2,π/2}],Line[{{0,radius},{0,-radius}}]]
Out[]=
要素メッシュを生成します:
In[]:=
mesh=ToElementMesh[Ω]
Out[]=
ElementMesh[{{0.,0.025},{-0.025,0.025}},{TriangleElement[<492>]}]
有限要素メッシュを可視化します:
In[]:=
mesh["Wireframe"]
Out[]=
Wolfram言語でPDEをモデル化するためには,「PDEコンポーネント」と呼ばれるものを使用します.これらはPDEを指定するための構成要素です.これらの関数は,変数やパラメータを受け取り,PDE演算子を返します.その演算子はNDSolveValueと一緒に使用できます.特定の分野に対応するすべてのPDEコンポーネントの一覧にアクセスできます.卵内部の温度変化を調べるために,
求めている解,すなわち従属変数は温度(ケルビン)です.円筒対称性の場合,空間的な独立変数が2つあります:r は半径方向の座標(中心軸からの距離)であり, zは軸方向の座標(円筒の軸に沿った距離)です.どちらもメートルで表されます.さらに,卵の沸騰は時間に依存するため,時間変数t (秒)も必要です.
T
従属変数と独立変数のモデルを設定します:
In[]:=
vars={T[t,r,z],t,{r,z}};
また,パラメータも必要です.最初の簡略化したバージョンでは,領域の対称性のみを軸対称として指定します.これにより,注目している卵の断面に対して正しい項が得られます.HeatTransferPDEComponentは,残りの必要なパラメータをデフォルト値で補完しますので,今はそれで問題ありません.
パラメータを設定します:
In[]:=
pars=<|"RegionSymmetry"->"Axisymmetric"|>;
PDE演算子を設定します:
In[]:=
op=HeatTransferPDEComponent[vars,pars]
Out[]=
-,0.T[t,r,z]+·{{-1,0},{0,-1}}.T[t,r,z]+[t,r,z]
1
r
∇
{r,z}
∇
{r,z}
∇
{r,z}
(1,0,0)
T
In[]:=
op==0
Out[]=
-,0.T[t,r,z]+·{{-1,0},{0,-1}}.T[t,r,z]+[t,r,z]0
1
r
∇
{r,z}
∇
{r,z}
∇
{r,z}
(1,0,0)
T
境界条件
境界条件
初期条件と境界条件は,PDEモデリングにおいて重要な要素です.なぜなら,これらは進行しているプロセスの基礎となる物理情報を含んでいるからです.境界条件は,境界における解の挙動を指定します.このシナリオでは,重要な境界条件は卵の外側であり,これは周囲の沸騰したお湯の温度にさらされています.初期条件は卵の開始時の温度,つまり卵が冷蔵されていたのか室温で保管されていたのかによって決まります.
卵が室温にあると仮定し,初期条件として時間0で20°Cを使用し,それをSI基本単位であるケルビンに変換します.
初期温度をケルビン単位で定義します:
卵の内部の温度に対する初期条件を設定します:
境界条件と初期条件は必ず整合していなければならないことを覚えておくことが重要です.領域の外側に100°Cという温度を単に設定して,沸騰したお湯との接触を表現したくなるかもしれません.この方法は全く間違いではありませんが,問題があります.初期条件によると,領域全体は外部境界も含めて初期温度が20°Cに設定されています.そのため,外部境界を100°Cに設定すると,それと矛盾してしまいます.
このシナリオでそのようなことを避けるために,初期温度から始まり,すぐに100°Cの沸騰温度に到達する滑らかなステップ関数を使用します.これは実質的に,卵を沸騰したお湯に入れるという,独自の過程をモデル化していることになります.
具体的には,境界条件を,1/10秒より大きい時間では1となる区分的関数を用いて定義します.より短い時間の場合は,コサイン関数に従って徐々に増加します.この関数をプロットすると,この挙動がより明確になります.これは,卵を沸騰したお湯に入れた直後の最初の1/10秒以内に,境界温度が急速かつ滑らかに上昇する様子をモデル化しています.
境界温度をケルビン単位で指定してください:
境界条件の時間的挙動を指定します:
境界条件の時間的挙動を可視化します:
しかし,今後のこのモデルのバージョンで役立つことを,今回は別の方法で行います.
点要素マーカーを視覚化します:
メッシュの可視化では,点要素および要素マーカーを表示します(詳細については,「Element Mesh Visualization」チュートリアルをご参照ください).この機能を活用して境界条件を定義します.
関数HeatTemperatureConditionを使用する利点の一つは,定義したパラメータを考慮することができる点です.
境界条件をケルビン単位で指定します:
対称性の条件を指定します:
これはNeumann零値に評価されることにご注意ください.Neumann値について復習が必要な場合は,ドキュメント内の次の情報源をご確認ください.
Neumannのゼロ値とは,簡単に言いますと,対称軸に垂直な方向での温度の導関数がゼロであることを示します.また,Neumannのゼロ値は,特に指定がない場合のデフォルト値です.したがいまして,次のこのモデルのバージョンではこれを省略いたします.しかし,このことを念頭に置いておくと有用です.
最後に,シミュレーション時間をどのくらいにするかを定義します.これは最初のシンプルなバージョンであり,PDEのデフォルトのパラメータを使用しているため,まず1秒間と定義しますが,後ほどこの値を更新します.
シミュレーション時間を定義します:
解を得るために,NDSolveValue を使用します.
解を計算し,変数に保存します:
解答を確認します:
この方法で,解を表す補間関数を得ることができ,後でこれをプロットして解を可視化できます.
"ValuesOnGrid" を呼び出すと,各メッシュ座標での関数値を取得できます.その後,解から最小値と最大値を抽出できます.
これをより解釈しやすくするために,温度を摂氏度に変換しますが,モデルに入力するパラメータはケルビンで定義する必要があります.そのため,摂氏度とケルビンの間をより簡単に変換できるよう,オフセットを定義します.
オフセットを定義します:
このプロットを見ると,領域全体がシミュレーション時間の半分(0.5秒)で温度100°Cになるように見えます.このことを確認するために,時間とともに温度がどのように変化するかのアニメーションを作成しましょう.
後で時間を節約するために,先にプロット用の補助関数を作成しておきます.
さまざまな時刻における温度の等高線プロットを作成する補助関数を作成します:
変数timeの関数を作成しました.この関数は特定の時刻に対するプロットを与えます.その後,この関数を時刻のリストに対してマップすることができます.
ヘルパー関数のDensityPlotバージョンを定義します.
いくつかの等高線プロットフレームを作成します:
ノートブックをより軽量にするためにフレームをラスタライズします(これは任意です):
フレームをアニメートします:
このアニメーションは,全体の領域が初期温度から始まり,ほとんど瞬時に100°C(水の温度)まで加熱されることを示していますが,これはあまり理にかなっていません.しかし,これが最終的な解決策になるとは予想していませんでした.現時点では,シミュレーションにエラーや警告がないことが分かっています.このステップは,ワークフローを開始するためのものであり,モデルのこの段階で本当に重要なのはその点です.
バージョン2:問題が発生する可能性がある例です
バージョン2:問題が発生する可能性がある例です
ワークフローが確立されましたので,次のバージョンでは,HeatTransferPDEComponent(PDEを生成するために使用された関数)がデフォルトでは提供するものではなく,より現実的なパラメータを定義することから始めます.
熱方程式には,モデル化する材料の種類によって異なる3つの物理量,すなわち密度,比熱容量,および熱伝導率が含まれます.
卵白と卵黄の両方についての平均値の推定値を用いて,材料パラメータを設定します.
パラメータを確認します:
現実的なシミュレーション終了時間として10分を設定します:
新しいパラメーターで方程式を再生成します:
先ほどと同様に,PDE を解きます.
ここまで,順調です!
解における最小値と最大値を摂氏度で確認します.
これは明らかに間違っています.ここでの氷点下の温度は意味を成しません.氷点下の温度に達した特定の時刻で解をプロットしてみます.
解データ内で最小値の位置を見つけます.
解に対して"Coordinates"を呼び出し,最初の部分を取得することで,解の間に取られた時刻ステップを得ます.
最小値が保存されている時刻ステップを見つけます:
こちらのbadPosition変数({{18, 386}})において,ペアの最初の数値は氷点下の値が発生する時刻のインデックスを表しています.問題となる時刻を抽出するには,時刻ステップの18番目の位置を取得します.
それでは,その正確な時刻で解を可視化できます.
極値を持つ時刻ステップでの解を可視化します:
このシミュレーションの初期段階では,境界付近で多くの振動が発生していることが明確に分かります.そして,そのために–2°C付近の極端な温度が観測されているのです.3Dプロットには,過剰な増加や減少が表示されます.
解のプロットにズームインし,基礎となる有限要素メッシュを可視化します.
2次元メッシュを3次元として可視化する方法についてさらに詳しく知るために,ElementMeshProjection のドキュメントをご覧になることをおすすめします.
境界から数要素以内でアンダーシュートが発生していることが,これで明らかになりました.しかし,これらのオーバーシュートやアンダーシュートは何が原因なのでしょうか.
メッシュの最小エッジ長(メートル単位)を取得します:
要素間の距離は約1 mmですので,境界付近の最初の要素の辺の長さも約1 mmとなります.一方,初期温度が20°Cで境界温度が100°Cであることを考えますと,境界と領域内の他の部分との温度変化は約80°Cになります.
その大きさの要素で境界における約80°Cの急激な変化を解決しようとするのは十分ではありません.それが,熱の逃げ場がない場所で初期温度を大きく下回る温度が得られている理由である可能性が高いです.したがって,メッシュを細かくする必要があります.
メッシュ全体を細かくすれば,おそらく問題は解決しますが,卵の中心部分での細分化は必要ありません.今メッシュ全体を細かくすると,シミュレーションが遅くなりますので,境界(卵殻)付近のみを細分化するのが最適です.
卵全体の領域の符号付き距離関数を作成します:
符号付き距離関数を可視化し,境界に近づくにつれてその値が線形に減少する様子を確認できます.この結果は理にかなっています.というのも,シェル上での距離は予測どおりゼロになるためです.
メッシュ化されたジオメトリ上でrefinement functionを可視化します.
メッシュ細分化関数を作成します:
次に,ToElementMesh の中で MeshRefinementFunction オプションを,refinementFunction に設定します.
メッシュを再定義します:
PDE を解き,進行状況を監視します:
新たに見つかった解の温度範囲を調べます:
素晴らしいです!これで,期待される温度範囲を得ることができました.もちろん,この方法には数値的な誤差が本質的に含まれているため,範囲が正確に{20, 100}にはなりません.しかし,これは有望です.
通常通り,定義した補助関数を新しい解で呼び出し,フレームをラスタライズしてからアニメーション化することで,アニメーションを作成できます.
解答のためのフレームを作成します:
フレームをラスタライズします:
フレームをアニメーション表示します:
オーバーシュートおよびアンダーシュートの問題が解決されたことが分かります.これは,解が妥当であることを確認する重要性と,必要に応じてメッシュを修正できることを示しています.
この時点で,卵を茹でる際に温度がどのように広がるか,おおよそのイメージをつかむことができます.6分間加熱した後の卵の中心の温度を計算してみます.
結果は41.6°Cです.この値は,実際の卵黄の温度を正確に表しているわけではないと思われます.なぜなら,現在多くの仮定を置いているためです.しかし,現時点では良い参考値になります.
このバージョンでは,卵黄と卵白という2つの異なる領域の存在を無視した基本的なモデルを使用しました.次のバージョンでは,その点について対応いたします.
バージョン3:複数の材料領域を表現します
バージョン3:複数の材料領域を表現します
複数素材のジオメトリを作成します:
ここでは,後で二つの領域のパラメータを区別するために,物質領域のマーカーを定義します.これにより,コードがより明確になりますので,ご覧いただけます.
材料領域マーカーを指定します:
"RegionMarker" オプションを使用して,各部分領域に異なるマーカーを割り当てることで,2つの領域を区別します.このオプションはリストのリストとして指定し,各内側のリストには部分領域内の点と,その点に対応する領域マーカーが含まれます.
材料マーカー付きのメッシュを作成します:
これを異なる色で視覚化するには,"Wireframe" オプション内で "MeshElementStyle" を指定します.
複数素材のメッシュを可視化します:
点要素のマーカーを可視化します:
領域を変更したため,境界の要素マーカーも変更されました.そのため,マーカー2と5に対する境界条件を再定義できます.
これらは境界のための要素マーカーであることにご注意ください.これは先ほど定義した領域マーカーとは異なりますので,混同しないようにしてください.
境界条件を再生成します:
点要素マーカーを可視化した後,以前と同様にメッシュ細分化関数を指定してメッシュを再定義します.
マテリアルマーカー付きのメッシュおよびメッシュ細分化を作成します:
洗練されたマルチマテリアルメッシュを可視化します:
素晴らしい結果です!しかし問題を解決する前に,各材料のパラメータを再定義いたします.これは,「Density, Heat Capacity and Thermal Conductivity of Liquid Egg Products」の研究結果に基づいて,卵黄および卵白の質量密度,熱伝導率,比熱容量の平均値の推定値を用いるためです.
それでは,先ほどと同様にPDEを再生成し,解くことができます.
偏微分方程式を再生成します:
このPDEを解きます:
新しい溶液の温度範囲を確認します:
得られた温度範囲は妥当です.解を可視化しましょう.
解答のフレームを作成します:
フレームをラスタライズします:
フレームをアニメーションします:
ここでは,前回のバージョンで単一の領域のみを使用した場合と同様に,卵が均一に加熱されている様子が見られます.しかし,このバージョンは卵の実際の構造をより正確にモデル化しているため,より高い精度を実現しています.
この解答において,38°Cという温度の結果は,与えられた時間に対してやはり低すぎるように思われます.実際,この温度は前回の約42°Cという結果よりもさらに低いです.
一定時間経過後に卵のどの部分が調理されたかを知る方法が必要です.
バージョン4:モデルをさらに洗練します
バージョン4:モデルをさらに洗練します
このバージョンでは,モデルに関与する物理量に使用するデータをより精密にします.また,結果を卵黄および卵白の変性温度(すなわち,卵のタンパク質が変性し固まり始める温度)と比較することで,卵が調理されたかどうかも予測します.
物理量(密度,比熱容量および熱伝導率)をより正確にするために,これら二つの研究によって提供されたデータに依拠します.
前のバージョンでは,卵黄と卵白の両方について物理量の良い推定値を持っていました.しかし実際には,これらの量は卵が加熱されるにつれて時間とともに変化する場合があります.より良いモデルは,これらの量が温度によってどのように変化するかを考慮するものです.
先ほど述べた研究からのデータを用いて,温度と対応する物理量(SI基本単位で表される)とのペアを設定します.目標は,それぞれの物理量(密度,熱伝導率,比熱容量)を温度の関数としてモデル化することです.
まずは密度から始めます.卵黄と卵白の両方について密度データの定義があることにご注目ください.
まずは密度から始めます.卵黄と卵白の両方について密度データの定義があることにご注目ください.
卵黄および卵白の質量密度の測定データを設定します:
"ExtrapolationHandler" オプションは,利用可能なデータ範囲外にある点を処理するために使用します.私たちのデータは完璧ではなく,すべての温度範囲を網羅していません.密度のような特性はほぼ直線的に変化しますので,既存のデータから外挿することで,かなり正確な近似値を得ることができます.
測定したデータおよび補間した関数を可視化します:
温度が上昇すると密度が減少する様子をご覧いただけます.335ケルビンを超える温度データはありませんが,このデータから妥当な外挿が得られております.
次に,同じ手順を伝導率について繰り返します.
卵黄および卵白の熱伝導率の測定データを準備します:
熱伝導率データのInterpolatingFunctionを作成します.
測定されたデータと補間関数を可視化します:
導電率が温度とともにわずかに増加する様子が確認でき,さらにデータの外挿も妥当であることがわかります.
比熱容量については,少し異なるアプローチを取ります.データ点が少なく,他の2つの量と比べてデータの質もあまり良くありません.最適な方法は,データに線形フィットを行うことであり,これによってデータに最もよく合う線形関数を得ることができます.
卵黄と卵白の熱容量の測定データを準備します.
この関数の動作をより詳しく調べるには,LinearModelFit のドキュメントを参照してください.
測定データとそれに適合する線形関数を可視化します:
このプロットは,データの妥当な近似を示し,高温域で良好な外挿も行っております.
物理量の関数を既に用意しましたので,あとは以前と同様にそれらを区分的な関数として指定するだけです.
偏微分方程式を再生成します.
重要な考慮点として,質量密度・熱伝導率・比熱容量は現在,温度Tの関数となっております.同時に,Tは解を求めるための従属変数です.この相互依存関係により,方程式の係数が解そのものによって変化し,PDEが非線形となります.非線形モデルは通常,解くためにより多くの時間と計算労力を必要とします.そのため,初期設定の段階では未精細なメッシュから始めることがしばしば良い考えとなります.合理的な解が得られたら,より高精度を目指して精細なメッシュへ切り替えることができます.しかしながら,今回の場合は最初から精細なメッシュを使用いたします.
これからPDEを解き,進捗状況を監視し,計算にかかる時間とメモリを測定します(通常のノートパソコンで約3分かかります).
この偏微分方程式(PDE)を解きます:
この計算時間の増加は,ほとんどの非線形モデルにおいて予想されます.
新たに発見された解の温度範囲を確認します:
解を可視化する前に,卵黄と卵白の変性温度を定義いたします.これによって,ある時点で卵の二つの領域が調理されているかどうかを把握しやすくなります.これらの温度は,先に引用したJournal of Food Measurement and Characterizationの研究から得られたものです.
卵白の変性温度を設定します:
卵黄の変性温度を設定します:
溶液のアニメーションとともに,変性温度を可視化することができます.そのための最適な方法は,領域内の変性温度の点を示す等高線を含む ContourPlot を用いることです.
私は密度プロットを好みますので,ここでは新しい関数TemperatureDenatureDensityPlotを定義します.この関数は,以前定義した補助関数TemperatureDensityPlotを呼び出します.以前と同様にShowを使ってプロットを表示しますが,今回は卵黄と卵白の変性温度の等高線も表示されます.等高線は追加のContourPlotでプロットされます.
温度分布の密度プロットを作成し,変性温度を可視化するための補助関数を作成します:
解答のためのフレームを作成します.
フレームをラスタライズします:
フレームをアニメーションします:
卵白の変性温度は緑色の破線で示されており,卵黄の変性温度はマゼンタ色の破線で示されています.このプロットは,予想される通り,卵の内部での熱の拡散を示しています.
6分経過時点で,卵白全体が変性温度を超えていることが分かります.これは,モデルにとって有望な結果です.
しかし,10分の時点を見ますと,全体の領域が卵白の変性温度を超えていますが,卵黄の一部は卵黄の変性温度を超えていませんので,10分経過しても卵黄は十分に加熱されていないことが分かります.
ゆで卵は通常10〜12分間加熱しますので,10分加熱すると卵黄が完全に固まるはずです.そのため,私たちのモデルはまだ完全ではありません.なぜなら,このモデルでは10分後でも卵黄が完全に調理されていないことを示しているからです.
卵黄の変性温度を調べます:
モデルのバージョン2および3では,それぞれ38°Cと42°Cの温度を得ました.ここバージョン4では,6分経過時点で卵の中心部の温度がおよそ40°Cとなります.黄身の変性温度と比べると,これでもまだ低すぎるように思われます.次のバージョンでこの点を改善する必要があります.
バージョン5:リアルな卵の形状を作成します
バージョン5:リアルな卵の形状を作成します
モデルが期待した答えを出さない場合,さらに洗練させることが良い方法となる場合があります.そのための一つの方法は,自分たちが立てた仮定を見直すことです.重要な仮定の一つは,卵が円形であるというものですが,実際にはこれは正しくありません.このバージョンでは,卵の形状をより現実的に再現するジオメトリを作成します.
卵黄の形状については,依然として円でよく近似できると考えます.しかし,卵殻の形状については,より現実的な方法で近似することができます.
卵形状の座標を計算する補助関数を作成します:
パラメータを操作すると,さまざまな形状を表示できます.
半径の値を確認します:
高さ5cm,幅4cmの卵に対して新しいeggShapePoints関数を呼び出します.
パラメータを指定することによって卵の形状の座標を作成します:
スプラインを作成します:
次に,RegionUnion を使い,半円と卵殻曲線を結合して,領域の骨組みを作成できます.これは,以前と同様ですが,今回は卵殻のためにスプライン曲線を使用して行います.
サブリージョンを持つ卵形ジオメトリを作成します:
次に,先ほどと同様に領域マーカー付きのメッシュを作成します.
マテリアルマーカー付きのメッシュを作成します:
点要素マーカーを視覚化します:
ご覧のとおり,卵の外側のマーカーは依然として2と5ですので,境界条件を再定義する必要はありません.
境界条件を確認します:
refinementFunction をご確認ください.
すでに良い段階に来ていますが,少し問題があります.前回のバージョンでは,領域が円形であることに基づいた精緻化関数を作成しました.新しい領域のために,新しい精緻化関数を作成する必要があります.
要素は,その半径方向の距離に基づいて,殻に対して洗練化することができます.この殻は画像中の緑色の線で示されています.赤い点は要素を表し,その点は対称軸(青色の線)までの半径方向の距離と,殻までの半径方向の距離を持っています.
緑色で示された距離に基づいた細分化関数を定義したいと思います.殻までの距離が大きいほど細分化は低くなり,殻までの距離が小さいほど細分化は高くなります.さらに,緑色の線の距離を得るためには,対称軸から卵殻までの距離(マゼンタで示されています)から青色の線の距離を引きます.
まず,卵殻曲線において,与えられた z に対する円筒座標の値 r を返す関数を定義します.
殻の座標を補間します:
単純にシェル座標の点の順序を反転させていま す—{r , z} を {z , r} に置き換え,その後で補間を行っています.ここで,rShell は z を引数として受け取り,r を返す補間関数です.
rShell 関数を与える関数をプロットします:
この関数のグラフは予想どおりのものです.卵殻の座標rは,zの関数として表されます.
それから,残りに必要な作業は,各点ごとに軸からシェルまでの距離から座標r の絶対値を引いて計算することです.これにより,各点からシェルまでの距離(緑色で示されています)を求めることができます.
シェルまでの動径距離を求める関数を定義し,それをプロットします.
このプロットは,各点からシェルまでの距離が,予想通りほぼ直線的に減少することを示してます.
シェルまでの距離の三乗をプロットします.
距離の三乗をプロットすると,シェルまでの距離が小さくなるにつれて,その値はより速く減少します.
精緻化関数の正しい動作を見つけるには,試行錯誤のプロセスになることを念頭に置いてください.
細分化関数を定義します:
ここでは,refinementFunction を,要素の面積がシェルまでの距離の三乗より大きい場合にその要素が細分化されるように定義しました.また,過度な細分化を防ぐためにオフセットも加えています.言い換えますと,要素の大きさはシェルまでの距離の三乗に比例します.シェルに近づくほど,要素は小さくなります.
複数の素材のメッシュを可視化します:
ご覧のとおり,卵の中心付近の要素は全く細分化されていませんが,卵殻の近くでは非常に細かく細分化されていることがわかります.
これまでと同様にPDEを解くことができますが,その際にかかった時間と消費されたメモリを測定します.このメッシュを用いた非線形モデルの計算は,終了までに約6分かかります(ちょうど私が好きな卵料理の調理時間と同じです!).モデルを作成する際には,計算時間に注意してください.このため,Monitor が役立つのであり,最初は非細分化メッシュを使用することが良いアイデアなのです.
PDE を解き,進行状況を監視し,計算時間とメモリ使用量を測定します.
新たに見つかった解の温度範囲を調べます:
今回も妥当な範囲の温度が得られましたので,アニメーションを続けます.
解答のためのフレームを作成します:
フレームをラスター化します:
フレームをアニメーションします:
重要な点が二つあります.まず第一に,6分では卵黄がまだ固まっていないことが予想通りであるという点です.次に,10分では卵白と卵黄の両方の変性温度を全体が超えているため,この時点で卵が完全に調理されたと考えることができます.これでモデルは卵を調理する現実的な時間枠を反映しています.素晴らしいですね!
6分後には,卵の中心温度がほぼ60°Cという,より現実的な温度になっております.これは,以前のモデルでは中心が約40°Cまでしか上昇しなかったことと比較できます.このことから,加熱のダイナミクスにおいて幾何学的形状が重要な役割を果たしていると結論付けることができます.
バージョン6:Modelを使用します
バージョン6:Modelを使用します
モデルが満足のいく結果を出すようになりましたので,このモデルを使って予測を行うことができます.
以前は,卵を室温から沸騰したお湯に入れる状況を紹介しました.ここでは,より現実的なケースとして,冷蔵庫から直接取り出した卵を考えます.この場合,卵の初期温度は8°Cであると仮定します.これをモデル化するためには,初期条件と境界条件を修正する必要があります.特に,境界条件は新しい初期温度から始まり,素早く100°Cに上昇する必要があります.
まず,初期温度を8°Cに設定し,それをケルビンに変換します.
卵の初期温度を設定します:
次に,新しい初期条件を定義します.
卵の内部の温度に対する初期条件を設定します:
また,新しい初期温度から始めるため,境界条件関数も修正する必要があります.
境界条件の時間的な挙動を指定します.
境界条件の時間的挙動を可視化します:
境界条件を再生成します:
次に,前と同様にPDEを解きます.ここでは新しい変数solutionFridgeに解を保存しますので,以前の解と比較することができます.
PDEを解き,その進行状況を監視し,計算にかかる時間とメモリを測定します.
新たに見つかった解の温度範囲を確認します:
冷蔵された卵に対して妥当な範囲が得られます.
以前と同じ方法で解をプロットすることができます.
解答のためのフレームを作成します:
フレームを Rasterize します:
フレームをアニメーションします:
6分では,卵黄はまだ加熱されていません.これは前回のバージョンと同様です.10分では,卵全体がしっかりと加熱されているように見えます.アニメーションを見るだけでは,冷蔵された卵と常温の卵の間に意味のある違いがあるかどうかを判断するのは難しいです.
良い選択肢の一つは,黄身の中心から卵殻を通る直線上の温度を,冷蔵された卵と室温の卵の両方についてプロットすることです.より単純なプロットによって,密度プロットでは見えにくい微妙な側面が明らかになるかもしれません.
この直線上にあるrの値について,温度をプロットしています.右側の境界におけるrの値が必要です.先ほど定義した関数rShellを使うことで,座標rの値をzの値が0.005,つまり卵黄の中心で取得できます.
この直線に対する r の値を取得します:
これで,注目している直線の右端のrの値が分かりました.次に,その直線について,特にrが0.から0.0179までの値の場合に,温度がどのように上昇するかを知りたいです.例えば6分後に,これらのrの値に対して,簡単なプロットを作成できます.
卵黄の中心から卵の外側までの放射線に沿って,6分後の解をプロットします.
6分後の温度は,r のすべての値において,冷蔵庫から直接取り出した卵の方が,常温の卵よりも低いことが分かります.特に,卵の中心部では約5°Cの差があることが分かります.
卵が調理されているかどうかを予測したいので,変性温度を可視化する方法も必要です.その一つの方法は,卵が変性温度に達するrの値を示す垂直線をプロットすることです.
溶液が変性温度と等しくなるr の値を取得します:
冷蔵庫から取り出した卵の黄身の変性温度は,中心から約10 mm(1 cm)の半径で到達することが分かります.これをプロット内の線で表現することができます.
卵黄の中心から卵の外側まで,半径方向の直線に沿った解をプロットします.
卵黄に加えて,室温の卵と冷蔵された卵の卵白の変性温度も必要です.これを手動で行う代わりに,補助関数を定義することができます.
この関数は,解と関心のある変性温度の値を受け取り,FindRoot を使用して,変性温度が発生するプロット上の線に沿った半径を決定します.最後に,プロット上に表示できる垂直線を返します.
溶液の変性温度の線を作成するための補助関数を作成します:
次に,実際のプロット用の関数を作成します.この関数では,Show を使用してプロットと変性温度の線を表示し,denatureLine 補助関数を呼び出すとともに,卵黄の半径を示す線も表示します.
さまざまな時刻における2つの解の変性位置をプロットする補助関数を作成します:
6分経過した時点でプロットを見ると,二つの点に気付きます.まず最初に,卵白の変性がすでに卵の中心部まで進んでいることが明らかです.次に,室温の卵の場合,卵黄の変性温度が冷蔵卵と比べて中心にずっと近づいていることが分かります.
これはどういう意味でしょうか.室温の卵の場合,黄身はちょうど加熱され始めています.しかし,冷蔵した卵の場合,黄身はまだ半熟のままであり,変性温度は黄身の端にかろうじて達しただけです.
したがって,とろとろの黄身にしたい場合,2つの選択肢があります.
1
.冷蔵された卵を使用し,ちょうど6分で沸騰したお湯から取り出します
2
.室温の卵を使用し,水から数秒早めに取り出します
同じプロットを作成しますが,今回は8分間の場合について見てみますと,興味深いことが分かります.約8分経過した時点で,室温の卵全体が黄身の変性温度を超えており,これは卵白の変性温度もすでに超えていることを意味します.しかし,冷蔵庫で冷やした卵の黄身の一部には,まだその温度に達していない部分があります.
これはどういう意味ですか?もし完全に火を通した卵黄をご希望で,冷蔵庫から出したばかりの卵を使う場合は,8分を過ぎてから数秒間沸騰したお湯に入れておいてください.常温の卵を使う場合,そしてお腹が空いていて待ちたくない場合は,8分でちょうど良いかもしれません.
結局のところ,要点はこれです
結局のところ,要点はこれです
まず,簡単な最初のバージョンを作成して,ワークフローを始動させました.いきなり複雑なモデルのコーディングに取り掛かると,間違いを犯す可能性が高くなります.さらに悪いことに,そのような間違いは特定するのが難しい場合があります.もし結果が正しく見えない場合,モデル自体に問題があるのか,それともコードに単純なミスがあるのか,必ずしも明確ではありません.
シンプルに始めて徐々に複雑さを加えていくことで,一歩一歩バグを導入する可能性を減らします.また,モデルがどのように振る舞うかについても,より良く理解できます.
先ほどご説明した通り,卵の外側と内側の温度が急激に変化することで,初期のモデルの一つで数値的な問題が発生していました.この問題を解決するために,境界部分でリファインメント関数を用いてメッシュを細かくしました.この方法は,全体のメッシュを細かくするよりもはるかに優れており,計算速度が大幅に遅くなることを防ぎます.
最後に,実験データを2つの方法で使用しました.まず,実験データを用いて熱方程式のパラメータを設定することによりモデルを構築しました.次に,データを用いて実際の卵の加熱時間と我々の結果を比較しました.
いくつかのことを結論付けることができます:
◼
まず簡単なバージョンを作成することは良い考えです.
◼
モデルの複雑さを段階的に高めていくことで,物事をより簡単にすることができます.
◼
解の変数において急激な変化や不連続性がある場合,数値的不安定や精度の低下が生じることがあります.そのような領域では,局所的なメッシュの細分化や平滑化を行ってください.
◼
実際のデータと結果を比較できる方法があると便利です.