GrasshopperでFEM解析をイチから実装する その3平面応力要素 四角形1次非適合要素編

はじめに

前回までの記事では、四角形1次要素の理論やGrasshopper上での実装方法について説明しました。今回はその続きとして、四角形1次要素のせん断ロッキングを解決する手法である非適合要素について、同様に理論を説明して実装、検証を行います。

コンポーネント群は　Github　で公開しています。

コメントなどリアクションいただけると嬉しくて記事を書くモチベーションになるので、ぜひお気軽にお願いします。

非適合要素とは

非適合要素とは、曲げ変形を生じさせたときに発生するせん断ロッキングと呼ばれる現象を解決するために曲げ変形を表現できるようにした要素です。
せん断ロッキングとは、1次要素において曲げ変形が生じた際、要素内で変位を線形で仮定するため、変形が台形となり、剛性が過度に高くなる現象です。これを解決するため、1次要素の変位場に曲げによる項を追加し、より現実的な変形を仮定します。

要素ごとの連続性は節点位置でのみ保証され、曲げ変形の曲率は要素ごとに異なり、要素辺は隣接する要素で適合しません。このことから、非適合要素と呼ばれています。1次要素の拡張的な要素のため節点数4・積分点4のままで、2次要素として曲げ変形を表現するよりも、計算負荷を抑えることができます。２次要素については、後日記事にする予定です。

ちなみに、建築系の方は使用している方も多いと思いますが、汎用構造解析ソフト「midas iGen」の2次元要素では、この非適合要素が用いられています。

定式化

非適合要素では、前述の通り変位に曲げ変形による項を加えます。式は以下の通り。
$$
\begin{gather}
u(\xi,\eta) = \sum_{i=1}^{4}Ni(\xi,\eta)ui+\alpha_1(1-\xi^2)+\alpha_2(1-\eta^2)\\
v(\xi,\eta) = \sum_{i=1}^{4}Ni(\xi,\eta)vi+\alpha_3(1-\xi^2)+\alpha_4(1-\eta^2)
\end{gather}
$$

ここで、形状関数$N_i$は、四角形１次用要素で登場したものと同じで、下記のとおり。
$$
\begin{gather}
N_1 = \frac{1}{4}(1-\xi)(1-\eta) \\
N_2 = \frac{1}{4}(1+\xi)(1-\eta) \\
N_3 = \frac{1}{4}(1+\xi)(1+\eta) \\
N_4 = \frac{1}{4}(1-\xi)(1+\eta)
\end{gather}
$$

この$\alpha$のついた項が、曲げ変形を表すために追加された項で、曲げ変形を２次式で表しています。前述のとおり、この$\alpha$は要素ごとに異なる値を持つので、要素辺の連続性が担保されず非適合となります。
これを整理すると、
$$
\begin{gather}
&\mathbf{u} = \mathbf{Nu_e} + \mathbf{P\alpha} \\
&\mathbf{N} =
\begin{bmatrix}
N1 & 0 & N2 & 0 & N3 & 0 & N4 & 0 \\
0 & N1 & 0 & N2 & 0 & N3 & 0 & N4
\end{bmatrix}
\hspace{10pt}
&\mathbf{ue} =
\begin{bmatrix}
u1 \\
v1 \\
u2\\
v2\\
u3\\
v3\\
u4\\
v4
\end{bmatrix} \\
&\mathbf{P} =
\begin{bmatrix}
1-\xi^2 & 1-\eta^2 & 0 & 0 \\
0 & 0 & 1-\xi^2 & 1-\eta^2
\end{bmatrix}
\hspace{10pt}
\mathbf{\alpha} =
\begin{bmatrix}
\alpha_1 \\
\alpha_2 \\
\alpha_3 \\
\alpha_4 \\
\end{bmatrix}
\end{gather}
$$

よって歪み$\epsilon$は、以下のように表せる。

