有限要素法の実装(1) ひずみ変位行列(2)

　前回は平面応力・平面ひずみの$B$マトリックスの式を展開したので、今回は軸対称固体要素の説明を行う。形状関数は前回同様パラメトリック変換要素を仮定するが、形状関数については特に仮定を設けないことにする。

　軸対称問題の座標系は直交座標系の$x,y,z$に対して$r,\theta,z$になる。それぞれ、半径、角度、高さである。モデル自身は$r,z$平面上で作成される。ソルバーでの座標系の取り方は$z$軸をを正直に全体座標系の$z$軸とするか、$x,y$平面上に作成するかのいずれかである。この場合、$x$軸と$y$軸は円筒座標系の$r$軸と$z$軸になる。前者のメリットは座標系の定義がそのまま使えることだが、後者は二次元問題のコードをそのまま再利用できる(これが一番大きい)ためだが、最終的に3次元形状への投影を行うにしても大した手間がかかるわけではないので、恐らくは$x,y$平面上でモデルを作るようにしている方が多いと思われる。
　前者の$z$軸をそのまま使用した例としてはI-deasがあるが他にもあるかもしれない。

　$\mathbf B$マトリックスの導出に入る。初めに軸対称体のひずみベクトルは次のような形で表される。

\mathbf \epsilon
=
\left [
\begin{matrix}
\epsilon_r \\
\epsilon_\theta \\
\epsilon_z \\
\gamma_{rz}
\end{matrix}
\right ]
=
\left [
\begin{matrix}
\frac {\partial u} {\partial r} \\
\frac {u} {r} \\
\frac {\partial w} {\partial z } \\
\frac {\partial u} {\partial z} +
\frac {\partial w} {\partial r}
\end{matrix}
\right ]

これを前回同様$\partial u_i/ \partial x_j$を一列に並べた形に変形する。このためには平面問題同様符号行列$\mathbf \Lambda$を使用する。

\left [
\begin{matrix}
\frac {\partial u} {\partial r} \\
\frac {u} {r} \\
\frac {\partial w} {\partial z} \\
\frac {\partial u} {\partial z} +
\frac {\partial w} {\partial r}
\end{matrix}
\right ]
=
\left [
\begin{matrix}
1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 \\
0 & 0 & 1 & 0 & 0 \\
0 & 1 & 0 & 1 & 0
\end{matrix}
\right ]

\left [
\begin{matrix}
\frac {\partial u} {\partial r} \\
\frac {\partial u} {\partial z} \\
\frac {\partial w} {\partial z} \\
\frac {\partial w} {\partial r} \\
\frac {u} {r}
\end{matrix}
\right ]

=\mathbf \Lambda \mathbf d

　$\mathbf \Lambda$と$\mathbf d$は前回と同じ記号を用いているが、項の中身は異なる。$\mathbf d$に半径方向のひずみが追加されるので、これに伴って$\mathbf \Lambda$も変更されるからである。まず$\mathbf d$のうち4行目迄は前回と同じになるからこの項の展開は省略することにする。今回の場合、座標系は二次元平面の$x,y$を$x,z$に変更し、変位は$u,v$を$u,w$に変更するだけである。
　$\mathbf d$の4行目までの最終的な形は次のようになる。

\left [
\begin{matrix}
\frac {\partial u} {\partial r} \\
\frac {\partial u} {\partial z} \\
\frac {\partial w} {\partial z} \\
\frac {\partial w} {\partial r}
\end{matrix}
\right ]
=
\left [
\begin{matrix}
\mathbf A & 0 \\
0 & \mathbf A 
\end{matrix}
\right ]

ここで上の各項は次のようになる。なお、ヤコビアンは二次元平面問題と同じである。

\mathbf A = \mathbf J^{-1} \mathbf H = 
\left [
\begin{matrix}
\mathbf A_1 \\
\mathbf A_2
\end{matrix}
\right ] \\
\mathbf H =
\left [
\begin{matrix}
\frac {\partial N^i}{r_1} \\
\frac {\partial N^i}{r_2}
\end{matrix}
\right ] \\
\mathbf J = \mathbf H \mathbf X \\
\mathbf X =
\left [
\begin{matrix}
r^i & z^i
\end{matrix}
\right ] \\

　$\mathbf d$の最期の項$u/r$は次のようにして求める。まず$r$は半径方向の積分点座標となるから$r=N^ir^i$となる。ここで$N^i$は$i$番目の形状関数であり、$r^i$は同じく$i$番目のコネクティビティに従った節点の$r$方向座標である。
分母の$u$は座標値と同じように$u$=$N^iu^i$の形になる。ここで注意することは分母は半径での総和を必要とするが、分子は各項ごとに離散化するため、分母と異なり総和を取らない。
　よって$u/r$は次のようにになる。ここで、$\mathbf R = N^i / N^j r^j$で、記法の略記に為に用いる。

\frac {u} {r} = \frac {N^i u^i} {N^j r^j} = 
\left [
\begin{matrix}
\frac {N^i} {N^j r^j}
\end{matrix}
\right ]

