LoginSignup
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@udi_

FEMの基礎からPINNsまで①

Last updated at Posted at 2025-09-09

はじめに

本記事では1次元の棒の引っ張りの問題を理論解・有限要素法・PINNsでそれぞれ解いて出力を比較します。
簡単な問題設定で有限要素もPINNsの門を叩こうという記事です。
PINNsについての説明は過去の記事や他のページを参考にしてください。
間違い等あればご指摘ください。

想定読者

  • FEMを基礎から学びたい人
  • PINNsやってみたい人

有限要素法

対象とする問題

前述の通り１次元の棒の引っ張りの問題を考えます。
図のように左端が固定されていて、右端に荷重（集中荷重）がかかるような状態です。
それぞれの記号と単位は下記の通りです。
（値がおかしいですが、簡単のためにこのように設定しています。）

  • 長さ：$L=3\ \text{[m]}$
  • ヤング率：$E=2\ \text[{N/m}^2]$
  • 断面積：$A=1\ \text[{m}^2]$
  • 荷重：$P=10\ \text{[N]}$
  • 変位：$u(L) $ ※ これを求める

image.png

材料力学の基礎

ここで材料力学の基礎を押さえておきます。
力と変位や応力とひずみの関係はフックの法則の 剛性 × 変位で表されます。
応力 $σ$ - ひずみ $ε$ ：

$$\sigma = E \varepsilon$$

ひずみ $ε$は変位 / 元の長さで表されます。

$$\quad \varepsilon = \frac{u(L)}{L}$$

また、応力$σ$は荷重 / 面積です。

$$\quad \sigma = \frac{P}{A}$$

後で使うので、荷重を求めるような式(1)と変位を求める式(2)を作っておきます。
$$
\frac{P}{A} = E \cdot \frac{u(L)}{L} ・・・(1)
$$

$$
u(L) = \frac{PL}{EA} ・・・(2)
$$

理論解

まずは理論解を求めてしまいましょう。
式(2)に代入するだけです。15 mも伸びるのはおかしいですが、今回は設定上こうなってしまいます。有限要素法やPINNsの結果とあえばいいのでこのまま進めます。

$$
u(L) = \frac{10 \cdot 3}{2 \cdot 1} = 15 \ \text{[m]}.
$$

有限要素法の基礎

1次元の有限要素法は対象の物体を点（節点）と線（要素）に分割して、それぞれの要素で計算して合わせるようなイメージです。今回は簡単のために節点を3つ、要素を2つで考えてみます。
要素 $a$と要素 $b$はそれぞれが剛性 ${k}^e$を持っています（今回は同じ剛性）。また、今回は1次元を考えるので節点1、2、3の自由度は1で水平方向にしか動きません。
例えば、節点1だけが節点2の方向にxだけ動いたとすると、要素aは$f = {k}^ex$ の力を出します。
このようにそれぞれの節点の動きとその時のそれぞの要素が出す力を考えていくのが有限要素法です。

image.png

有限要素法で求める

支配方程式

支配方程式を確認しておきます。棒の変形の支配方程式は

$$
\frac{d}{dx}\left( EA \frac{du}{dx} \right) + f(x) = 0.
$$

有限要素法では形状関数 $N(x)$ というものを用いて

$$
u(x) \approx N_1(x) u_1 + N_2(x) u_2,
$$

と表します。要素が足しあわされているのがわかるかと思います。

解いていく

要素$a$に注目してみていきます。節点$i$の変位を$u_i$、要素$a$に発生する力を$f^a_i$とします。
節点2で発生する力は、
$$
f^a_1 = k^e(u_1-u_2)
$$
同様に節点2で発生する力は
$$
f^a_2 = k^e(u_2-u_1)
$$

節点は3つあるので$u_3$もありますが、要素$a$に影響がないので$u_3$がかかる部分はすべて0になります。
また、節点1、2の変異によって要素$a$から発生る力の向きは逆であると考えると（ここら辺はうまく理解できているか微妙です。）下記のように書けます。