$$
\mathbf{\epsilon} =
\begin{bmatrix}
\epsilon_x \\
\epsilon_y \\
\tau_xy
\end{bmatrix}=
\begin{bmatrix}
\frac{\partial u}{\partial x} \\
\frac{\partial v}{\partial y} \\
\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \\
\end{bmatrix}=
\mathbf{B_c}\mathbf{u_e} +
\mathbf{B_I}\mathbf{\alpha}
$$
ここで、
$$
\mathbf{Bc} =
\begin{bmatrix}
\frac{\partial N_1}{\partial x} & 0 &
\frac{\partial N_2}{\partial x} & 0 &
\frac{\partial N_3}{\partial x} & 0 &
\frac{\partial N_4}{\partial x} & 0 \\
0 &
\frac{\partial N_1}{\partial y} & 0 &
\frac{\partial N_2}{\partial y} & 0 &
\frac{\partial N_3}{\partial y} & 0 &
\frac{\partial N_4}{\partial y} \\
\frac{\partial N_1}{\partial y} & \frac{\partial N_1}{\partial x} &
\frac{\partial N_2}{\partial y} & \frac{\partial N_2}{\partial x} &
\frac{\partial N_3}{\partial y} & \frac{\partial N_3}{\partial x} &
\frac{\partial N_4}{\partial y} & \frac{\partial N_4}{\partial x}
\end{bmatrix}
$$

$$
\mathbf{B_I} =
\begin{bmatrix}
-2\xi\frac{\partial \xi}{\partial x} &
-2\eta\frac{\partial \eta}{\partial x} &
0&0 \\
0&0&-2\xi\frac{\partial \xi}{\partial y} &
-2\eta\frac{\partial \eta}{\partial y} \\
-2\xi\frac{\partial \xi}{\partial x} &
-2\eta\frac{\partial \eta}{\partial x} &
-2\xi\frac{\partial \xi}{\partial y} &
-2\eta\frac{\partial \eta}{\partial y} \\
\end{bmatrix}
$$
$B_c$は、四角形１次要素の$B$マトリクスと同じで、$B_I\alpha$が非適合要素の追加分となります。

仮想仕事の原理を用いて、要素剛性マトリクスを計算すると、以下４つのマトリクスが得られる。

$$
\begin{gather}
K_{cc} = t\int{B_C^{\intercal}DB_C}dA,K_{ci} = t\int{B_C^{\intercal}DB_I}dA
\\
K_{ic} = t\int{B_I^{\intercal}DB_C}dA,K_{ii} = t\int{B_C^{\intercal}DB_i}dA
\end{gather}
$$
この４つの剛性マトリクスは以下の関係を持つ。
（内力による仮想仕事を０とみなしている。)
$$
\begin{bmatrix}
K_{cc} & K_{ci} \\
K_{ic} & K_{ii}
\end{bmatrix}
\begin{bmatrix}
ue \\
\alpha
\end{bmatrix} =
\begin{bmatrix}
fe \\
0
\end{bmatrix}
$$
詳細な式展開は下記の記事を参照されたい。
https://qiita.com/Altaka4128/items/e0ccbcf7447f3c92219c
ここで上式には、$\alpha$が含まれますが、式変形して削除します。
内部自由度に関する項は、求めることもできるが求める必要がないので、節点自由度の項に代入してしまって、解くべき式を単純化して計算負荷を減らそうというわけです。このアイデアを静的縮約と呼びます。
$$
\begin{align}
K_{cc}ue + K_{ci} \alpha = fe \\
K_{ic}ue + K_{ii}\alpha = 0
\end{align}
$$
これを解いて、$\alpha$を削除すると
$$
\begin{gather}
K_e u_e = f_e \\
K_e = K_{cc} - K_{ci}K_{ii}^{\intercal}K_{cc}
\end{gather}
$$
各剛性マトリクスは数値積分の用いて下記で求まります。
(積分点は、四角形１次要素同じ)
$$
\begin{gather}
K_{cc} = t\int{B_C^{\intercal}DB_C}dA = t\sum_{i=1}^{2}\sum_{j=1}^{2}w_iw_jB_C(\xi,\eta)^{\intercal}DB_C(\xi,\eta)detJ(\xi,\eta) \\
K_{ci} = t\int{B_C^{\intercal}DB_I}dA
= t\sum_{i=1}^{2}\sum_{j=1}^{2}w_iw_jB_C(\xi,\eta)^{\intercal}DB_I(\xi,\eta)detJ(\xi,\eta) \\
K_{ic} = t\int{B_I^{\intercal}DB_C}dA
= t\sum_{i=1}^{2}\sum_{j=1}^{2}w_iw_jB_I(\xi,\eta)^{\intercal}DB_C(\xi,\eta)detJ(\xi,\eta) \\
K_{ii} = t\int{B_I^{\intercal}DB_I}dA
= t\sum_{i=1}^{2}\sum_{j=1}^{2}w_iw_jB_I(\xi,\eta)^{\intercal}DB_I(\xi,\eta)detJ(\xi,\eta)
\end{gather}
$$

