controlDictにて熱流束を与えたが温度が変化しない (Openfoam ESI版 v2412 )

35 views

Skip to first unread message

fujimoto

unread,

Aug 14, 2025, 11:25:51 PMAug 14

to OpenFOAM

お世話になっております。

fujimotoです。

以前、質問させていただいた内容( fvOptionsの使用でおきたエラー Entry 'name' not found in dictionary )（https://groups.google.com/g/openfoam/c/_jk5Olh4NXI）に関連するものです。

添付しているよう写真のように、上側に気体(topAir)、下側に固体(bottomSolid)を想定したモデルで、その境界(bottomSolid_to_topAir)に面する固体側のいくつかのセルに熱量を与えたいと考えています。(ガウス分布のレーザを想定している）

そのために、controlDictに上記の内容を実行するコードを書きました。(添付しているケースディレクトリ内にあります）

しかし，実行したところ、温度が初期値のまま一切変化しない状態です。

controlDictの内容・流れを簡単に記しておきます。

内容

熱量を与える時間は0 ~ 0.05秒 (計算するのは0.1秒まで）

ビーム半径の1.5倍の領域まで熱量を与えることを想定

流れ

1. if文でシミュレーション内での経過時間がレーザを与える時間か確認

↓

2. もしそうであれば、境界bottomSolid_to_topAirに面するセルを列挙

↓

3. そのセルのうちビーム半径の1.5倍の領域内のセルを列挙し、熱量を与える。

上の1, 2, 3の手順それぞれにおいて、予想した挙動（望んだ時間にif文の中に入る，意図したセルが列挙されている）であることは確認しています。

お手数おかけしますが、なぜ全く温度変化がないのか、ご教授いただければ幸いです。

（添付しているケースディレクトリはchtMultiRegionFoamのコマンドで実行できる状態となっています）

fujimoto

--------------controlDict----------------------------

customHeat

{

type coded;

libs ("libutilityFunctionObjects.so");

name customHeat;

region bottomSolid;

writeControl timeStep;

writeInterval 1;

codeExecute

const Time& myTime = mesh().time();

const scalar myDeltaT = myTime.deltaT().value(); //時間step⊿tを取得

volScalarField& h = const_cast<volScalarField&>(mesh().lookupObject<volScalarField>("h"));

const vectorField& C = mesh().C();

//const scalarField& V = mesh().V();

const scalar p0 = 1000; //レーザ出力

const scalar absoptionRate = 0.4; //吸収率

const scalar x0 = 0.0; //光軸中心座標

const scalar y0 = 0.0;

const scalar radius = 0.01; //ビーム半径

const scalar endTime = 0.05; //レーザ照射を終了する時刻

const scalar pi = 3.1415; //円周率

scalar p = 0;

scalar r = 0;

scalar x = 0;

scalar y = 0;

OFstream logFile("output.txt",

IOstream::ASCII,

IOstream::UNCOMPRESSED,

IOstreamOption::APPEND); //IOstreamOption::APPEND　IOstreamOption::NO_APPEND

word patchName = "bottomSolid_to_topAir";

logFile << "\nTime: " << myTime.value() << endl;

if (myTime.value() < endTime)

{

label patchI = mesh().boundaryMesh().findPatchID(patchName); //境界"bottomSolid_to_topAir"に属するパッチid

// パッチの開始face番号とサイズ

const polyPatch& pp = mesh().boundaryMesh()[patchI];

label startFace = pp.start();

label nFaces = pp.size();

DynamicList<label> cellsOnPatch;

label cellI;

DynamicList<scalar> areaOnPatch;

scalar area;

for (label faceI = startFace; faceI < startFace + nFaces; ++faceI)

{

cellI = mesh().faceOwner()[faceI]; // 該当するパッチを保有するセル

cellsOnPatch.append(cellI); //セルのリスト

area = mag(mesh().Sf()[faceI]); //該当するパッチの面積

areaOnPatch.append(area); //パッチの面積のリスト

}

logFile << "=======================cellsOnPatch======================" << endl;

logFile << cellsOnPatch << endl;

logFile << "=======================q_location======================" << endl;

forAll(cellsOnPatch, celli)

{

x = C[cellsOnPatch[celli]].x(); //該当するセルの中心のx座標

y = C[cellsOnPatch[celli]].y();

r = sqrt(pow(x-x0, 2) + pow(y-y0, 2)); //光軸中心からの半径

if (r < 1.5*radius)

{

logFile << C[cellsOnPatch[celli]].x() << " " << C[cellsOnPatch[celli]].y() << " " << C[cellsOnPatch[celli]].z() << endl;

p = 2 * absoptionRate * p0 * exp(-2 * pow(r / radius, 2)) / (pi * pow(radius, 2)); //該当する位置のパワー密度(ガウス分布)

h[cellsOnPatch[celli]] += p * areaOnPatch[celli] * myDeltaT; // パワー密度 x 面積 x 時間step で熱量を算出し加える

}

#};

}

cu_melting_laser.zip

画像1.png

画像2.png

fujimoto

unread,

Aug 17, 2025, 8:19:25 PMAug 17

to OpenFOAM

申し訳ありません、chtMultiRegionFoamを実行するだけで良いと書いてますが、Allrun-serialを実行してもらわないと動きません。