レーザーのシミュレーション設定について

[amadai]

お世話になります．OpenFOAM初心者のamadaiといいます．penguistisさんのサイトや本でOpenFOAMを勉強中です．今回レーザーのシミュレーションについて質問があり投稿致します．

チュートリアル(tutorials/multiphase/icoReactingMultiPhaseInterFoam/solidMelting2D)を参考に金属の相変化を伴うレーザーの照射シミュレーションをしてみたところ，計算結果が不自然なものになりました．

OpenFOAM-v2006のケースファイルを添付いたしますので，どのような情報でも結構ですのでアドバイスを頂けないでしょうか．

アルミの金属板を空気で挟み(図1)，上からレーザーを照射する計算設定です(図2)．

金属の熱伝導率は空気と比べて大きいので，金属が先に熱くなると思うのですが，長時間計算した温度分布を見ると空気のほうが熱くなっており，金属温度は融点に達する温度よりもはるかに低いままでした(図3)．

また，金属部分の温度を可視化してみると，金属上部とともに金属下部も熱くなっていました(図4)．

金属上部から融点に達して釣鐘状に溶けていくイメージだったのですが，想像と全く違う結果でした．

設定ファイルが不正で，空気にレーザーの熱が過剰に伝わっているのでしょうか？

OpenFOAMの経験が浅いため色々と馬鹿な計算設定をしていると思われますが，もし可能でしたら，レーザー照射のような計算の正しい(もしくは不正に環境へ熱が伝わらない)設定を教えて頂けますと幸いです．

どうぞよろしくお願いいたします．

[M. Tamura]

amadai様

田村と申します。

解決済みかもしれませんが、参考になれば幸いです。

問題は、光が通過した領域で気相・固相に関係なく発熱が生じていることだと思われます。

このことは発熱量Q（amadai様の図2）から理解できます。

解決法として、radiationProperties内で、

alpha.solidとalpha.liquidのみに吸収係数などの物性値を設定する必要があります。

absorptionEmissionModelを変更する必要があり、例えば、

absorptionEmissionModel localDensityAbsorptionEmission;  //constantAbsorptionEmission;

localDensityAbsorptionEmissionCoeffs{

    alphaNames (alpha.solid alpha.liquid);

    aCoeff     (10 10);

    eCoeff     (0 0);

    ECoeff     (0 0);

}

のような書式です。

radiationPropertiesに元々あったabsorptivityなどはコメントアウトして問題ないと思います。

ちなみに、focalLaserRadiusは、

DTRM法においてrayを極座標で配置する際の最大半径を与えるものと認識しています。

（Gaussianでなくuniformだと、その強度を計算するためにも利用されるみたいですが）

Gaussianのスポット半径sigmaに対し、2倍程度の値がないと、

Gaussianのプロファイルを正しく評価できないと思うのですが、大丈夫でしょうか？

このあたりを修正すると、中心から溶ける様子を確認できましたが。

温度場を正しく計算できていないためか、αの振る舞いも怪しく、

正しく計算するためには、計算領域の広さなど、他にも修正が必要かと思います。

2021年1月2日土曜日 19:25:50 UTC+9 amadai:

[amadai]

田村様

いつもお世話になっております．雨大と申します．

メールの確認が遅くなり申し訳ございませんでした．

radiationProperties，およびfocalLaserRadiusの内容について

詳細にご教示下さりありがとうございました．

戴いたアドバイスを元に両ファイルを変更し計算してみたところ，

無事に固相のみが中心から発熱していることが確認できました．

未解決でしたので大変助かりました．

ご助言頂き誠にありがとうございました．

取り急ぎお礼申し上げます．

体積率などの計算結果が不自然なのはおっしゃる通りで，引き続き

調査を続けてみます．固相と液相の吸収係数を変えるだけでも，

劇的に固体の融解挙動が変わってしまいます．

例えば固相と液相の吸収係数をそれぞれ250, 500としたときの

照射初期における固相の体積率が図1で，図2が温度です．

レーザーが金属を貫通してしまうのは不自然なので，本来ならば

金属上部に熱が吸収され，上部から徐々に熱伝導で温度が上がって

いってほしいのですが，なかなか上手くいきません．

また，熱伝達率や熱伝導率の値には具体的なステンレスやアルミの

数値を入れているのですが，constant/phasePropertiesで指定する

固液共存領域のモデリング(Voller Prakashモデル)定数Cuや，

融解モデル(Leeモデル)のモデル定数Cの値によっても，融解挙動が

相当変わってしまいます．

constant/phaseProperties

interfacePorous

(

    (solid and liquid)

    {

        type            VollerPrakash;

        solidPhase      alpha.solid;

        Cu              1e7;　<---------------

    }

);

massTransferModel

(

    (solid to liquid)

    {

        type            Lee;

        C               420;　<---------------

        Tactivate       1400;

    }

);

もし，各相の吸収係数と，2つのモデル定数の選び方についての

ご知見がございましたら，お手すきの際にご教示頂くことは

