OpenFOAMのESI版で粘弾性流体のCFDができるみたい

### 概要

OpenFOAMの主要なフォークには，ESI版/Foudation版/extend版の3つがあります．
これまで粘弾性流体のCFDができるのはextend版のみだと聞いていました．また，Web検索で見つかる情報も同じでした．調べたところ非常に基礎的な粘弾性流体モデルがESI版に入っていたことが判ったため投稿します．
粘弾性流体は非常に複雑な物性でそのモデルには各種あります．非常に基礎的な粘弾性流体モデルのみが実装されているESI版は“粘弾性流体のCFDができない”という表現は大局的に正しいと思います．

### 粘性と弾性
非常に簡単な説明をすると，粘性とは液体と気体を含む流体の性質で，弾性とは固体の性質です．
物体に何か外力を加えた時に物体は変形しその外力に対する抵抗が生じます．ある一定の大きさの外力を加え続ける時に，変形に抵抗する力の大きさが一定でいつまでも変形し続けるのが流体で，変形に抵抗する力の大きさが変形した量に応じて大きくなるため加えた外力と抵抗する力が釣り合った位置で変形が止まるのが固体です．

粘性という性質によって生じる抵抗力の方向は外力によって変形する方向のみで同じです．空間に３次元直交座標系(x1-x2-x3)をとり，外力による変形の方向がx1方向だとすると，抵抗力の成分はx1方向のみが非零の値でx2とx3の成分の値は零です．
弾性という性質によって生じる抵抗力は外力によって変形する方向以外にも生じます．３次元直交座標系(x1-x2-x3)において，外力による変形の方向がx1方向だとすると抵抗力のx2とx3の成分の値が常に零とは限りません．

### 粘弾性流体
一般的に流体は粘性流体として扱われます．（非粘性流体として扱う場合もありますが，ここでは粘弾性流体についての説明のため粘弾性流体に対する粘性流体の特別な場合と思ってください．）
粘弾性流体とは粘性と弾性の両方の性質を持った流体です．実は水も粘弾性流体ですが日常的な感覚では粘性の性質のほうが卓越しているため粘性流体として扱うことが一般的です．流体解析シミュレーション(CFD)でも一般的に水を粘性流体として扱って計算します．

粘弾性流体の特徴的な現象として，<a href="https://ja.wikipedia.org/wiki/%E3%83%AF%E3%82%A4%E3%82%BB%E3%83%B3%E3%83%99%E3%83%AB%E3%82%AF%E5%8A%B9%E6%9E%9C"> ワイセンベルク効果(Wikipedia) </a>や<a href="https://ja.wikipedia.org/wiki/%E3%83%90%E3%83%A9%E3%82%B9%E5%8A%B9%E6%9E%9C"> バラス効果(Wikipedia) </a>があります．気が付かないだけで，水でもワイゼンベルグ効果やバラス効果が起きます．水のワイゼンベルグ効果を利用した例に，<a href="https://youtu.be/P0gaRnFaScE"> 回転式のアロマディフューザー(YouTube) </a>があります．

### OpenFOAMのTutorialCase
チュートリアルケースは multiphase/compressibleInterFoam/laminar/climingRod で，回転する丸棒を液体が上昇する現象すなわちワイゼンベルグ効果を２次元軸対称モデルで計算しています．

compressibleInterFoamというカテゴリーは，名前からすれば圧縮性流体でのVOF法気液２相流です．このカテゴリーに含まれる他のケースは圧縮性の粘性流体のケースでした.（depthCharge2DとdeptheCharge3D．waterCoolerはコメント行によると開発中のテストケースです．調べたのはOpenFOAM V2206です．）

<dl>
<dd>
compressibleInterFoamという名前から，粘弾性流体に対応しているなんて気が付きませんでした．
</dd>
</dl>

climingRodの<b>constant/tubulentPropeteis</b>は下記のようになっています．depthCharge2DとdeptheCharge3DのケースではsimulationTypeがlaminarです．

```
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      turbulenceProperties;
}
simulationType  twoPhaseTransport;
```
climingRodには<b>constant/tubulentPropeteis.air</b>と<b>constant/tubulentPropeteis.liquid</b>もあります．この２つのファイルはdepthCharge2DとdeptheCharge3Dのケースにはありまません．

<b>constant/tubulentPropeteis.air</b>
```
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      turbulenceProperties;
}
simulationType  laminar;
```

<b>constant/tubulentPropeteis.liquid</b>
```
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      turbulenceProperties;
}
simulationType  laminar;
laminar
{
    model       Maxwell;

    nuM         0.01476;
    lambda      0.018225;
}
```