\begin{bmatrix}
f^a_1 \\
f^a_2 \\
f^a_3
\end{bmatrix}
=
k^e
\begin{bmatrix}
1 & -1 & 0 \\
-1 & 1 & 0 \\
0 & 0 & 0 \\
\end{bmatrix}
\begin{bmatrix}
u_1 \\
u_2 \\
u_3
\end{bmatrix}

要素$b$に関しても同様に求めます。要素$b$には節点1の影響はないので、$u_1$がかかる部分はすべて0です。

\begin{bmatrix}
f^b_1 \\
f^b_2 \\
f^b_3
\end{bmatrix}
=
k^e
\begin{bmatrix}
0 & 0 & 0 \\
0 & 1 & -1 \\
0 & -1 & 1 \\
\end{bmatrix}
\begin{bmatrix}
u_1 \\
u_2 \\
u_3
\end{bmatrix}

これらの$k^e$[行列]を**要素剛性マトリクスと呼びます。
要素合成マトリクスをすべて足し合わせたものが全体剛性マトリクスになります。
有限要素法ではこのマトリクスを考え、ここに境界条件や荷重条件を入れていきます。

\begin{bmatrix}
f_1 \\
f_2 \\
f_3
\end{bmatrix}
=
k^e
\begin{bmatrix}
1 & -1 & 0 \\
1 & 2 & -1 \\
0 & -1 & 1 \\
\end{bmatrix}
\begin{bmatrix}
u_1 \\
u_2 \\
u_3
\end{bmatrix}

変位を求める

$k^e$を求めておきます。式(1)を変形すると、下記の式のようになり、力=剛性×変位の形になります。（$L^e$は要素の長さ）
$$
P = \frac {EA} {L^e}\cdot{u(L)}
$$
$$
k^e = \frac {EA} {L^e}
$$

今回の問題設定の境界条件は、
・左端固定 ： $u_1 = 0$
・右端に集中荷重 ： $[f_1, f_2, f_3] = [0,0,P]$
となります。
これらをさっきの式に入れていきます。

\begin{bmatrix}
0 \\
0 \\
P
\end{bmatrix}
=
\frac {EA} {L^e}
\begin{bmatrix}
-1 & 1 & 0 \\
1 & 2 & -1 \\
0 & -1 & 1 \\
\end{bmatrix}
\begin{bmatrix}
0 \\
u_2 \\
u_3
\end{bmatrix}

あとは数字を入れて代入するだけです。
$u_3 = 15$となり理論解と一致します。

PINNsで解く

