前回の記事「Elmerを使ったパラメトリック解析」はこちら
異方性材料の物性値設定
伝熱解析、構造解析、電磁界解析など、Elmerのモジュールの中には異方性材料をサポートしているものがあります。今回は伝熱解析で異方性の物性値を与える例で、パラメトリック解析を行いたいと思います。具体的にはFORTRANのプログラムで熱伝導率テンソルの主軸を回転させて、材料内で熱の通りやすい方向を色々変えてみるというものです。
今回のモデルはL字型の部品を考えます。ジオメトリは前回の記事のようにチュートリアルのものを使用できればよかったのですが、境界の設定が少々面倒だったため、Salomeを使ってUNVファイルで作成した以下のようなモデルを用いました。Elmerが対応しているメッシュ、ジオメトリファイル形式で作成できれば何でも大丈夫です。
ついでに境界条件を説明すると、L字の下部に一定の熱流束（100.0）を、それ以外の端部には室温（20.0）空気中での熱伝達（熱伝達率10.0）を与えました。（数値はテキトーです）
各種解析設定の方法は前回の記事をご覧いただくとして、今回の記事の目玉となるFORTRANサブルーチンを使ったパラメトリック解析を行っていきます。
まず、今回の解析のパラメータである熱導電率テンソルの主軸の角度は以下のように、1°から180°へ1°ずつ変えていきます。
[Timestep Intervals] -> 180
[Timestep sizes] -> 1.0
次に、熱伝導率を以下のように記述します。
[Heat Conductivity] -> Variable Time; size 2 2; Real Procedure "cond" "getHeatCond"
このsize 2 2;は以降の値が2×2の配列であることを示すElmerのキーワードで、Real Procedure "cond" "getHeatCond"はcond.dll（cond.f90をコンパイルしたもの）のgetHeatCondサブルーチン（または関数）を呼び出すことを示しています。
次に、熱伝導率テンソルを計算するサブルーチンをプログラムしていきます。
SUBROUTINE getHeatCond(model, n, t, cond)
USE DefUtils
IMPLICIT None
TYPE(Model_t) :: model
INTEGER :: n
REAL(KIND=dp) :: t
REAL(KIND=dp),pointer,contiguous :: cond(:, :)
REAL(KIND=dp) :: ref_c(2, 2)
REAL(KIND=dp) :: rad
REAL(KIND=dp) :: R(2, 2)
ref_c = reshape((/100.0, 0.0, 0.0, 1.0/), shape(ref_c))
rad = t / 180.0 * acos(-1.0)
R = reshape((/cos(rad), -sin(rad), sin(rad), cos(rad)/), shape(R))
cond = matmul(matmul(TRANSPOSE(R), ref_c), R)
END SUBROUTINE getHeatCond
細かい説明は省きますが、パラメータt（角度）を使って、主軸を回転させた熱伝導率テンソルを計算し、condとして返しています。condのcontiguous属性は無しで書く人もいるみたいですが、自分の場合invalid memory referenceのエラーが出てしまったので、対応としてつけたら動くようになりました。
次に、cond.f90のあるフォルダでコマンドラインを開き、elmerf90 cond.f90 -o cond.dllを実行すると、FORTRANのプログラムがコンパイルされます。（もしかしたら最初はelmerf90のディレクトリにPATHが通っていないかもしれないです）
cond.dllが作成されたら準備はOKなので解析を実行します。
解析結果
角度を変えると、定常状態はこんな感じに変わります。
せっかくなので、最高温度をプロットしてみましょう。
角度を20°、150°くらいにすると最高温度を抑えられるらしい。（テキトーな設定なのでテキトーな結果）
最後に
今回はElmerを使ったパラメトリック解析の応用編として、FORTRANの自作サブルーチンを使った解析を行いました。正直に言えば今回の解析だとFORTRANを使わずElmerのMATCプログラムを使っても記述できましたが、より複雑な数学的計算を含む場合は今回のような方法が役に立つと思います。