\left [
\begin{matrix}
u^i
\end{matrix}
\right ]
= \mathbf R \left [
\begin{matrix}
u^i
\end{matrix}
\right ]

これらを合成すると次のようになる。

\mathbf d =
\left [
\begin{matrix}
\mathbf A & 0 \\
0 & \mathbf A \\
\mathbf R & 0
\end{matrix}
\right ]

　最終的な$B$マトリックスは次のようになる。

\mathbf B = \mathbf \Lambda \mathbf d\\
=
\left [
\begin{matrix}
1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 \\
0 & 0 & 1 & 0 & 0 \\
0 & 1 & 0 & 1 & 0
\end{matrix}
\right ]

\left [
\begin{matrix}
\left [
\begin{matrix}
\mathbf A_1 \\
\mathbf A_2
\end{matrix}
\right ] & 0\\
0 & \left [
\begin{matrix}
\mathbf A_1 \\
\mathbf A_2
\end{matrix}
\right ] \\
\mathbf R & 0
\end{matrix}
\right ]

\\ =
\left [
\begin{matrix}
\mathbf A_1 & 0 \\
\mathbf R & 0 \\
0 & \mathbf A_2 \\
\mathbf A_2 & \mathbf A_1
\end{matrix}
\right ]

　実際のプログラムを示す。ヤコビアンの計算は次のようになる。

fem.c

void calc_Jmat (double *J, const double *H, const Element *e)
{
	int i, j, k;
	int ndim = etinfo [prob].ndim, nnode = e->info->info1->nnode;
	memset (J, 0, sizeof (double) * ndim * ndim);
	for (i = 0; i < ndim; i++)
		for (j = 0; j < ndim; j++)
			for (k = 0; k < nnode; k++)
				J [i * ndim + j] += H [i * nnode + k] * node [e->conn [k]].coord [j];
}

　ndimは要素が定義される次元であり、ここでは2を返している。nnodeは要素辺りの節点数で、これは要素の形状と字数によって決められる。これらを管理するのが、要素タイプ構造体ElementTypeInfoと要素情報構造体ElementInfoである。ElementInfoはBマトリックスの計算に必要なデーターの他、表面力や体積力の情報も有している。これらはどれも形状関数の情報と積分計算に関する情報から成っているので、構造体とするのが一番便利である。この構造体が、ElemIntegralInfoというもので、計算に必要な諸量を格納している。
　これらの構造体の定義は次のようになっている。

fem.c

typedef struct {
	int number;
	char *name;
	int prob, ndim, ndof, nvstat, ijh_scale, ijh_add;
} ElementTypeInfo;

typedef struct {
	int number;
	int nnode, norder, ntint, ndim; // 要素あたりの節点数、次数、全積分点数、次元数
	SF sf, dsf;   // 形状関数へのポインター
	double *N, *dN, *ipt, *wt; // 形状関数と形状関数の微分形、積分点情報と重みの配列
} ElemIntegralInfo;

typedef struct {
	int number;
	char *name;
	int nvert, nface;  // 要素あたりの節点数、頂点(隅節点)の数、表面の数
	ElemIntegralInfo *info1, *info2, *info3, *info4; // 要素全体と表面の計算情報
	int *edge_list;
} ElementInfo;

　形状関数の値はElemIntegralInfoのNとdNに格納される。これがElementInfoから引用されるのである。ElementInfoはElemIntegralInfoを4個持っているが、Bマトリックスの計算のために2個を確保している。実際には2個目にあたるinfo2は未使用だが、将来の拡張のために保持している。3個目以降は面荷重(2次元要素では辺荷重)と体積力のためにあり、辺荷重用は線要素の形状関数データーを保持している。体積力は完全積分要素の形状関数データーを保持するが、荷重計算に関しては辺荷重と共に完全積分を行うこととされているので、Bマトリックス用と別にデーターを持たせている。これによって次数低減積分要素の実装はそれ用のデーターを別途用意してinfo1に入れておけば完成する。
　これらの構造体はプログラムの開始に先立って初期化されるので、要素剛性マトリックスの計算中に形状関数値を計算する必要はなくなる。このため以前示したように要素剛性の計算の際にはどの要素も積分点の個数だけループを回せばよく、四角形要素や六面体の実相例でよくあるような各座標ごとのループを回す必要はなくなる。

nodeは節点の情報を持つ構造体の配列で、メンバーのcoordは座標値となる。これらの構造体を要ししておけば後は式を立てたとおりに書くだけの作業になる。なお、この関数は3次元固体要素にそのまま転用可能である。座標値をひとまとめにしておけばループ本体は1行で書くことが出来る。こうすれば確実に動くものが得られる可能性が上がるのである。

　mat_mul2は2次元配列同士の積を計算する関数である、結果はAに入る。Cの二次元配列は扱いづらいので1次元配列の形にしている。この関数はここではAマトリックスの計算に用いられる。

fem.c

void mat_mul2 (double *A, const double *B, const double *C, int sz1, int sz2, int sz3)

