後編はこちらです。
https://www.shirokurohitsuji.studio/2020/08/02/adjointoptimisationfoam-usermanual/
adjointOptimisationFoam とは
OpenFOAM でトポロジー最適化に使うソルバーです。
v1906から追加されました。
この記事はユーザマニュアルを読んで今後の自分のために分かる限りで日本語にしたもの。長くなってしまったのでこの記事は「optimisationDict」についてだけ。
でもよくわかっていないところもたくさんあるので解釈が違うとか嘘が書いてあったりするかもしれないです。
見つけたら教えてください。
https://www.openfoam.com/documentation/files/adjointOptimisationFoamManual.pdf
optimisationDict
adjointOptimisation を定義するDict。
「system」フォルダの中に作成する。
OptimisationManager
optimisationManager singleRun ;操作モードの設定を行う。
singleRun は、交互に計算するのではなく primal 方程式と adjoint 方程式を1度ずつだけ計算する手法。
v1906では singleRun 以外の選択肢は用意されていない。
PrimalSolvers
primalSolvers
{
p1
{
active true;
type incompressible;
solver simple;
solutionControls
{
nIters 3000;
averaging
{
average true;
startIter 1000;
}
residualControl
{
"p.*" 1.e-7;
"U.*" 1.e-7;
}
}
fvOptions { }
}
}primal方程式の設定を行う。
- active (デフォルト値:true）この primal ソルバーを有効にするかどうか
- type primalソルバーの解析対象のタイプ。現時点では「incompressible」しか使えない。
- solver ソルバーのアルゴリズム。simple か RASTurbulenceModel を選択する。simple は simpleFoam と同じ働きをする。RASTurbulenceModel は constant/turbulenceProperties に定義された値を用いて乱流場を解く。
- useSolverNameForFields (デフォルト値:false）「true」にすると、変数名にこの primal ソルバーの名前が追記される。
- 例えば今回の例では primal ソルバーの名称が「p1」なので、「U」として定義されていた変数は、「Up1」となる。このオプションを有効にした場合には、fvSolution と fvScheme を手動で対応させなければならない。場合によっては境界条件の設定においても変数名を指定する場合があるのでその部分の調整も必要になる。
- 「multi-point runs」を行うときに有効、「single-point runs」のときには無効のままにすることが望ましい。
- このオプションが有効なときでも「Up1」が見つからなければ「U」から境界条件を読み取る。
- solutionControls resiidualControl など、通常 fvOptions の中の「SIMPLE」などで定義する項目の設定ができる。
- niters primal 方程式を解くときの最大反復回数
- averaging 平均化オプション。指定しなくても良い。simpleFoam のような定常ソルバーを使用して非定常問題を解くときなどに有効にする。
- average （デフォルト値:false） 平均化するかどうか。trueにすると、「U」「p」や乱流場を定義する変数などが平均化される。平均化された変数は「UMean」のように名前の末尾に「Mean」が付く。
- startIter 平均化を開始するイテレーション数。
- fvOptions ソース項の設定に使用する。指定しなくても良い。
adjointManagers
adjointManagers
{
am1
{
primalSolver p1;
operatingPointWeight 1;
adjointSolvers
{
as1
{
// choose adjoint solver
//----------------------
active true;
type incompressible;
solver adjointSimple;
useSolverNameForFields false;
computeSensitivities true;
// manage objectives
//------------------
objectives
{
type incompressible;
objectiveNames
{
losses
{
weight 1;
type PtLosses;
patches (Inlet Outlet);
}
}
}
// ATC treatment
//--------------
ATCModel
{
ATCModel standard;
extraConvection 0;
zeroATCPatchTypes ();
nSmooth 0;
maskType faceCells;
}
// solution control
//------------------
solutionControls
{
nIters 3000;
printMaxMags true;
residualControl
{
"pa.*" 1.e-7;
"Ua.*" 1.e-7;
}
}
}
}
}
}各 adjointManagers は、PrimalSolvers それぞれに対応するように作成する。
- primalSolver PrimalSolvers で定義した primal ソルバーの名前
- operatingPointWeight (デフォルト値:1） 複数の目的関数を取り扱う際にここで指定した値を重みとして使って結合することで最終的な一つの目的関数を得ることができる。
- adjointSolvers adjoint ソルバーの定義。複数の目的関数をそれぞれ別々に解くときには複数の adjoint ソルバーを設定する。この中で設定する項目名はすべての adjointManagers の中で重ならないようにつける必要がある。
- active (デフォルト値:true）この adjoint ソルバーを有効にするかどうか
- type adjoint ソルバーの解析対象のタイプ。現時点では「incompressible」しか使えない。
- solver adjoint 方程式を解くアルゴリズム。現時点では「adjointSimple」しか使えない。 ※元のマニュアルには、この項目も「type」と記載されていたがおそらくコードを見る限りでは「solver」を指していると思われるためこのように表記している。
- useSolverNameForFields (デフォルト値:false）primal ソルバーの useSolverNameForFields に対応する項目で、複数の adjoint ソルバー を定義する場合にはこのオプションを true にしなくてはならない。
- computeSensitivities (デフォルト値:true） adjoint 方程式を解いたあとに感度解析をするかどうか
- objectives 目的関数の設定
- type 目的関数のタイプ。現時点では「incompressible」しか使えない。
- objectiveNames 最小化したい目的関数のリスト。関数の中身は解析対象ケースの「objective/time_name/objective_name+adjoint_solver_name」内の「optimisation」フォルダの中に書かれる。この目的関数名はすべての adjointManagers の中で重ならないようにつける必要がある。
- type 最小化したい目的関数のタイプ。「force」「moment」「PtLosses」の中から選択する。詳細は下記。
- weight 目的関数の重み
- ATCModel adjoint 運動方程式に存在する Adjoint Transpose Convection 項を設定する項目。ATC は、刻み幅を極めて小さくしないと数値的に不安定になる stiff equation である。「standard」「UaGradU」「cancel」のいずれかもしくは(cancel, standard, UaGradU) のように複数を選択する。「standard」では、primal 速度ベクトル $\textbf{v}$ と、adjoint 速度ベクトル $\textbf{u}$ を用いて ATC 項を $u_j \dfrac{\delta v_j}{\delta x_i}$ と表す。primal 運動方程式における非保存形の移流項を微分することで定式化している。「UaGradU」では、ATC 項を$-v_j \dfrac{\delta u_j}{\delta x_i}$ と表す。こちらはprimal 運動方程式における保存形の移流項を微分している。「cancel」を指定すると、adjoint 運動方程式から ATC 項を除外して計算する。もちろん「cancel」を設定すると精度は落ちる。
- extraConvection (デフォルト値:0） adjoint 方程式を収束させるために、この値を adjoint 移流項から加減算する。
- zeroATCPatchTypes (デフォルト値:()）文字列のリストを指定する。ここで指定されたタイプの面の ATC 項をゼロにする。
- zeroATCZones (デフォルト値:()）zeroATCPatchTypes に似ているが、面ではなくセル領域を指定する。
- nSmooth (デフォルト値:0）zeroATCPatchTypes と zeroATCZones で指定されたセルを平滑化するためのラプラシアンフィルタを適用する回数
- maskType (デフォルト値:faceCells）ATC 項を平滑化するための隣接項の選び方。「faceCells」か「pointCells」を選択する。「faceCells」では境界面を共有するセルが選択され、「pointCells」では点を共有するセルが選択される。
- solutionControls primal ソルバーにおける「solutionControls」と同じ作用をする。平均化された adjoint 場は感度係数を計算するのに使う。
- printMaxMags (デフォルト値:false） ログに adjoint 場の最大値を出力するかどうかのオプション。シミュレーションの安定性を見たいときに有用。
目的関数の種類
force
$$J = \dfrac{\int_{S_W} \rho(-\tau_{ij}n_j+pn_i)r_idS}{\dfrac{1}{2} \rho AU_\infty^2}$$
- $\tau_{ij}$応力テンソルの要素
- $p$ 圧力を密度 $\rho$ で割ったもの
- $\vec{n}$ 単位法線ベクトル
- $\vec{r}$ 方向ベクトル。例えば、抗力を最小化するときには速度ベクトル方向に平行に設定する。
- $S_W$ 「force」が定義されている壁境界面
- $A$ 境界面の面積
- $U_\infty$ 速度の大きさ
objectives
{
type incompressible;
objectiveNames
{
drag
{
weight 1.;
type force;
patches ("wall.*" wallGroup);
direction (0.9993908270189999 0.034899496703 0);
Aref 2.;
rhoInf 1.225;
UInf 1;
}
}
}patches で指定しているように、ワイルドカードやグループ名も使える。
この目的関数は「最小化」を目指すものなので、最大化する目的関数を設定したいときには逆方向の値を与えなくてはならない。
moment
$$J = \dfrac{\int_{S_W} \rho r_i^M e_{ijk} (x_j- x_j^C) (-\tau_{kl} n_l + p n_k)dS}{\dfrac{1}{2} \rho AlU_\infty^2}$$
- $\vec{r^M}$ 最小化したいモーメント
- $\vec{x}$ 各境界面の位置ベクトル
- $\vec{x^C}$ 回転の中心となる位置ベクトル
- $l$ 基準長さ
- $e_{ijk}$ 交代記号（エディントンのイプシロン）
上記以外の項目は force の項と同じである。
moment
{
weight 1.;
type moment;
patches ("wall.*" wallGroup);
direction (0 0 1);
rotationCenter (0 0 0);
Aref 1.;
lRef 1.;
rhoInf 1.225;
UInf 6.;
}PtLosses
Pt は prestress ？
$$J = \int_{S_{I,O}} \left( p + \dfrac{1}{2} v_k^2 \right) v_i n_i dS$$
- $S_I, S_O$ 流入面と流出面
losses
{
weight 1;
type PtLosses;
patches (Inlet Outlet);
}optimisation
感度係数を導出する必要があるときにだけ設定する。
optimisation
{
sensitivities
{
type multiple;
patches (wall);
sensTypes
{
faces
{
type surface;
patches (wall);
includeSurfaceArea true;
includeObjectiveContribution true;
includeMeshMovement true;
adjointMeshMovementSolver
{
iters 300;
tolerance 1.e-7;
}
}
points
{
type surfacePoints;
patches (wall);
includeSurfaceArea true;
includeObjectiveContribution true;
includeMeshMovement true;
adjointMeshMovementSolver
{
iters 300;
tolerance 1.e-7;
}
includeDistance true;
adjointEikonalSolver
{
iters 300;
tolerance 1.e-7;
epsilon 0.1;
}
}
}
}
}- sensitivities adjoint ソルバーの設定で「computeSensitivities」をtrueにしたとき、全ての adjoint 方程式が解かれてから行う感度解析の設定。「type」と「patch」の指定が必須で、そこで選択した項目に応じて他のオプションを設定する必要がある。
- type 「surface」「surfacePoints」「sensitivityMultiple」の3つのうちから選択する。以下に追加で設定する項目について記載する。
- surface 境界壁面に対して垂直方向の変位に応じた目的関数の微分係数である「感度マップ」を計算する。計算時に faceSensNormal という名前の volScalarField に各 adjoint ソルバーの名前が追加されて、各タイムステップフォルダに書き込まれる。
- includeSurfaceArea (デフォルト値:true） 感度の値にローカルな面領域を含めるかどうか。大きな形状変化を起こしうる大きな面の移動があり、メッシュの解像度を考慮する必要がある場合には、true に設定する。逆に、正規化された感度分布が必要な場合は、感度が面領域の変化が一定になるように移動し場合、目的関数の変化はどのようになるかを示すように false に設定する。
- includeObjectiveContribution (デフォルト値:true） 例えば、抗力を最小化するときの法線ベクトルの変更のように、目的関数が感応度に対して直接的に影響することがある。これを true に設定することでこのような影響を計算に含む。
- includeMeshMovement (デフォルト値:true） グリッド変位スキームに adjoint 法を適用する際に生じる感度変化を考慮するかどうか。 true では拡張面積分、false ならば、面積分が使用される。感度の精度を上げるためには true に設定しておくと良い。 true では、グリッド変位変数の volVectorField「ma」に隣接する、周期境界条件やシンメトリーなどを除いたすべての境界条件を fixedValue の0に設定する。「ma」は基本的に自動的に生成される。また、「ma」用のソルバーを fvSolution に追加し、laplacian 離散化スキームを fvSchemes/laplacianSchemes に追加する。この方程式を解く際には緩和は必要ない。adjointMeshMovementSolver includeMeshMovement を true にした際に使用するソルバー。隣接するグリッド変位の収束を制御できる。
- iter (デフォルト値:1000） 最大反復回数
- tolerance (デフォルト値:1.e-06） 計算を終了する残差
- includeDistance (デフォルト値:true） Spalart-Allmaras 乱流モデルを微分するときに設定する項目。アイコナール方程式を解くかどうか。true にした場合には volScalarField「da」の境界条件を設定する必要がある。inlet と outlet は fixedValue の 0 に、 壁境界には zeroGradient に設定する。「da」は自動生成され、 fvSolution に「da」を解くためのソルバーと緩和係数を追加する必要がある。fvSchemes/divSchemes には
div(-yPhi, da)を追加する。adjointEikonalSolver includeDistance が true のときに設定。 iters と tolerance には adjointMeshMovementSolver と同じものを指定する。 epsilon は、デフォルト値が 0.1 だが、wallDist/advectionDiffusionCoeffs で method advectionDiffusion が指定されているときには fvSchemes で指定されているのと同じものを使う。不安定になっているときにはより大きい値を使うと良い。
- surfacePoints 境界面ではなく境界部分の点における法線方向の変位を元に感度解析をする。法線ベクトルの動きのような項は点を基準にして微分するときの方が適切な結果が出るため、感度マップを知りたいときは surfacePoints を選択する方が良い。計算時に pointScalarField「pointSensNormal」に各 adjoint ソルバーの名前が追加されて、各タイムステップフォルダに書き込まれる。設定項目は surface と同じ。
- multiple 複数のタイプの感度係数を解きたいときに指定する。sensTypes 以下で「surface」「surfacePoints」どちらかの項目を選んで設定する。
- surface 境界壁面に対して垂直方向の変位に応じた目的関数の微分係数である「感度マップ」を計算する。計算時に faceSensNormal という名前の volScalarField に各 adjoint ソルバーの名前が追加されて、各タイムステップフォルダに書き込まれる。
- patches どの境界面で感度解析をするか。ワイルドカードやグループ名を使うことができる。
- type 「surface」「surfacePoints」「sensitivityMultiple」の3つのうちから選択する。以下に追加で設定する項目について記載する。
2020/08/02 -
[…] adjointOptimisationFoam ユーザマニュアルを読んだ （前編） optimisationDictについて […]