「Javaによるはじめての有限要素法」を読んでみた

概要

Javaによるはじめての有限要素法(長岐滋 著)
https://www.coronasha.co.jp/np/isbn/9784339046069/

有限要素法の入門に当たり、上記の書籍を参考に勉強してみました。
その内容を当記事でまとめてみます。

特徴

この書籍では1次元、2次元有限要素法を取り扱っており、非常に分かり易く
原理を説明しているという印象を受けました。
また、実際にプログラミングするところまで説明されており、1次元に関しては
ソースコードまで丁寧に解説されています。

内容

第2章

有限要素法に必要な1次元での仮想仕事の原理と材料力学の基礎知識を解説しています。
図は少ないですが1次元なので数式がすっきりしていて分かりやすいです。

第3章 ~ 第4章

Javaに関する基本的な知識について書かれており、Javaを初めて扱う人向けに書かれている印象でした。逆にJavaが分かっている人は読み飛ばしても構わない項目だと思います。(私はほぼ読み飛ばしました)

第5章

行列とベクトルの基礎知識に関して書かれています。基本的な演算のみなので、知識がある方は読み飛ばしても構わない項目です。

第6章

いよいよ、有限要素法の原理の解説に入ります。ここでは、1次元でプログラムを組み、課題を解くところまで取り扱っています。

6章で扱っている単純な引張り変形モデル

形状関数から要素剛性方程式の解説まで多くの図で丁寧に記載されており、
具体例により解説が進んでいくのでイメージが付きやすいです。
また、境界条件をどうプログラミングで扱うか等細かく記載されており、
この内容は重宝しました。

第7章～第8章

ここで2次元の有限要素法の原理について解説になります。
ここでは、原理のみの解説でソースコードの解説は無しでした。
また、形状は三角形のみですが、付録のCDにプログラムがあり、
そちらを使ってCAE解析を体験することが可能です。
実際にメッシュを切るところまで用意されており、その辺りは非常に丁寧な印象を受けました。ただ、2次元の三角形要素に関しては有限要素法のつくり方!(石川博幸、青木伸輔、日比学 著)でより詳しくソースコード付きで書かれているので、そちらの参照をお勧めします。

実践

第6章で作られていたJavaのプログラムを参考にC#(WPF)で実際に作ってみました。
C#で作ったのは可視化がJavaに比べて容易だったためです。
入力にcsvファイルを読み込む形式にしたので、色々な1次元の解析に使えると思います。

ソースコードは下記に置きました。
https://github.com/Altaka4128/Simple1DFEM
使用には入力csvファイルを読み込んで、解析開始ボタンを押すと、
exeの場所にoutput.csvファイルが出力され、その中に結果が書き込まれています。
SampleDataにcsvファイルのサンプルをいくつか置いておきます。(実際のテキストの課題)

また、有限要素法の計算部のソースコードをこちらに記載しておきます。

Beam1DFEM.cs

