はじめに
Javaと材料力学の復習をかねて作成した。
原理
たわみの基礎式の離散化
分布荷重$p(x)$が加わるヤング率$E$,断面二次モーメント$I$の梁のたわみ$y$は以下の式から求まる。
\frac{d^4y}{dx^4} = \frac{p(x)}{EI}
ここで4階の数値微分を以下の式で計算する。
\frac{d^4y}{dx^4} \sim \frac{y_{i-2}-4y_{i-1}+6y_{i}-4y_{i+1}+y_{i+2}}{{\Delta x}^4}
よってたわみの基礎式の離散化は以下の式となる。
\frac{y_{i-2}-4y_{i-1}+6y_{i}-4y_{i+1}+y_{i+2}}{{\Delta x}^4} = \frac{p(x_i)}{EI}
\frac{1}{{\Delta x}^4}
\begin{bmatrix}
1 & -4 & 6 & -4 & 1\\
\end{bmatrix}
\begin{bmatrix}
y_{i-2} \\
y_{i-1} \\
y_{i} \\
y_{i+1} \\
y_{i+2} \\
\end{bmatrix}
= \frac{p(x_i)}{EI}
節点が$N$個のとき、$i=2$から$i=N-2$までを行列にまとめると以下の式になる。
（※端点は微分の計算で$y_{-2}$や$y_{N+1}$が出てきてしまうので計算できない。）
\frac{1}{{\Delta x}^4}
\begin{bmatrix}
1 & -4 & 6 & -4 & 1 & 0 & &\cdots && 0\\
0 & 1 & -4 & 6 & -4 & 1 & 0 & \cdots && 0\\
&&&&\vdots &&&&&\\
0 & \cdots && 0 &1 & -4 & 6 & -4 & 1 & 0 \\
0 && \cdots && 0 &1 & -4 & 6 & -4 & 1 \\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
=
\frac{1}{EI}
\begin{bmatrix}
p(x_{1}) \\
p(x_{2}) \\
\vdots \\
p(x_{N-1}) \\
p(x_{N}) \\
\end{bmatrix}
境界条件の離散化
前提として自由端、単純支持、固定端の境界条件は以下の通りとなる。
| 条件 | たわみ$$y$$ | たわみ角$$y'$$ | 曲げモーメント$$EIy''$$ | せん断応力$$EIy'''$$ |
|---|---|---|---|---|
| 自由端 | * | * | 0 | 0 |
| 単純支持 | 0 | * | 0 | * |
| 固定端 | 0 | 0 | * | * |
節点$i=0$では前進差分、節点$i=N$では後退差分で離散化する。
以下に節点$i=0$のときの離散化を示す。
節点$i=N$では左辺の行列の非0成分を右端の列にスライドさせればよい。
y(0)=0 \Rightarrow \begin{bmatrix}
1 & 0 & \cdots & 0\\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
= 0
y(0)'=0 \Rightarrow \frac{y_{2}-y_1}{\Delta x} = 0\Rightarrow \begin{bmatrix}
1 & -1 & 0 & \cdots & 0\\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
= 0
EIy(0)''=0 \Rightarrow \frac{y_{3}-2y_{2}+y_1}{\Delta x ^2} = 0\Rightarrow \begin{bmatrix}
1 & -2 & 1 & 0 & \cdots & 0\\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
= 0
EIy(0)'''=0 \Rightarrow \frac{y_{4}-3y_{3}+3y_{2}-y_1}{\Delta x ^3} = 0\Rightarrow \begin{bmatrix}
1 & -3 & 3 & -1 & 0 & \cdots & 0\\
\end{bmatrix}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
= 0
たわみの基礎式+境界条件の行列計算
ここまでの式をまとめて以下の式で表す。
\boldsymbol{A}
\begin{bmatrix}
y_{1} \\
y_{2} \\
\vdots \\
y_{N-1} \\
y_{N} \\
\end{bmatrix}
= \boldsymbol{b}
この式をLU分解を用いた手法で計算したわみ$y_1,y_2,\cdots,y_{N}$ を求める。
コード
package MechDesign.StructureStudy.Beam;
import MechDesign.StructureStudy.StructureStudy;
import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
import org.apache.commons.math4.legacy.linear.ArrayRealVector;
import org.apache.commons.math4.legacy.linear.DecompositionSolver;
import org.apache.commons.math4.legacy.linear.LUDecomposition;
import org.apache.commons.math4.legacy.linear.RealMatrix;
import org.apache.commons.math4.legacy.linear.RealVector;
public class Beam {
double length;
int nodeNum;
double I;
double E;
char startBoundaryCondition;
char endBoundaryCondition;
public static final class Constraint {
public static final char FIX = 'F';
public static final char PIN = 'P';
public static final char FREE = 'e';
}
public Test(double I, double E, double length, int nodeNum, char startBoundaryCondition,
char endBoundaryCondition) {
this.I = I;
this.E = E;
this.length = length;
this.nodeNum = nodeNum;
this.startBoundaryCondition = startBoundaryCondition;
this.endBoundaryCondition = endBoundaryCondition;
}
public double[] solve(double q) {
// パラメータの設定
double dx = length / (nodeNum - 1); // メッシュ間隔
// 係数行列と右辺ベクトルの設定
RealMatrix A = new Array2DRowRealMatrix(nodeNum, nodeNum);
RealVector b = new ArrayRealVector(nodeNum);
// 内部点の設定
for (int i = 2; i < nodeNum - 2; i++) {
A.setEntry(i, i - 2, 1 / Math.pow(dx, 4));
A.setEntry(i, i - 1, -4 / Math.pow(dx, 4));
A.setEntry(i, i, 6 / Math.pow(dx, 4));
A.setEntry(i, i + 1, -4 / Math.pow(dx, 4));
A.setEntry(i, i + 2, 1 / Math.pow(dx, 4));
b.setEntry(i, q / (E * I));
}
// 端点の境界条件
switch (startBoundaryCondition) {
case Constraint.FIX:
A.setSubMatrix(new double[][] { { 1 } }, 0, 0);
b.setEntry(0, 0.0);
A.setSubMatrix(new double[][] { { 1 } }, 1, 1);
b.setEntry(1, 0.0);
break;
case Constraint.PIN:
A.setSubMatrix(new double[][] { { 1 } }, 0, 0);
b.setEntry(0, 0.0);
A.setSubMatrix(new double[][] { { 1, -2, 1 } }, 1, 0);
b.setEntry(1, 0.0);
break;
case Constraint.FREE:
A.setSubMatrix(new double[][] { { 1, -2, 1 } }, 0, 0);
b.setEntry(0, 0.0);
A.setSubMatrix(new double[][] { { 1, -3, 3, -1 } }, 1, 0);
b.setEntry(1, 0.0);
break;
}
switch (endBoundaryCondition) {
case Constraint.FIX:
A.setSubMatrix(new double[][] { { 1 } }, nodeNum - 2, nodeNum - 1);
b.setEntry(nodeNum - 2, 0.0);
A.setSubMatrix(new double[][] { { 1 } }, nodeNum - 1, nodeNum - 2);
b.setEntry(nodeNum - 1, 0.0);
break;
case Constraint.PIN:
A.setSubMatrix(new double[][] { { 1 } }, nodeNum - 2, nodeNum - 1);
b.setEntry(nodeNum - 2, 0.0);
A.setSubMatrix(new double[][] { { 1, -2, 1 } }, nodeNum - 1, nodeNum - 3);
b.setEntry(nodeNum - 1, 0.0);
break;
case Constraint.FREE:
A.setSubMatrix(new double[][] { { 1, -2, 1 } }, nodeNum - 2, nodeNum - 3);
b.setEntry(nodeNum - 2, 0.0);
A.setSubMatrix(new double[][] { { 1, -3, 3, -1 } }, nodeNum - 1, nodeNum - 4);
b.setEntry(nodeNum - 1, 0.0);
break;
}
// 連立方程式を解く
DecompositionSolver solver = new LUDecomposition(A).getSolver();
RealVector solution = solver.solve(b);
return solution.toArray();
}
public static void main(String[] args) {
Beam beam = new Beam(1e-4, 10., 50, Beam.Constraint.FIX, Beam.Constraint.FIX, 210e-9);
double[] yArr = beam.solve(1000);
double max = yArr[0];
for (var y : yArr) {
System.out.println(y * 1e3);
if (y > max)
max = y;
}
System.out.println("max:" + max);
}