わざわざPINNsを使って解いてみます。
問題設定と節点と要素に分けた図を思い出します。
境界条件は

  • 左端固定 : 節点1の変位は0 ($u_1 = 0$)
  • 右端に集中荷重 : $EA \frac{du}{dx}=P$
    偏微分方程式として
  • 節点2には荷重がかからない : $EA \frac{d^2u}{dx^2}=0$ (PINNsは節点の概念がないので仮想断面に荷重がないの方が正確かもしれません）

これらを損失関数として扱ってニューラルネットワークで解いていきます。

image.png

まず、節点を作ります。問題設定では3点でしたが、3点だと結果が合わなかったので50点作っています。
x_pdeで節点をx_bc_pで右端を作っています。
左端固定はニューラルネットワークの構造を作る部分で実装します。

n_pde = 50
x_pde = torch.linspace(0, L, n_pde).view(-1, 1).to(device) # L = 3
x_bc_p = torch.tensor([[L]], dtype=torch.float32).to(device)

損失関数の部分のコードはこのようになります。
上側で$EA \frac{d^2u}{dx^2}=0$を
下側で$EA \frac{du}{dx}=P$を計算しています。

def loss_func(x_pde, x_bc_p):
    # PDE残差 (EA * u'' = 0)
    x_pde.requires_grad = True
    u = model(x_pde)
    du_dx = torch.autograd.grad(u, x_pde, torch.ones_like(u), create_graph=True)[0]
    d2u_dx2 = torch.autograd.grad(du_dx, x_pde, torch.ones_like(du_dx), create_graph=True)[0]
    pde_loss = torch.mean((EA * d2u_dx2)**2)

    # Neumann BC (x=L で EA*u'(L)=P)
    x_bc_p.requires_grad = True
    u_bc_p = model(x_bc_p)
    du_dx_bc_p = torch.autograd.grad(u_bc_p, x_bc_p, torch.ones_like(u_bc_p), create_graph=True)[0]
    bc_p_loss = torch.mean((EA * du_dx_bc_p - P)**2)

    return pde_loss, bc_p_loss

ニューラルネットワークの出力の部分で$x * self.net(x)$とすることで$x = 0$での変位$u_1=0$を担保しています。

class PINN(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(1, 20),
            nn.Tanh(),
            nn.Linear(20, 20),
            nn.Tanh(),
            nn.Linear(20, 20),
            nn.Tanh(),
            nn.Linear(20, 1)
        )

    def forward(self, x):
        # u(0)=0 を保証
        return x * self.net(x)

結果

境界条件の損失（bc_p_loss）と偏微分方程式の損失（pde_loss）に重みをかけています。ここら辺はチューニングが必要かと思います。（うまくチューニングすればepoch数すくなくてもいける・・・？）
結果は3mの部分で15mのびるという理論と一致する解が得られました。

全コード 
import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt

# ===============================
# 物理パラメータ
# ===============================
E = 2   # ヤング率 [Pa]
A = 1    # 断面積 [m^2]
P = 10    # 集中荷重 [N]
L = 3.0     # 棒の長さ [m]
EA = E * A

device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')

# ===============================
# PINN モデル定義
# ===============================
class PINN(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(1, 20),
            nn.Tanh(),
            nn.Linear(20, 20),
            nn.Tanh(),
            nn.Linear(20, 20),
            nn.Tanh(),
            nn.Linear(20, 1)
        )

    def forward(self, x):
        # Dirichlet条件 u(0)=0 を保証
        return x * self.net(x)

# ===============================
# モデル・最適化
# ===============================
model = PINN().to(device)

# ===============================
# 損失関数
# ===============================
def loss_func(x_pde, x_bc_p):
    # PDE残差 (EA * u'' = 0)
    x_pde.requires_grad = True
    u = model(x_pde)
    du_dx = torch.autograd.grad(u, x_pde, torch.ones_like(u), create_graph=True)[0]
    d2u_dx2 = torch.autograd.grad(du_dx, x_pde, torch.ones_like(du_dx), create_graph=True)[0]
    pde_loss = torch.mean((EA * d2u_dx2)**2)

    # Neumann BC (x=L で EA*u'(L)=P)
    x_bc_p.requires_grad = True
    u_bc_p = model(x_bc_p)
    du_dx_bc_p = torch.autograd.grad(u_bc_p, x_bc_p, torch.ones_like(u_bc_p), create_graph=True)[0]
    bc_p_loss = torch.mean((EA * du_dx_bc_p - P)**2)

    return pde_loss, bc_p_loss

# ===============================
# 学習データ
# ===============================
n_pde = 50
x_pde = torch.linspace(0, L, n_pde).view(-1, 1).to(device)
x_bc_p = torch.tensor([[L]], dtype=torch.float32).to(device)

# ===============================
# 学習 
# ===============================
# Step 1: Adam
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
epochs = 5000
for epoch in range(epochs):
    optimizer.zero_grad()
    pde_loss, bc_p_loss = loss_func(x_pde, x_bc_p)
    loss = pde_loss + 10 * bc_p_loss
    loss.backward()
    optimizer.step()

    if (epoch + 1) % 1000 == 0:
        print(f"[Adam] Epoch {epoch+1}/{epochs}, Loss={loss.item():.4e}, PDE={pde_loss.item():.4e}, BCp={bc_p_loss.item():.4e}")

output.png

最後に

1次元の棒の引っ張りという簡単な問題を題材に理論～有限要素法～PINNsの実装まで実施しました。損失の境界条件をいれる問題だったので勉強になりました。（まだまだ理解が甘い部分だらけですが。）
できれば、3次元のものまでやりたい気持ちだけはあります。

参考文献

1
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?