using MathNet.Numerics.LinearAlgebra.Double;
using System;
using System.Collections.Generic;
using System.Globalization;
using System.IO;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Simple1DFEM
{

    public struct InputData
    {

        public List<double> Node;
        public List<BeamElement> Elem;
    }

    public class Beam1DFEM
    {
        private int NodeNum = 0;   // 節点数
        public List<BeamElement> BeamElems   // 要素の集合
        {
            get;
            private set;
        }
        public DenseVector DispVector   // 変位の境界条件
        {
            get;
            private set;
        }

        public DenseVector ForceVector   // 荷重の境界条件
        {
            get;
            private set;
        }

        public List<bool> Rest   // 拘束の境界条件
        {
            get;
            private set;
        }
        bool AnalysisFlg = false;

        public Beam1DFEM(int nodenum, List<BeamElement> beamelems)
        {
            NodeNum = nodenum;
            BeamElems = beamelems;
        }

        public Beam1DFEM(_1DFEMData data)
        {
            // 要素の形式を変換して格納する
            List<BeamElement> elems = new List<BeamElement>();
            {
                for (int i = 0; i < data.elems.Count; i++)
                {
                    Node[] nodes = new Node[2];
                    nodes[0].No = data.elems[i].NodeNo1;
                    nodes[0].Point = data.nodes[nodes[0].No - 1].Point;
                    nodes[1].No = data.elems[i].NodeNo2;
                    nodes[1].Point = data.nodes[nodes[1].No - 1].Point;

                    int materialNo = data.elems[i].MaterialNo - 1;
                    double area = data.materials[materialNo].Area;
                    double young = data.materials[materialNo].Young;

                    elems.Add(new BeamElement(nodes, area, young));
                }
            }
            NodeNum = data.nodes.Count;
            BeamElems = elems;

            // 拘束条件の形式を変換して格納する
            // 変位
            List<double> disp = new List<double>();
            List<double> force = new List<double>();
            List<bool> constraint = new List<bool>();
            for (int i = 0; i < (data.nodes.Count * 1); i++)
            {
                disp.Add(data.nodes[i].Displacement);
                force.Add(data.nodes[i].Force);
                constraint.Add(data.nodes[i].Constraint);
            }
            DenseVector dispVector = DenseVector.OfArray(disp.ToArray());
            DenseVector forceVector = DenseVector.OfArray(force.ToArray());
            setBoundaryCondition(dispVector, forceVector, constraint);
        }

        // Kマトリックスを作成する
        private DenseMatrix makeKMatrix()
        {
            // 例外処理
            if (Rest == null)
            {
                return null;
            }
            if (NodeNum <= 0 || BeamElems == null || Rest.Count != NodeNum)
            {
                return null;
            }

            DenseMatrix kMatrix = DenseMatrix.Create(NodeNum, NodeNum, 0.0);

            // 各要素のKeマトリックスを計算し、Kマトリックスに統合する
            for (int i = 0; i < BeamElems.Count; i++)
            {
                Console.WriteLine("要素" + (i + 1).ToString());
                DenseMatrix keMatrix = BeamElems[i].makeKeMatrix();

                for (int r = 0; r < 2; r++)
                {
                    int rt = BeamElems[i].Nodes[r].No - 1;
                    for (int c = 0; c < 2; c++)
                    {
                        int ct = BeamElems[i].Nodes[c].No - 1;
                        kMatrix[rt, ct] += keMatrix[r, c];
                    }
                }
            }

            Console.WriteLine("Kマトリックス");
            Console.WriteLine(kMatrix);

            // 境界条件を考慮して修正する
            ForceVector = ForceVector - kMatrix * DispVector;
            for (int i = 0; i < Rest.Count; i++)
            {
                if (Rest[i] == true)
                {
                    for (int j = 0; j < kMatrix.ColumnCount; j++)
                    {
                        kMatrix[i, j] = 0.0;
                    }
                    for (int k = 0; k < kMatrix.RowCount; k++)
                    {
                        kMatrix[k, i] = 0.0;
                    }
                    kMatrix[i, i] = 1.0;

                    ForceVector[i] = DispVector[i];
                }
            }

            Console.WriteLine("Kマトリックス(境界条件考慮)");
            Console.WriteLine(kMatrix);
            Console.WriteLine("荷重ベクトル(境界条件考慮)");
            Console.WriteLine(ForceVector);

            return kMatrix;
        }

        // 境界条件を設定する
        public void setBoundaryCondition(DenseVector dispvector, DenseVector forcevector, List<bool> rest)
        {
            DispVector = dispvector;
            ForceVector = forcevector;
            Rest = rest;
        }

        public void Analysis()
        {
            DenseMatrix kMatrix = makeKMatrix();

            if(kMatrix == null)
            {
                return;
            }

            // 変位を計算する
            DispVector = (DenseVector)(kMatrix.Inverse().Multiply(ForceVector));
            Console.WriteLine("変位ベクトル");
            Console.WriteLine(DispVector);

            // 各要素の応力を計算する
            DenseVector dispElemVector = DenseVector.Create(2, 0.0);
            for (int i = 0; i < BeamElems.Count; i++)
            {
                for (int j = 0; j < 2; j++)
                {
                    dispElemVector[j] = DispVector[BeamElems[i].Nodes[j].No - 1];
                }

                Console.WriteLine("要素" + (i + 1).ToString());
                BeamElems[i].makeStrainVector(dispElemVector);
                BeamElems[i].makeStressVector();
            }

            AnalysisFlg = true;
        }

        // 結果を出力する
        public void outputReport()
        {
            if (AnalysisFlg != true)
            {
                return;
            }

            StreamWriter sw = new StreamWriter("output.csv");

            // 節点変位を書き込む
            List<Node> nodes = new List<Node>();
            for(int i = 0; i < BeamElems.Count; i++)
            {
                for (int j = 0; j < 2; j++)
                {
                    nodes.Add(BeamElems[i].Nodes[j]);
                }
            }
            nodes = nodes.Distinct().ToList();
            nodes.Sort((a,b) => a.No - b.No);

            sw.WriteLine("Node, Coordinate, Ux");
            for (int i = 0; i < NodeNum; i++)
            {
                sw.WriteLine((i + 1).ToString() + ", " + nodes[i].Point + ", " + DispVector[i]);
            }

            // 要素情報を書き込む
            sw.WriteLine("Element, X-Stress");
            for(int i = 0; i < BeamElems.Count; i++)
            {
                sw.WriteLine((i + 1).ToString() + "," + BeamElems[i].StressVector[0]);
            }
            sw.Close();
        }
    }
}