{
	int i, j, k;
	memset (A, 0, sizeof (double) * sz1 * sz3);
	for (i = 0; i< sz1; i++)
		for (j = 0; j < sz2; j++)
			for (k = 0; k < sz3; k++)
				A [i * sz2 + j] += B [i * sz3 + k] * C [k * sz2 + j];
}

　ヤコビアンの行列値とAマトリックスの値は最終的にひずみと応力の計算に用いるため格納しておく必要がある。このための構造体、要素情報構造体は次のようになる。この構造体は要素データー内で積分点数の配列として割り当てられる。

fem.c

typedef struct {
	double detJ, sigy, *iJH;  // ヤコビアンの行列値とJ-1*Hの結果。
	double *stress, *strain, *FT; // 応力とひずみ、変形勾配
	double *stat;   // 状態変数(温度など)
} ElementCalcInfo;

　最後にヤコビアンの逆行列を計算する関数を用意する。ここでは行列値を関数の中で求めている。これも2次元配列ではなく1次元配列の形にしている以外は公式通りである。

fem.c

double invJ_2D (double *iJ, const double *J)
{
	double detJ;
	detJ = J [0] * J [3] - J [2] * J [1];
	iJ [0] =  J [3] / detJ;
	iJ [1] = -J [1] / detJ;
	iJ [2] = -J [2] / detJ;
	iJ [3] =  J [0] / detJ;
	return detJ;
}

　ここまで用意できれば、あとはBマトリックスの計算を行うだけになる。まず平面要素用の関数を示す。以前に説明したとおりになっていることがわかるだろう。これらの関数は積分点ごとに呼び出される。mat_mul2は先にも書いたとおり、Aマトリックスの計算のために呼び出される。

fem.c

void calc_PL_B_mat (double *B, const Element *e, const double *N, const double *H, int nint)
{
	int i;
	int ndim = etinfo [prob].ndim, nnode = e->info->info1->nnode;
	ElementCalcInfo *cp = e->cinfo + nint;
	int size = ndim * nnode;
	double J [ndim * ndim], iJ [ndim * ndim];

	calc_Jmat (J, H, e);
	cp->detJ = invJ_2D (iJ, J);
	mat_mul2 (cp->iJH, iJ, H, ndim, nnode, ndim);
	// Bマトリックスの作成。一つ飛ばし(つまり自由度分)に値を代入している。
	for (i = 0; i < nnode; i++) {
		B [0 * size + i * 2 + 0] = B [2 * size + i * 2 + 1] = cp->iJH [i];
		B [2 * size + i * 2 + 0] = B [1 * size + i * 2 + 1] = cp->iJH [i + nnode];
		B [1 * size + i * 2 + 0] = B [0 * size + i * 2 + 1] = 0;
	}
}

　次が軸対称要素用の関数である。計算情報の大きさを弄ってこのような形にしているが、iJHはnnode*2個までが二次元平面要素と同じであり、この後のnnodeが式中の$\mathbf R$の分である。最後に積分点における半径の座標が格納される仕組みになっている。途中までは二次元平面要素の計算と同じだが、この後に半径を計算し、BマトリックスにAマトリックスと形状関数値を半径で割った値を格納している。

fem.c

void calc_AXSOL_B_mat (double *B, const Element *e, const double *N, const double *H, int nint)
{
	int i;
	int ndim = etinfo [prob].ndim, nnode = e->info->info1->nnode;
	int size = ndim * nnode;
	ElementCalcInfo *cp = e->cinfo + nint;
	double J [ndim * ndim], iJ [ndim * ndim];
	double cx;

	calc_Jmat (J, H, e);
	cp->detJ = invJ_2D (iJ, J);
	mat_mul2 (cp->iJH, iJ, H, ndim, nnode, ndim);

	// 積分点での半径を求める。
	cx = 0;
	for (i = 0; i < nnode; i++)
		cx += N [i] * node [e->conn [i]].coord [0];
	cp->iJH [nnode * 3] = cx;
	for (i = 0; i < nnode; i++) {
		B [0 * size + i * 2 + 0] = B [3 * size + i * 2 + 1] = cp->iJH [i];
		B [3 * size + i * 2 + 0] = B [1 * size + i * 2 + 1] = cp->iJH [i + nnode];
		B [2 * size + i * 2 + 0]                            = cp->iJH [i + nnode * 2] = N [i] / cx;
		B [1 * size + i * 2 + 0] = B [0 * size + i * 2 + 1] = B [2 * size + i * 2 + 1] = 0;
	}
}

　これらの関数ではループを多用しているが、こうすれば初めから正しい結果が得られる公算が大きいからである。これらの関数を書くにあたりご丁寧に一つ一つばらしている例があるが、そのような手法は関数がいたずらに長くなり、書くのも大変になる。当然間違える確率も上がり、理解も遠のくだけである。最も、その前に集中力が切れそうだが。実際のプロダクトでは効率も重視されるから、のような手法を採用することも考えられるが、初心者に教える場合はこのような形にした方がよい。