これで、非適合要素の剛性マトリクスが求められるようになり、$u_e$が求まります。
また、$\epsilon,\sigma$は下記により求められます。
$$
\begin{gather}
\epsilon = B_cu_e + B_I\alpha=B_Cu_e - B_IK_{ii}^{-1}K_{ic}u_e \\
\sigma = D\epsilon
\end{gather}
$$

実装

ここからは、Grasshopper上で非適合要素を実装していきます。
基本的な構成は前回の記事と同じなので、前回の記事を参照してください。

今回は、ソルバーの部分のコードが一部書き換わるだけになります。

$B,K$マトリクスが書き換わりますが、そのほかの記述はそのままとなります。
まずは、$B_C,B_I$を算出するメソッドを用意します。

calcBc

    public Matrix<double> calcBc(Matrix<double> x,Matrix<double> y, double xi, double eta)
    {
        Matrix<double> dndxi = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+eta)/4,(1-eta)/4,(1+eta)/4,(-1-eta)/4 }
        });
        Matrix<double> dndeta = Matrix<double>.Build.DenseOfArray(new double[,]
        {
                {(-1+xi)/4,(-1-xi)/4,(1+xi)/4,(1-xi)/4 }
        });
        for(int i =0 ; i <4 ; i++)
        {

        }
        Matrix<double> dxdxi = dndxi * x;
        Matrix<double> dydxi = dndxi * y;
        Matrix<double> dxdeta = dndeta * x;
        Matrix<double> dydeta = dndeta * y;
        Matrix<double> J = calcJ(x,y,xi,eta);

        Matrix<double> Bc = Matrix<double>.Build.Dense(3, 8);

        for (int i = 0; i < 4; i++)
        {
            Matrix<double> Bi = J.Solve(Matrix<double>.Build.DenseOfArray(new double[2, 1] { { dndxi[0, i] }, { dndeta[0, i] } }));
            Bc[0, 2 * i] = Bi[0, 0];
            Bc[1, 2 * i + 1] = Bi[1, 0];
            Bc[2, 2 * i] = Bi[1, 0];
            Bc[2, 2 * i + 1] = Bi[0, 0];
        }
        return Bc;
    }

calcBi

    //Biマトリクス(3x4)を計算するメソッド
    public Matrix<double> calcBi(Matrix<double> x,Matrix<double> y, double xi, double eta)
    {
        Matrix<double> J = calcJ(x,y,xi,eta);
        Matrix<double> invJ =J.Inverse();
        Matrix<double> Bi = Matrix<double>.Build.DenseOfArray(new double[3, 4]
        {
            {-2*xi*invJ[0,0],-2*eta*invJ[0,1],0,0},
            {0,0,-2*xi*invJ[1,0],-2*eta*invJ[1,1]},
            {-2*xi*invJ[1,0],-2*eta*invJ[1,1],-2*xi*invJ[0,0],-2*eta*invJ[0,1]}
        });
        return Bi;
    }

$K_e,B_e$を算出するメソッド

calcBe

    private Matrix<double> calcBe(Matrix<double> x,Matrix<double> y,double _xi, double _eta,Matrix<double> D,double Section)
    {
        Matrix<double> Kcc = Matrix<double>.Build.Dense(8, 8);
        Matrix<double> Kci = Matrix<double>.Build.Dense(8, 4);
        Matrix<double> Kic = Matrix<double>.Build.Dense(4, 8);
        Matrix<double> Kii = Matrix<double>.Build.Dense(4, 4);

        for (int i = 0; i < gps.GetLength(0); i++)
        {
            double xi = gps[i, 0];
            double eta = gps[i, 1];
            double wi = gps[i, 2];
            double wj = gps[i, 3];
            Matrix<double> J = calcJ(x,y,xi,eta);
            Matrix<double> Bc = calcBc(x,y,xi,eta);
            Matrix<double> Bi = calcBi(x,y,xi,eta);
            Matrix<double> BcT = Bc.Transpose();
            Matrix<double> BiT = Bi.Transpose();

            Kcc += wi * wj * BcT * D * Bc * J.Determinant() * Section;
            Kci += wi * wj * BcT * D * Bi * J.Determinant() * Section;
            Kic += wi * wj * BiT * D * Bc * J.Determinant() * Section;
            Kii += wi * wj * BiT * D * Bi * J.Determinant() * Section;
        }
        var Be = calcBc(x,y, _xi, _eta) - calcBi(x,y, _xi, _eta) * Kii.Solve(Kic);
        return Be;
    }

