2.2 穴あき板の応力解析
2.2.1.4 制御
通常どおり，解法の制御に関する情報は
controlDict
ディクショナリから読み込まれます．本 ケースでは，startTime
は0
です．本ケースは定常状態ですので，時間刻みは重要ではありま せん．このような状況では，定常状態のケースにおける反復回数カウンタとして働くよう，時間刻み
deltaT
を1
に設定するのが最善です．このようにした場合，本ケースで100
に設定した
endTime
は反復回数の上限として働きます．writeInterval
は20
に設定します．controlDict
のエントリは以下のようになります．17
18 application solidDisplacementFoam;
19
20 startFrom startTime;
21
22 startTime 0;
23
24 stopAt endTime;
25
26 endTime 100;
28 deltaT 1;
29
30 writeControl timeStep;
31
32 writeInterval 20;
33
34 purgeWrite 0;
35
36 writeFormat ascii;
37
38 writePrecision 6;
39
40 writeCompression off;
41
42 timeFormat general;
43
44 timePrecision 6;
45
46 graphFormat raw;
47
48 runTimeModifiable true;
49 50
51 // ************************************************************************* //
2.2.1.5
離散化スキームおよび線形方程式ソルバ制御次は
fvSchemes
ディクショナリについて見てみましょう．まず，この問題は定常状態ですので，
timeScheme
における時間微分としてはsteadyState
を選択します．これによって時間微 分項がオフの状態になります．全てのソルバが定常状態および過渡的状態の双方に対して適用 可能な訳ではありませんが，solidDisplacementFoam
は基本的なアルゴリズムが双方のシミュ レーションともに共通であるため，双方に適用可能となっています．線形弾性応力解析における運動方程式には，変位の勾配を含む陽な項がいくつか存在します．
この勾配を正確かつ滑らかに評価できれば，良い計算結果が得られます．通常，有限体積法にお ける離散化は，ガウスの定理に基づいています．ガウス法は大抵の目的においては十分に正確 ですが，本ケースにおいては最小二乗法を使用することとします．したがって
fvSchemes
ディ クショナリを開き，grad(U)
勾配離散化スキームとしてleastSquares
を選択してください．17
18 d2dt2Schemes
19 {
20 default steadyState;
21 }
22
23 ddtSchemes
24 {
25 default Euler;
26 }
27
28 gradSchemes
29 {
30 default leastSquares;
31 grad(D) leastSquares;
32 grad(T) leastSquares;
33 }
34
35 divSchemes
36 {
37 default none;
38 div(sigmaD) Gauss linear;
39 }
40
41 laplacianSchemes
42 {
43 default none;
44 laplacian(DD,D) Gauss linear corrected;
45 laplacian(DT,T) Gauss linear corrected;
46 }
47
48 interpolationSchemes
49 {
50 default linear;
51 }
52
53 snGradSchemes
54 {
55 default none;
56 }
57
58 fluxRequired
59 {
60 default no;
61 D yes;
62 T no;
63 }
64 65
66 // ************************************************************************* //
system
ディレクトリにあるfvSolution
ディクショナリでは，求解に使用される線形方程式のソルバおよびアルゴリズムを設定します．まず
solvers
サブディクショナリを見ると，D
のsolver
がGAMG
になっていることがわかります．tolerance
で表されるソルバ許容値（ソルバ名の次の 数値）は，本問題では10−6を設定します．relTol
で表されるソルバの相対許容値（さらにそ の次の数値）には各反復ごとの残差の所要低減量を設定します．本問題においては多くの項が 陽であり，また個別の反復的手順の一部としてアップデートされるため，各反復において厳し い相対許容値を設定することは非効率的です．したがって相対許容値として合理的な値は0.01， もしくはさらに高めの0.1，あるいはせいぜいこのケースのように0.9程度にしておきます．17
18 solvers
19 {
20 "(D|T)"
21 {
22 solver GAMG;
23 tolerance 1e-06;
24 relTol 0.9;
25 smoother GaussSeidel;
26 cacheAgglomeration true;
27 nCellsInCoarsestLevel 20;
28 agglomerator faceAreaPair;
29 mergeLevels 1;
30 }
31 }
32
33 stressAnalysis
34 {
35 compactNormalStress yes;
36 nCorrectors 1;
37 D 1e-06;
38 }
39 40
41 // ************************************************************************* //
fvSolution
ディクショナリは，アプリケーションソルバに特有の制御パラメータを含むstress-Analysis
サブディクショナリを含みます．まず，各時刻ステップ内での表面力境界条件処理を 含めた，全方程式系に関する外側ループの数を指定するnCorrectors
があります．本問題は定 常状態を扱いますので，「時刻ステップ」を反復回数カウンタとして使い収束解へと向かう反復 を実行することになります．したがってnCorrectors
を1
に設定します．D
キーワードには外側反復ループにおける収束許容値，すなわち初期残差に対して反復計算 によって消去されるべきレベルを設定します．本問題では前述において設定したソルバ許容値 の10−6に設定します．2.2.2
コードの実行以下に示すようなコマンドによって，実行後にログファイルに記録された収束状況を見るこ とができるよう，バックグラウンドでコードを実行します．
cd $FOAM_RUN/tutorials/stressAnalysis/solidDisplacementFoam/plateHole solidDisplacementFoam > log &
実行後には生成されたログファイルを見て，反復回数および解を求める各方向変位の初期・最 終残差などの収束状況を確認できます．本ケースの反復許容回数設定では，最終残差は必ず初 期残差の0.1倍以下となるはずです．いったん両初期残差ともに10−6の収束許容残差以下とな れば，その計算は収束したとみなしバッチジョブを
kill
することによって止めることができ ます．2.2.3
後処理後処理は2.1.4項と同様に行うことができます．
solidDisplacementFoam
ソルバは，応力場σ を対称テンソル場として出力します．したがって例えば，σxxをparaFoam
で見ることができます．
OpenFOAM
ソルバにおける変数名は通常，それらを表す数学記号にならって名付けられることは，ここで再度述べるに値するでしょう．ギリシア記号の場合は，変数は発音どおりに 名付けられます．例えば，σxxは
sigmaxx
と名付けられます．独立したスカラ成分のσxxやσxy などは，2.1.5.7で述べた
foamCalc
をsigma
に関して実行 して求めます．foamCalc components sigma
sigmaxx
やsigmaxy
などと名づけられた成分のファイルが時間のディレクトリに生成されます．図2.18のように応力を
paraFoam
で見ることができます．式 (2.14) の解析解とここで得られた数値解を比較しましょう．そのためには，解析領域左
端の対称面に沿ってのデータを出力しなければなりません．このようなグラフのために必要な データは，
sample
ユーティリティによって作成することができます．sample
の設定はsystem
ディレクトリ内のsampleDict
で行い詳細は表6.3に要約しています．データのサンプリングを 行う座標区間は，sets
によって(0.0,0.5,0.25)から(0.0,2.0,0.25)に指定されています．物理 量はfields
に指定します．0 5 10 15 20 25 30
σ
xx(kP a)
図2.18 穴あき板における応力場
17
18 interpolationScheme cellPoint;
19
20 setFormat raw;
21 22 sets
23 (
24 leftPatch
25 {
26 type uniform;
27 axis y;
28 start ( 0 0.5 0.25 );
29 end ( 0 2 0.25 );
30 nPoints 100;
31 }
32 );
33
34 fields ( sigmaxx );
35 36
37 // ************************************************************************* //
通常通り
sample
を実行してください．writeFormat
はraw
形式で2
列のフォーマットとなっ ています．データはsets
ディレクトリの時刻サブディレクトリの中のファイルに書き込まれ ます．たとえば，時刻t= 100 sのデータはsets/100/leftPatch_sigmaxx.xy
に書き込まれます．GnuPlot
のようなアプリケーションでは，コマンドプロンプトで以下を入力することで，数値解および解析解の両方をプロットすることができます．
plot [0.5:2] [0:] ’sets/100/leftPatch_sigmaxx.xy’, 1e4*(1+(0.125/(x**2))+(0.09375/(x**4)))
プロット例を図2.19に示します．
2.2.4
演習以下は
solidDisplacementFoam
に習熟していただくための演習課題です．0 5 10 15 20 25 30 35
0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
Stress ( σ
xx)
x=0(kP a)
Distance, y (m)
Numerical prediction Analytical solution
図2.19 垂直方向対称面における法線方向応力
2.2.4.1
メッシュ解像度の増加x，y方向それぞれのメッシュ解像度を増やしてみましょう．2.2.3項の最終的な粗メッシュ の結果を，
mapFields
を使って密メッシュの初期条件にマッピングしてください．2.2.4.2
非等間隔メッシュの導入穴に近いセルが遠いセルより密になるよう，メッシュ幅を変化させてください．隣接するセ ルの大きさの比率が1.1以上にならないように，またブロック間のセルの大きさの比率がブロッ ク内の比率と同様となるようメッシュを作成してください．メッシュの非等間隔化については
2.1.6項で述べました．ここでも，2.2.3項の最終的な粗メッシュでの結果を，
mapFields
を使って非等間隔メッシュの初期条件としてマッピングします．結果を解析解および非等間隔化する 前の結果と比較してみましょう．等間隔メッシュと同一のセル数を使用した場合，解の精度は 改善されるでしょうか？