OpenFOAM の adjoint Shape Optimization Foam を考える #04境界条件の実装

前回に引き続き、adjoint法を用いてトポロジー最適化を行う際の実装を見てみます。
今回は境界条件がどのように設定されているのかを見ます。

コードの行数と内容はOpenFOAMv1812を参考にしています。

前回は「OpenFOAM の adjoint Shape Optimization Foam を考える #03実装」でした。

1. OpenFOAM の adjoint Shape Optimization Foam を考える #01手法
2. OpenFOAM の adjoint Shape Optimization Foam を考える #02境界条件
3. OpenFOAM の adjoint Shape Optimization Foam を考える #03実装
4. OpenFOAM の adjoint Shape Optimization Foam を考える #04境界条件の実装 （今回）

境界条件の実装

adjointShapeOptimizationFoamの実装では、圧力損失の最小化を行っています。
OpenFOAM の adjoint Shape Optimization Foam を考える #02境界条件 にて、流出部分での流量の均一化について述べましたが、その実装はされていません。

今回は、特にupdateCoeffs()に注目して、どのように値を更新するかを見ていこうと思います。

流入部分の境界条件は、通常の速度や圧力の方程式と同じように設定されています。
流出部分での境界条件は、adjoint法に固有のものが追加で実装されています。

adjointOutletPressureFvPatchScalarField.C

adjoint pressureを導出するupdateCoeffs()は以下のように実装されています。

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::adjointOutletPressureFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const fvsPatchField<scalar>& phip =
        patch().lookupPatchField<surfaceScalarField, scalar>("phi");

    const fvsPatchField<scalar>& phiap =
        patch().lookupPatchField<surfaceScalarField, scalar>("phia");

    const fvPatchField<vector>& Up =
        patch().lookupPatchField<volVectorField, vector>("U");

    const fvPatchField<vector>& Uap =
        patch().lookupPatchField<volVectorField, vector>("Ua");

    operator==((phiap/patch().magSf() - 1.0)*phip/patch().magSf() + (Up & Uap));

    fixedValueFvPatchScalarField::updateCoeffs();
}

90行目から100行目では、lookupPatchField関数で変数を検索することによって、必要なfvsPatchFieldsとfvPatchFieldsを定義しています。
102行目で演算子の定義をして、adjoint pressure $pa(q)$ に値を割り当てます。
ここでは計算対象の面領域ベクトルの大きさを patch().magSf の定義を使用して得られる

[1]　$q=(u_n-1)v_n+\textbf{u}\cdot\textbf{v}$

の値を用いています。

流出部分の値をadjoint pressureに割り当てるために解く方程式は、

[2]　$q = \textbf{u}\cdot\textbf{v} + u_n v_n + \nu (\textbf{n}\cdot\nabla ) u_n + \frac{\partial J_\Gamma}{\partial v_n}$

です。

現在の実装では、OpenFOAM の adjoint Shape Optimization Foam を考える #02境界条件 で見たような、

[3]　$q = \textbf{u}\cdot\textbf{v} + u_n v_n + \nu (\textbf{n}\cdot\nabla ) u_n + c(v_n – v_n^d)$

[4]　$0 = v_n \textbf{u}_t + \nu (\textbf{n}\cdot\nabla ) \textbf{u}_t + c(\textbf{v}_t – \textbf{v}_t^d)$

で与えられる境界条件とは異なる境界条件の式を用いているようです。

adjointOutletVelocityFvPatchVectorField.C

adjoint velocity を導出するupdateCoeffs()は以下のように実装されています。

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

// Update the coefficients associated with the patch field
void Foam::adjointOutletVelocityFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const fvsPatchField<scalar>& phiap =
        patch().lookupPatchField<surfaceScalarField, scalar>("phia");

    const fvPatchField<vector>& Up =
        patch().lookupPatchField<volVectorField, vector>("U");

    scalarField Un(mag(patch().nf() & Up));
    vectorField UtHat((Up - patch().nf()*Un)/(Un + small));

    vectorField Uan(patch().nf()*(patch().nf() & patchInternalField()));

    vectorField::operator=(phiap*patch().Sf()/sqr(patch().magSf()) + UtHat);
    // vectorField::operator=(Uan + UtHat);

    fixedValueFvPatchVectorField::updateCoeffs();
}

adjoint velocity の境界条件は adjoint pressure と同じように実装されています。
速度の接線方向成分の計算には、

[5]　$0 = v_n \textbf{u}_t + \nu (\textbf{n}\cdot\nabla ) \textbf{u}_t + \frac{\partial J\Gamma}{\partial \textbf{v}_t}$

が用いられ、面全体のadjoint velocity は $\textbf{u}=\textbf{u}_t+\textbf{u}_n$ として求まります。
これはコードの102行目に対応しています。

UtHatを数式で表すと下記のようになります。

[6]　$\textbf{u}_{p,t} = \frac{\textbf{v}_p-\textbf{v}_{p,n}}{u_{p,n} + SMALL}$

ここで、添字${}_p$ は現在の中心ノード、${}_n$ は法線要素ノード、${}_t$ は接線要素を示し、微小な値 SMALL は、分母がゼロにならないように足す値です。

このように実装されている境界条件に対応するコスト関数がどのようなものなのかについては、式[1] や式[6] を見ただけでは簡単には分かりません。

100行目で定義しているものの使われていない変数 $Uan$ を見ると、patchInternalField()関数を使用して、隣接ノードから得られる値にアクセスし、tmp<Field<Type>>を生成します。
tmpクラスを使用すると、値をコピーせずに関数から返すことができるため、メモリの節約ができます。

隣接ノードにおける速度は、$\nu(\textbf{n}\cdot\nabla)\textbf{u}_t$ 項が、対象領域と隣接ノード間の速度差を距離で割った値として近似されるときに必要になります。
このあたりの情報は今後コスト関数を変更し、境界条件を設定する際に使います。

まとめ

今回の実装は、OpenFOAM の adjoint Shape Optimization Foam を考える #02境界条件 に対応しています。

今後は圧力損失以外の境界条件の設定について見ていこうと思います。