Gridap.jlで実装する非線形構造解析
Juliaの有限要素法ライブラリ Gridap.jl を使って、
有限変形（非線形）構造解析を実装してみます。
非線形構造解析、理屈を理解するのは本当に難しいですよね。和書でちゃんとしたものあまりない気がします。翻訳書では、「非線形有限要素法のための連続体力学 Javier Bonet (著), Richard D. Wood (著)」はわかりやすかったです。しかし、実装となるとまたハードルが高くなります。
Gridapは非線形方程式最適化を割と扱いやすくしてくれます。
非線形の残差とヤコビアンを記述すれば、Newton法による非線形解法はGridapが自動的に処理してくれるためです。
そのため、複雑なソルバ実装を意識せずに、数式に対応したコードを書くことに集中できます。
この記事では：
- なぜGridap / Juliaが良いのか
- どんな数式を解いているのか
- 実装がどう数式と対応しているのか
をざーっくり書きます。時間もないし、詳細はあまり書かないです。
（なんだか便利そう！が伝わるといいですが・・・
線形問題を解いてみたのはこちら。シンプルに試してみたい方はこちらのほうがいいかも？
https://eye.kohei-kevin.com/2025/08/21/solve-linear-elasticity-problem-with-gridap/
具体的な実装はGitHubにおきました。
https://github.com/kevin-tofu/gridap-test/blob/main/neohookean.jl
なぜGridap / Juliaが良いのか
① 数式に近いコードが書ける
Gridapの最大の特徴はこれです。
例えば弱形式（線形の剛性行列の組み立てですね）：
\int_\Omega \nabla v : \mathbb{C} : \nabla u \, d\Omega
が、そのまま：
∫( inner( ∇(v), C ⊙ ∇(u) ) )dΩ
と書けます。
抽象度高くて、慣れないプログラマーには、何これ、、となるかもしれません。
（PythonのFEMライブラリに比べ、利用可能な演算子が多いので書き味が良いのだと思います。）
② 非線形問題が自然に書ける
通常のFEMコード：
- 残差
- 接線剛性（ヤコビアン）
を別々に書く必要があります。
Gridapでは：
op = FEOperator(res, jac, U, V0)
数学的な「非線形問題そのもの」を記述していけます。
③ Juliaの利点
- 高速（Cに近い）
- 自動微分と相性が良い
- 数値計算に特化
研究用途だけでなく、幅広い用途に使える。
非線形構造解析
非線形PDEも、残差とヤコビアンを記述するだけで比較的容易に解けます。
通常はヤコビアンの導出やNewton法の実装が必要になりますが、Gridapではそれらを意識せずに進めます。
cantilever
- 左端固定
- 右端に表面力
- 材料：Neo-Hookean（非線形弾性）
基本式 と Gridapコードの対応
変形勾配
F = I + \nabla u
F(∇u) = I + ∇u
コーシー グリーンテンソル
C = F^T F
C(Fx) = transpose(Fx) ⋅ Fx
体積変化
J = \det(F) = \sqrt{det(C)}
J(Fx) = sqrt(det(transpose(Fx) ⋅ Fx))
応力（2nd Piola-Kirchhoff）
Neo-Hookean：
S = \mu (I - C^{-1}) + \lambda \log(J) C^{-1}
function S(∇u)
Fx = F(∇u)
Cx = transpose(Fx) ⋅ Fx
Cinv = inv(Cx)
Jx = sqrt ∘ det(Cx)
μ*(I - Cinv) + λ * (log ∘ Jx) * Cinv
end
弱形式
仮想仕事原理：
\int_\Omega \delta E : S \space d\Omega
=\int_{\Gamma} t \cdot v \space d\Gamma
ここで：
\delta E = \frac{1}{2}(\nabla \delta u \cdot F + (\nabla \delta u \cdot F)^T)
res(uh, v) =
∫( dc( dE(∇(v), ∇(uh)), S(∇(uh)) ) )dΩ -
lf[] * ∫( inner(tvec, v) )dΓ_right
接線剛性（ヤコビアン）
非線形問題なのでNewton法を使います：
K = \frac{\partial R}{\partial u}
幾何非線形
幾何学非線形のみ考慮
jac(uh, du, v) = jac_geo(uh, du, v)
jac_geo(uh, du, v) =
∫( dc( ∇(v), S(∇(uh)) ⋅ ∇(du) ) )dΩ
非線形ソルバ
nls = NLSolver(
method = :newton,
linesearch = BackTracking()
)
荷重増分
いきなり全荷重を方程式にぶち込むと、うまく収束しません。
そこで、少しづつ、外力を加えていくという戦略が一般的に取られます。
load_factors = collect(0.05:0.05:1.0)
for f in load_factors
lf[] = f
uh, _ = solve!(uh, solver, op)
end
結果