BeamElement.cs

using MathNet.Numerics.LinearAlgebra.Double;
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Simple1DFEM
{
    public class BeamElement
    {
        public Node[] Nodes
        {
            get;
            private set;
        }
        private double Area;
        private double Young;
        private DenseMatrix KeMatrix;
        public DenseVector StrainVector
        {
            get;
            private set;
        }
        public DenseVector StressVector
        {
            get;
            private set;
        }

        public BeamElement()
        {
        }

        public BeamElement(
            Node[] nodes,
            double area,
            double young)
        {

            if (nodes.Length != 2)
            {
                return;
            }

            Nodes = nodes;
            Area = area;
            Young = young;

        }

        // Bマトリックスを計算する
        private DenseMatrix makeBMatirx()
        {
            // 例外処理
            if (Area <= 0)
            {
                return null;
            }

            double length = Nodes[1].Point - Nodes[0].Point;
            double[,] bmatrixArray = new double[1, 2];
            bmatrixArray[0, 0] = -1.0 / length;
            bmatrixArray[0, 1] = 1.0 / length;

            return DenseMatrix.OfArray(bmatrixArray);
        }

        // Keマトリックスを計算する
        public DenseMatrix makeKeMatrix()
        {
            // Bマトリックスを計算する
            DenseMatrix BMatrix = makeBMatirx();
            Console.WriteLine("Bマトリックス");
            Console.WriteLine(BMatrix);

            // 例外処理
            if (BMatrix == null || Young <= 0)
            {
                return null;
            }

            double Volume = Area * (Nodes[1].Point - Nodes[0].Point);
            var keMatrix = Young * Volume * BMatrix.Transpose() * BMatrix;
            DenseMatrix KeMatrix = DenseMatrix.OfColumnArrays(keMatrix.ToColumnArrays());
            Console.WriteLine("Keマトリックス");
            Console.WriteLine(KeMatrix);

            return KeMatrix;
        }

        // ひずみベクトルを計算する
        public void makeStrainVector(DenseVector dispvector)
        {
            DenseMatrix bMatrix = makeBMatirx();
            StrainVector = (DenseVector)bMatrix.Multiply(dispvector);
            Console.WriteLine("ひずみベクトル");
            Console.WriteLine(StrainVector);
        }

        // 応力ベクトルを計算する
        public void makeStressVector()
        {
            if (StrainVector == null || Young <= 0)
            {
                return;
            }

            StressVector = Young * StrainVector;
            Console.WriteLine("応力ベクトル");
            Console.WriteLine(StressVector);
        }

        public BeamElement ShallowCopy()
        {
            return (BeamElement)MemberwiseClone();
        }
    }
}

計算の流れ
要素のデータを入力(座標、ヤング率など)
↓
Dマトリックスを計算する
↓
Bマトリックスを計算する
↓
D、Bマトリックスから各要素のKマトリックス(ソースコード内のKeマトリックス)を計算する
↓
各要素のKマトリックスを計算し、全体のKマトリックスに足し合わせる(重ね合わせの原理)
↓
全体剛性方程式ができるが、Kの逆行列がこのままだと計算できないので、境界条件を加味したKマトリックス及び荷重ベクトルに修正する
↓
剛性方程式から変位を求める

K * \vec{u} = \vec{f} \\
\vec{u} = K^{-1} * \vec{f}

↓
変位から各要素のひずみ、応力が求まる

\vec{\epsilon} = B * \vec{u} \\
\vec{\sigma} = D * \vec{\epsilon}

参考

[1]Javaによるはじめての有限要素法(長岐滋 著)
[2]<解析塾秘伝>有限要素法のつくり方!(石川博幸、青木伸輔、日比学 著)
[3]有限要素法入門(三好俊郎 著)