できないでしょうか．

何卒よろしくお願いいたします．

2021年1月28日木曜日 22:21:58 UTC+9 M. Tamura:

[M. Tamura]

雨大様

田村です。

そもそも、物質内でどのように吸収を計算するかに関してですが、

DTRM法はLambert-Beerの法則などに関連して、光強度の変化を評価するようです。

https://www.afs.enea.it/project/neptunius/docs/fluent/html/th/node114.htm

が参考になり、

https://develop.openfoam.com/Development/openfoam/blob/OpenFOAM-v2006/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.C

のソースコードを見ると、そのことを理解できます。

eやEを0とする場合は話が非常に単純で、要するに、aの吸収係数を持つ物体内で光強度は、

I(x) = I0*exp(-a*x)

で減衰し、減衰した分の出力が発熱となります。（I0は物体入射前の光強度です）

aは吸収係数ですが、上述の式より物質内での光の侵入長（Iが1/eになる長さ）の逆数に対応します。

雨大様が設定された250という値ですと、少なくとも1/250 = 4[mm]は光が物質内を伝搬します。

具体的な吸収係数の値は、物質の複素屈折率の虚部、消衰係数と関連しますが、そのあたりの詳細は調べてみてください。

「吸収係数　消衰係数」や「refractive index Al」などで検索すると良いかと思います。

ただし、複素屈折率も温度で変化するはずですが、その具体的な情報は中々見つからないかと思います。

割り切って常温の値を用いるか、妥当だと思う値を適当に採用するなど、対応は様々だと思います。

また、別の注意点として、今のモデルでは反射を考えていないはずです。

物質の厚みが光の侵入長よりも十分に長ければ、投入したレーザーの出力全てが発熱に変換されます。

つまり吸収率100%となり、これは非現実的です。

反射を設定するモデルreflectionModelなども用意されているようですが、

簡単な対応方法としては、物質の吸収率をあらかじめ計算し、投入する出力に吸収率をかけて減らしておけば、ひとまずは良いかと思います。

フレネルの式から反射率を計算し、吸収率も計算できるかと思いますが、

上述の屈折率の温度依存性の話から、そこまで厳密な値を用いる意味もあまりない気がしています。

なお、深い穴を掘るような計算をする場合は、穴の側面での反射も重要になるかと思いますので、その点も注意が必要でしょう。

Leeモデルの係数に関してですが、論文によっては0.1～10^7と幅広く様々な値が用いられているようです。

単位は[1/s]で、まさに相変化の速度を与えるような係数です。

経験的に定める値のような気がしますが、その逆数ぐらいのオーダーで相変化が進行すると考えて良いのではないでしょうか？

例えば、dm/dt = C*ρ*α*(T-T_a)/T_aで与えられますが、

仮に、T_a=1400 K, T = 2*T_a = 2800 K, α=1, C=10の状況を考えると、dm/dt = 10*ρであり、

0.1秒でそのセル内の物質全ての相変化が進行することに相当するかと思います。実際の検証はしておらず、想像でしかないですが。

C=420であれば、1/C = 2.4[ms] 程度ですが、それを妥当と思われるかどうかではないでしょうか？

また、こちらの値を調整する前に、まずは温度場を正しく計算できるようにパラメータを定める方が重要な気がします。

Voller Prakashに関して、仕組みとしては、

運動量保存の式に対し、-Cu * α_s^2 / ( (1-α_s)^3 + 1e-3 ) * U

のようなソース項を加えているものと認識しています。

α_sは固相の体積分率ですが、要するに固相が存在する領域でUを減衰させる効果を持ちます。

正直なところ、私も適切な値に関して詳しくコメントできません。

（mushy regionの振る舞いについて詳しくないです、すみません）

これもまた妥当だと思われる値を設定していただくしかないかと思います。

ひとまずはデフォルトの値を採用し、温度場と相変化を正しく計算できるようになった後にご検討されるのが良いかと思います。

どうぞよろしくお願いします。

2021年2月6日土曜日 1:39:55 UTC+9 amadai:

[dai ama]

田村様

いつもお世話になっております．雨大です．

ご返信が遅れまして大変申し訳ございません．

浅学寡聞のためご指導を頂き大変ありがたく存じます．

DTRM法の実装について詳細にご教示下さり，誠にありがとうございました．

田村様のおかげで，吸収係数aの物理的な意味について理解させて頂くことができました．

吸収係数は侵入長の逆数と考えれば良いのですね．与え方が悩ましいですが，まずは

固相と液相の両方に常温の吸収係数を与えて計算し直してみます．

引き続きご紹介頂いたサイトやLambert-Beerの法則や吸収係数，消衰係数などを検索し，

現象の理解に努めたいと思います．

反射を考慮していない点についてもご指摘下さり，誠にありがとうございました．

ご教示頂いた方法にしたがい，まずは吸収率をレーザー出力に乗じてみたいと思います．

OpenFOAMでは反射のモデルも実装されているのですね．恥ずかしながらフレネルの式も