どうやら laminar {　} の中の Maxwelが粘弾性の設定のようです．
おそらく，粘弾性のマクスウェルモデルのことで，nuMは粘性係数，lamdaは緩和時間でしょう．

粘弾性モデルにはマクスウェルモデル，フォークトモデル，標準線形固体モデルがあります．参考資料のリンクを記載しておきます．
 * <a href="https://ja.wikipedia.org/wiki/%E7%B2%98%E5%BC%BE%E6%80%A7"> 粘弾性(Wikipwdia) </a>
 * <a href="https://en.wikipedia.org/wiki/Maxwell_material"> マクスウェル物性(en:Wikipwdia) </a>
 マクスウェル物性についての日本語のWikipediaのページはありませんが，連続体力学やレオロジーの教科書に載っています．

連続体力学やレオロジーの教科書の粘弾性の章には，ふつうは１次元モデルの図しか載っていません．３次元の計算を理解するには応力テンソルが必要です．この資料が良いかと思います
→　<a href="http://opencae.gifu-nct.ac.jp/pukiwiki/index.php?plugin=attach&refer=%C2%E8%A3%B1%A3%B7%B2%F3%CA%D9%B6%AF%B2%F1%A1%A7H240915&openfile=%C5%EC%B3%A417%B2%F3%40AY.pdf"> 粘弾性流体ソルバーviscoelasticFluidFoam～ 秋山善克</a>

粘弾性流体ソルバーviscoelasticFluidFoamを自分の使っているOpenFOAMのESI版では見たことがないため，この資料はextend版が前提なのでしょう．versionが2.1のようだからかなり古い資料だと思います．

この資料に記載されている“上対流微分マクスウェルモデル”が，ソースコードのコメントにも出てきました．

<b>scr/TurbulentModels/tubulentModels/lamina//Maxwell.H</b>から一部を抜粋
```
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2016 OpenFOAM Foundation
    Copyright (C) 2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
Description
    Maxwell model for viscoelasticity using the upper-convected time
    derivative of the stress tensor.
    See http://en.wikipedia.org/wiki/Upper-convected_Maxwell_model

    The model includes an additional viscosity (nu) from the transport
    model from which it is instantiated, which makes it equivalent to
    the Oldroyd-B model for the case of an incompressible transport
    model (where nu is non-zero).
    See https://en.wikipedia.org/wiki/Oldroyd-B_model

    Reference:
    \verbatim
        Amoreira, L. J., & Oliveira, P. J. (2010).
        Comparison of different formulations for the numerical calculation
        of unsteady incompressible viscoelastic fluid flow.
        Adv. Appl. Math. Mech, 4, 483-502.
        DOI:10.4208/aamm.10-m1010
    \endverbatim

SourceFiles
    Maxwell.C
\*---------------------------------------------------------------------------*/
```

Upper-convected_Maxwell_model が，さきほどの資料に出てきている上対流微分マクスウェルモデルのことです．
Oldroyd-B_model は，粘弾性モデルの一種です．（さきほどの資料にも出てきています．）

ソースコード中のコメントに記載されている<a href="https://en.wikiphttps://en.wikipedia.org/wiki/Upper-convected_Maxwell_modeledia.org/wiki/Upper-convected_Maxwell_model"> Upper-convected Maxwell model(en:Wikipedia) </a> のページに行くと，さきほどの資料と同じようなテンソル方程式が載っています．微妙に違っている理由は，せん断速度(ずり速度)をそのまま使って表記するか，せん断速度を変形速度を使って表記するかという違いです．

式に出てくる <b>*T*</b> は応力テンソルです．
scr/TurbulentModels/tubulentModels/lamina//Maxwell.Hのプログラム本体には
```
//- Return the Reynolds stress tensor
virtual tmp<volSymmTensorField> R() const;
```
と記載されています．

<dl>
<dd>
simulationTypeがlaminarで， RANS型乱流モデルの式で出てくる Reynolds stress (レイノルズ応力)が使われているのは不思議な感じがしますが……
</dd>
</dl>

<b>scr/TurbulentModels/tubulentModels/lamina//Maxwell.C</b> のほうにfvm法の数値計算の本体が記載されています．２相流のため計算式に表面張力 sigima も入ってきていて，修行不足で未熟者の私にはまだまだ理解できません．

### 実行結果
チュートリアルケースの multiphase/compressibleInterFoam/laminar/climingRod を./Allrunで実行するとワイゼンベルグ効果が発生します．depthCharge2Dのケースを参考にして改造するとワイゼンベルグ効果は発生しませんでした．

読者の方も“趣味レーション”してみてください．