calcKe

    public Matrix<double> calcKe(Matrix<double> x,Matrix<double> y,double Section,Matrix<double> D)
    {
        Matrix<double> Kcc = Matrix<double>.Build.Dense(8, 8);
        Matrix<double> Kci = Matrix<double>.Build.Dense(8, 4);
        Matrix<double> Kic = Matrix<double>.Build.Dense(4, 8);
        Matrix<double> Kii = Matrix<double>.Build.Dense(4, 4);

        for (int i = 0; i < gps.GetLength(0); i++)
        {
            double xi = gps[i, 0];
            double eta = gps[i, 1];
            double wi = gps[i, 2];
            double wj = gps[i, 3];
            Matrix<double> J = calcJ(x,y,xi,eta);
            Matrix<double> Bc = calcBc(x,y,xi,eta);
            Matrix<double> Bi = calcBi(x,y,xi,eta);
            Matrix<double> BcT = Bc.Transpose();
            Matrix<double> BiT = Bi.Transpose();

            Kcc += wi * wj * BcT * D * Bc * J.Determinant() * Section;
            Kci += wi * wj * BcT * D * Bi * J.Determinant() * Section;
            Kic += wi * wj * BiT * D * Bc * J.Determinant() * Section;
            Kii += wi * wj * BiT * D * Bi * J.Determinant() * Section;
        }
        var Ke = Kcc - Kci * Kii.Solve(Kic);
        return Ke;
    }

これで非適合要素の要素剛性マトリクスが算出でき、計算できるようになりました。

検証

1次要素との比較

非適合要素の解析結果を、1次要素と比較しながら検証していきます。

1次要素と非適合要素を比べると、メッシュが荒い時はかなり差が出ていて、非適合要素では曲げ変形をよく表現できていることが分かります。
かなりメッシュが荒い状態でも理論値に対して、0.2パーセントほどの誤差しかないのはすごいですね。
変形は精度よく再現できますが、応力はあくまで積分点での応力を外挿して求めているため、メッシュが荒いと理論値との誤差もそこそこある点は、設計者としては把握しておくべき特徴かと思います。

計算時間については1次要素に対して、20%程増加していました。

Kramba3Dとの比較

Grasshopperの構造解析コンポーネント「Karamba3D」の結果とも比較していきましょう。

かなり細かくメッシュを切ったKarambaモデルと、四角形非適合要素の比較しました。
要素数はKaramba：8218,非適合要素：34となっています。
それに対して、解析結果の変位の精度は梁理論による理論値に対して、

• Karmba：97.7％
• 非適合要素：100.5%

とはるかに少ない要素分割で精度も上回る結果となりました。
Karambaでは要素内で「定ひずみ」を仮定する三角形1次要素を採用しており、精度がそこまでよくないみたいですね。

解析時間も見てみましょう。メッシュが荒い分、非適合要素のほうが速度がかなり早くなっています。
三角形1次要素では、積分点が1のため基本的には計算負荷はかなり小さく高速で動作しますが、メッシュを細かく切らないと精度が得られないため、非適合要素同程度の解析結果を得ようとすると、どうしても要素数が増えてしまい、このような結果になってしまいます。

KarambaはRhino+Grasshopperの環境上で実装されていることから、高速にパラスタをすることが目的であるため、解析精度はそこそこに処理速度に主眼を置いて開発されているのでしょう。この特性を理解して使用することが設計者としては大事ですね。

おわりに

本記事では、平面応力要素の四角形1次非適合要素の定式化について説明し、Grasshopper上で実装し、その解析精度について検証しました。最後に、本記事で紹介したポイントをまとめます。

• 非適合要素は、1次要素に比べて荒いメッシュでも変位を精度よく表現できます。
  ただし、節点位置の応力は4つの積分点から補間して求めているため、応力勾配とメッシュの設定には注意が必要です。
• 非適合要素では、Karamba3D（三角形1次要素）と比べて、より高い精度の解析結果が得られます。
  荒いメッシュでもそこそこの解析精度を確保できるため、計算負荷を抑えながら解析することが可能です。
• Karamba3Dは、Grasshopper環境での高速なパラスタを重視して三角形1次要素を採用している（と思います）。こうした特性を把握して、各フェーズで適当なツールを使用していくのが大事ですね。