不勉強でございますので，こちらも勉強したいと思います．ご教示ありがとうございました．

また，Leeモデルの物理的解釈についても詳細にご教示頂きまして，誠にありがとうございました．

こちらも田村様が明快にご説明下さったお蔭で，大変良く理解させて頂くことができました．

係数Cは相変化速度の逆数なのですね．仮に液相の温度が融点付近にあると仮定しますと，

融解の写真などから係数Cの妥当な値を与えることができるのではないかと思いました．

貴重なご助言を頂きありがとうございました．

まずは温度場を正しく計算すべきとのこと，重要なご指摘を頂戴し誠にありがとうございました．

田村様がおっしゃる通り，現状の私の計算ですと，液相温度が数万度のような非現実的な温度に

上昇してしまい，全体として温度場が正しく計算できておりません．

このような場合の対処方法といたしましては，液相に注入されたレーザーの熱を逃がすために

液相から気相への相変化モデルを導入する，という理解で合っておりますでしょうか？

(単に計算格子や収束条件を厳しくするだけでは解決いたしませんでした．)

また，このような三相の相変化を検証するベンチマークなどはあったりするものなのでしょうか？

もしご知見がございましたら，お手すきの際で結構ですので，ご助言を頂戴できますと幸いです．

ご多忙のところ恐れ入りますが，ご検討のほど何卒よろしくお願いいたします．

また，Voller Prakashモデルに関しましても，式の意味を詳細にご教示頂きありがとうございました．

運動方程式にソース項を加えることで，相の境界速度を調整していると理解いたしました．

まずはご教示頂いた通りに，デフォルトの値を設定し計算してみたいと思います．

ご教示頂きましてありがとうございました．

2021年2月7日(日) 14:26 M. Tamura <simt...@gmail.com>:

--
このメールは Google グループのグループ「OpenFOAM」のトピックを登録しているユーザーに送られています。
このトピックの登録を解除するには https://groups.google.com/d/topic/openfoam/SXoTJN_sS3g/unsubscribe にアクセスしてください。
このグループを退会し、グループのすべてのトピックの登録を解除するには openfoam+u...@googlegroups.com にメールを送信してください。
このディスカッションをウェブ上で閲覧するには https://groups.google.com/d/msgid/openfoam/b7be7b25-faaa-4a8a-8d58-5c1f63b98ff5n%40googlegroups.com にアクセスしてください。

[M. Tamura]

雨大様

田村です。

最初に断っておくことを忘れていましたが、私もこの分野については勉強中でして、

間違って理解している可能性や、詳細まで把握していないことなどありますのであしからず。

訂正として、

＞係数Cは相変化速度の逆数なのですね

LeeモデルのCが速度のようなものと認識しています。

ただし、あくまで私の理解です。

文献などでそのように説明されているのを見たわけではありません。

（そもそも現象論的な式のためか、詳細な説明を見たことが無いのですが…）

蒸発についてチュートリアルがあり、

kineticGasEvaporationを使えば、何となく導入はできるかと思いますが、

係数の調整や金属のgasモデルの与え方に悩むかと思います。

私も挑戦したことがありますが、上手くパラメータを設定できなかったのであきらめました。

実際に再現したい状況に対し蒸発が大きく寄与していると考えてらっしゃるのであれば、蒸発を導入するのは良いと思うのですが、

現在のシミュレーションを成功させる目的としては良くない気がします。

温度が高くなりすぎる原因として、熱放射など他の要因もあるかもしれません。

温度を現実と合わせることは不可能でしょうし、まずは、思いきって出力を1桁程度減らすなどして、

シミュレーションを成功させることを優先すべきでしょう。

それでもシミュレーションが成功しないのであれば、設定内容に異常がある可能性を疑うことができます。

そもそも、設定を眺めていると幾つか気になる点があります。

致命的だと思うのは下記の１で、他は確認・提案程度です。

１．fvSolutionのcAlphasの設定ですが、全ての界面について0でしたが修正されましたか？

　　設定の意味をご存じでない場合は、詳細な説明がinterFoamと関連してwebに多数ありますので、お調べください。

　　ひとまず、liquid/airとsolid/airについては1が妥当でしょう。

２．gasのequationOfStateがincompressiblePerfectGasで、pRefが1e-4でした。

　　真空を想定されているということでしょうか？

　　そうでないのであれば、チュートリアルに倣って1e5が良いのではないでしょうか？

　　そこまで大きな影響はないかもしれませんが。

３．計算格子の条件を厳しくしたとのことですが、メッシュサイズを小さくするよりは、

　　計算領域を広げる方が、意味がある気がします。

４．先に温度を正しく計算できることを確認すべきといった手前恐縮ですが、

　　LeeモデルのC=1は、100ぐらいに上げておいても良い気がします。

　　細かい調整は後回しで良いでしょうが、

　　小さすぎると中々溶けず、温度が高くなりすぎる原因になる気がします。

どうぞよろしくお願いします。

2021年2月10日水曜日 22:31:25 UTC+9 amadai: