材端回転ばね剛性を考慮した梁要素剛性マトリクスの導出

はじめに

構造解析における材端の結合条件について、自作の解析ソルバーに組み込もうとした際に
参考になる記事があまり見当たらなかったので、備忘録も兼ねて記事にします。

参考

定式化

部材端部に軸方向ばねを挿入する方法は参考のページに記載があります。
本記事では、平面フレームの１節点３自由度（並進２成分、回転１成分）の条件で、梁の両端部に回転ばねを挿入した場合の剛性マトリクスの処理方法について記述します。

まず下図のような、モデルを想定します。

i端j端を両端に持つ梁用素(i,jは３自由度)の両端にそれぞれ節点s,t(s,tは回転１成分)を用意しi-s,t-jに回転ばねを挿入し、４節点８自由度のモデルを考えます

u=[ui,​​vi,​​θi​,​uj,​​vj,​​θj,​​θs,​​θt​​]^{\intercal}

この４節点８自由度における剛性マトリクスは下記のようになります。
端部回転ばねの回転剛性はそれぞれ$k_{rs},k_{rt}$とする

\mathbf{K}
=
\begin{bmatrix}
 \dfrac{EA}{L} & 0 & 0 & -\dfrac{EA}{L} & 0 & 0 & 0 & 0 \\[6pt]
 0 & \dfrac{12EI}{L^3} & \dfrac{6EI}{L^2} & 0 & -\dfrac{12EI}{L^3} & \dfrac{6EI}{L^2} & 0 & 0 \\[6pt]
 0 & \dfrac{6EI}{L^2} & \dfrac{4EI}{L}+k_{rs} & 0 & -\dfrac{6EI}{L^2} & \dfrac{2EI}{L} & -k_{rs} & 0 \\[6pt]
 -\dfrac{EA}{L} & 0 & 0 & \dfrac{EA}{L} & 0 & 0 & 0 & 0 \\[6pt]
 0 & -\dfrac{12EI}{L^3} & -\dfrac{6EI}{L^2} & 0 & \dfrac{12EI}{L^3} & -\dfrac{6EI}{L^2} & 0 & 0 \\[6pt]
 0 & \dfrac{6EI}{L^2} & \dfrac{2EI}{L} & 0 & -\dfrac{6EI}{L^2} & \dfrac{4EI}{L}+k_{rt} & 0 & -k_{rt} \\[6pt]
 0 & 0 & -k_{rs} & 0 & 0 & 0 & k_{rs} & 0 \\[6pt]
 0 & 0 & 0 & 0 & 0 & -k_{rt} & 0 & k_{rt}
\end{bmatrix}

端部の回転ばねを持つ要素だけ自由度を増やすのは、計算するのに都合がよくないので、
内部自由度となるs,tに作用する外力が０であることから、この回転自由度に対応する部分をマトリクスから除去し、6ｘ６のマトリクスに縮小することを考えます。
この処理を静的縮約といい、平面応力要素の非適合要素などでも似た処理を行っています。気になる方は、以下の記事をご覧になってください。

まずは、先ほどの$Ｋ$マトリクスを以下の４マトリクスに分解します。

\mathbf{K}
=
\begin{bmatrix}
k_{rr} & k_{ri} \\
k_{ir} & k_{ii}
\end{bmatrix}

\mathbf{K}_{ii}
=
\begin{bmatrix}
\dfrac{4EI}{L} + k_{rs} & \dfrac{2EI}{L} \\[6pt]
\dfrac{2EI}{L} & \dfrac{4EI}{L} + k_{rt}
\end{bmatrix}

\mathbf{K}_{ir}
=
\mathbf{K}_{ri}^{\intercal}
=
\begin{bmatrix}
0 & \dfrac{6EI}{L^2} & -k_{rs} & 0 & \dfrac{6EI}{L^2} & 0 \\[6pt]
0 & \dfrac{6EI}{L^2} & 0 & 0 & \dfrac{6EI}{L^2} & -k_{rt}
\end{bmatrix}

\mathbf{K}_{rr}
=
\begin{bmatrix}
\dfrac{EA}{L} & 0 & 0 & -\dfrac{EA}{L} & 0 & 0 \\[6pt]
0 & \dfrac{12EI}{L^3} & 0 & 0 & -\dfrac{12EI}{L^3} & 0 \\[6pt]
0 & 0 & k_{rs} & 0 & 0 & 0 \\[6pt]
-\dfrac{EA}{L} & 0 & 0 & \dfrac{EA}{L} & 0 & 0 \\[6pt]
0 & -\dfrac{12EI}{L^3} & 0 & 0 & \dfrac{12EI}{L^3} & 0 \\[6pt]
0 & 0 & 0 & 0 & 0 & k_{rt}
\end{bmatrix}

つり合い式は以下のようになる。

\begin{bmatrix}
K_{rr}&_​K_{ir}\\
​​K_{ri}&​K_{ii}​​
\end{bmatrix}
\begin{bmatrix}
u_r\\
​​u_i
\end{bmatrix}
=
\begin{bmatrix}
f_r\\
f_i
\end{bmatrix}

この時、
$$
\begin{gather}
u_r=[ui,​​vi,​θi​,​uj,​​vj,​​θj​​]^{\intercal} \\
u_i=[θs​,​​θt​​]^{\intercal}
\end{gather}
$$

ここで内部節点に外力は作用しないので　$f_i = 0$ より

K_{ri} u_r + K_{ii} u_i = 0

よって

u_i = -K_{ii}^{-1}K_{ir}​u_r

これを上段の式に代入すると

K_{eq}u_r=f_r

ここで

K_{eq} = K_{rr}-K_{ri}K_{ii}^{-1}K_{ir}

この$K_{eq}$を等価な要素剛性マトリクスとすれば、めでたく材端ばねを考慮した応力計算を行うことができます。

実装

現在私が開発しているGrasshopperの構造解析コンポーネント「FEMur」の中で実際に書いたプログラムの一部を抜粋して記載します。誰かの実装の参考になれば幸いです。
要素剛性マトリクスに対して、Jointの条件を適用するメソッドのを記載しています。
これまでの記事では、1節点３自由度の２次元平面で考えましたが、本プログラムは３次元なので１節点６自由度となっていますが、やっていることは同じです。

        private Matrix<double> ApplyJointStaticCondensation(Matrix<double> keLocal, Joint joint, ElementBase element)
        {
            // 両端とも剛接の場合は縮約不要
            if (joint.IsStartRigid() && joint.IsEndRigid())
            {
                return keLocal;
            }

            // 内部自由度のインデックスを特定
            var internalDofs = new List<int>();
            var retainedDofs = new List<int>();

            // Start端（i端）の自由度チェック: 0-5
            for (int i = 0; i < Joint.DOFS_PER_END; i++)
            {
                if (!joint.IsFixed[i])
                {
                    internalDofs.Add(i);
                }
                else
                {
                    retainedDofs.Add(i);
                }
            }

            // End端（j端）の自由度チェック: 6-11
            for (int i = Joint.DOFS_PER_END; i < Joint.TOTAL_DOFS; i++)
            {
                if (!joint.IsFixed[i])
                {
                    internalDofs.Add(i);
                }
                else
                {
                    retainedDofs.Add(i);
                }
            }

            if (internalDofs.Count == 0)
            {
                return keLocal;
            }

            // 拡張剛性マトリクスを構築（材端ばねを追加）
            // 元の12x12マトリクスに内部自由度分を追加
            int originalSize = 12;
            int internalCount = internalDofs.Count;
            int expandedSize = originalSize + internalCount;

            var kExpanded = Matrix<double>.Build.Dense(expandedSize, expandedSize);

            // 元の剛性マトリクスをコピー
            for (int i = 0; i < originalSize; i++)
            {
                for (int j = 0; j < originalSize; j++)
                {
                    int i2 = i;
                    int j2 = j;

                    if(internalDofs.Contains(i))
                    {
                        i2 = originalSize + internalDofs.IndexOf(i);
                    }
                    if(internalDofs.Contains(j))
                    {
                        j2 = originalSize + internalDofs.IndexOf(j);
                    }

                    kExpanded[i2, j2] = keLocal[i, j];
                }
            }

            // 材端ばねの剛性を追加
            for (int idx = 0; idx < internalCount; idx++)
            {
                int dofIdx = internalDofs[idx];
                double springStiffness = GetSpringStiffness(joint, dofIdx);

                if (springStiffness <= 0)
                {
                    springStiffness = 1.0E-5; // 数値安定性のための微小値
                }

                // 対角成分にばね剛性を加算
                kExpanded[dofIdx, dofIdx] += springStiffness;

                // 内部自由度の対角成分
                int internalIdx = originalSize + idx;
                kExpanded[internalIdx, internalIdx] += springStiffness;

                // 結合項
                kExpanded[dofIdx, internalIdx] = -springStiffness;
                kExpanded[internalIdx, dofIdx] = -springStiffness;
            }

            // K を [K_rr, K_ri; K_ir, K_ii] に分割
            var kRR = kExpanded.SubMatrix(0, originalSize, 0, originalSize);
            var kRI = kExpanded.SubMatrix(0, originalSize, originalSize, internalCount);
            var kIR = kExpanded.SubMatrix(originalSize, internalCount, 0, originalSize);
            var kII = kExpanded.SubMatrix(originalSize, internalCount, originalSize, internalCount);

            // K_eq = K_rr - K_ri * K_ii^(-1) * K_ir
            Matrix<double> kIIinv;
            try
            {
                kIIinv = kII.Inverse();
            }
            catch
            {
                // 逆行列が計算できない場合は警告を出して元のマトリクスを返す
                Warnings.Add($"Joint condensation failed for Element {element.Id}: K_ii is singular. Using original stiffness.");
                return keLocal;
            }

            var keq = kRR - kRI * kIIinv * kIR;

            return keq;
        }

するとこんな感じで無事材端のばねを考慮した応力解析ができるようになりました！

机上の計算を実装して、正しく計算結果が描画されるのはうれしいですね。

ちなみにこのGrasshopperの構造解析は近日中にβ版を公開できると思いますので、公開されたらぜひ触ってみてください。