はじめに
土は粒子とその間を埋める水もしくは空気で構成されています。よって、数値解析等で土の挙動を追うためには土粒子と水・空気の動きを同時に把握する必要があります。大多数の解析では簡単のため（安全側をみて？）土粒子間は水で飽和していると仮定し、これを土・水連成解析といいます。今回は有限要素法を用いた土・水連成解析をPythonで実装しました。
支配方程式の離散化
有限要素定式化は以下の条件で行います。
・土粒子は線形弾性体とする
・水はダルシー則に従う
・静的問題のみを考える
・微小変形理論に基づく定式化
・二次元平面ひずみ条件
支配方程式は以下に示すように運動量方程式，連続の式，ひずみの適合条件，有効応力の原理，ダルシー則，土骨格の構成式による連立偏微分方程式にまとめられます。
\frac{\partial \sigma_{ji}}{\partial x_i} = 0
\frac{\partial \varepsilon_v}{\partial t} = \frac{\partial v_{ji}}{\partial x_i}
\varepsilon_{ij} = \frac{1}{2} (\frac{\partial u_j}{\partial x_i} + \frac{\partial u_i}{\partial x_j})
\sigma_{ij}' = \sigma_{ij} - \delta_{ij} p_w
v_{ij} = - \frac{k_{ij}}{\gamma_w} \text{grad} p_w
\sigma_{ij} = D_{ijkl}^{e} \varepsilon_{kl}
有効応力は簡単にいうと，土全体が受け持つ応力のうち土粒子が受け持つ応力のことです。$k_{ij}$は土の透水係数（水の通りやすさ）です。
有限要素定式化はZienkiewicz and Bettesの手法を用います。
定式化には四角形一次要素を用います。なお，一般的に変位と間隙水圧を同一の形状関数で補完することは必ずしも論理的に整合性のある方法とは言い難く，多くの場合変位と間隙水圧に異なる形状関数が採用されますが，ここでは簡単のため変位と間隙水圧に同一の形状関数を用いることとします。
任意の要素の節点変位ベクトル$\boldsymbol{u}'$，節点間隙水圧ベクトル$\boldsymbol{p}_w'$を次式で定義します。
\boldsymbol{u}' = [u'_{x1} \quad u'_{y1} \quad u'_{x2} \quad u'_{y2} \quad u'_{x3} \quad u'_{y3} \quad u'_{x4} \quad u'_{y4}]^T
\boldsymbol{p}_w' = [p'_{w1} \quad p'_{w2} \quad p'_{w3} \quad p'_{w4}]^T
変位-ひずみ関係は$[B]$マトリクスを用いて表されます。
\{u\} = [N] \{u'\}
\{\varepsilon\} = [B] \{u'\}
[B] =
\begin{bmatrix}
\displaystyle \frac{\partial N_1}{\partial x} & 0 & \displaystyle \frac{\partial N_2}{\partial x} & 0 & \displaystyle \frac{\partial N_3}{\partial x} & 0 & \displaystyle \frac{\partial N_4}{\partial x} & 0 \\
& \\
0 & \displaystyle \frac{\partial N_1}{\partial y} & 0 &\displaystyle \frac{\partial N_2}{\partial y} & 0 & \displaystyle \frac{\partial N_3}{\partial y} & 0 & \displaystyle \frac{\partial N_4}{\partial y} \\
& \\
\displaystyle \frac{\partial N_1}{\partial x} & \displaystyle \frac{\partial N_2}{\partial y} & \displaystyle \frac{\partial N_3}{\partial x} & \displaystyle \frac{\partial N_3}{\partial y} & \displaystyle \frac{\partial N_4}{\partial x} & \displaystyle \frac{\partial N_4}{\partial y} & \displaystyle \frac{\partial N_4}{\partial x} & \displaystyle \frac{\partial N_4}{\partial y}
\end{bmatrix}
変位-体積ひずみ関係は$[B_{\Delta}]$マトリクスで表されます。
\{\varepsilon_v\} = [B_{\Delta}] \{u'\}
[B_{\Delta}] =
\begin{bmatrix}
\displaystyle \frac{\partial N_1}{\partial x} & \displaystyle \frac{\partial N_1}{\partial y} & \displaystyle \frac{\partial N_2}{\partial x} & \displaystyle \frac{\partial N_2}{\partial y} & \displaystyle \frac{\partial N_3}{\partial x} & \displaystyle \frac{\partial N_3}{\partial y} & \displaystyle \frac{\partial N_4}{\partial x} & \displaystyle \frac{\partial N_4}{\partial y}
\end{bmatrix}
要素内の間隙水圧勾配は$[B_q]$マトリクスで表されます。
\{p_w\} = [B_q] \{p_w'\}
[B_q] =
\begin{bmatrix}
\displaystyle \frac{\partial N_1}{\partial x} & \displaystyle \frac{\partial N_2}{\partial x} & \displaystyle \frac{\partial N_3}{\partial x} & \displaystyle \frac{\partial N_4}{\partial x} \\
& \\
\displaystyle \frac{\partial N_1}{\partial y} & \displaystyle \frac{\partial N_2}{\partial y} & \displaystyle \frac{\partial N_3}{\partial y} & \displaystyle \frac{\partial N_4}{\partial y}
\end{bmatrix}
表面力についても同様に形状関数で補完します．
\{\bar{t}\} = [N] \{\bar{t}'\}
仮想変位ベクトル$\delta \boldsymbol{u}$と仮想ひずみベクトル$\delta \boldsymbol{\varepsilon}$を定義すると，仮想仕事の原理は次式で表されます。
\int_V \{\delta \varepsilon\}^T \{\sigma\} dV - \int_V \{\delta u\}^T \{X\} dV - \int_{S_{1}} \{\delta u\}^T \{\bar{t}\} dS = 0
\{\delta u\} = [\delta u_x \quad \delta u_y]^T
\{\delta \varepsilon\} = [\delta \varepsilon_x \quad \delta \varepsilon_y \quad \delta \tau_{xy}]^T
\{X\} = [0 \quad -\gamma_t]^T
\sum_{m=1}^M \left(\int_{V_m} \{\delta \varepsilon\}^T \{\sigma\} dV - \int_V \{\delta u\}^T \{X\} dV - \int_{S_{1m}} \{\delta u\}^T \{\bar{t}\} dS \right) = 0
\sum_{m=1}^M \left(\int_{V_m} - \{u'\}^T [B]^T \{\sigma\} dV - \int_V \{\delta u'\}^T [N]^T \{X\} dV - \int_{S_{1m}} \{\delta u'\}^T [N]^T \{\bar{t}\} dS \right) = 0
\sum_{m=1}^M \left(\int_{V_m} \left[[B]^T [D] [B] - [B]^T [\sigma'_0] - [B]^T {u_w} - [N]^T \{X\} \right] dV - \int_{S_{1m}} [N]^T \{\bar{t}\} dS \right) = 0
ここで，${u_w}$を間隙水圧を用いて表すと，${u_w} = [1 \quad 1 \quad 0] [N] {p'_w}$と表され，$[B]^T$を乗じると次式で表される．
[B]^T \{u_w\} = [B]^T [1 \quad 1 \quad 0] [N] \{p'_w\} = [B_{\Delta}] [N] \{p'_w\}
解析領域全体の節点変位ベクトルを${u}$，節点間隙水圧ベクトルを${p_w}$とすると，次式でまとめられます。
[K] \{u\} + [C] \{p_w\} = - \{M_1\} + \{M_2\} + \{P_1\}
[K] = \sum_{m=1}^M \int_{V_m} [B]^T [D] [B] dV
[C] = - \sum_{m=1}^M \int_{V_m} [B_{\Delta}]^T [N] dV
\{M_1\} = \sum_{m=1}^M \int_{V_m} [B]^T [\sigma'_0] dV
\{M_2\} = \sum_{m=1}^M \int_{V_m} [N]^T \{X\} dV
\{P_1\} = \sum_{m=1}^M \int_{S_{1m}} [N]^T [N] \{\bar{t}'\} dS
一方，間隙水の流れはダルシー則と連続式より導かれる次式を基礎式とします。
\frac{\partial}{\partial x} \left( \frac{k}{\gamma_w} \frac{\partial p_w}{\partial x} \right) + \frac{\partial}{\partial y} \left\{ \frac{k}{\gamma_w} \left(\frac{\partial p_w}{\partial y} + \gamma_w \right) \right\} - \frac{\partial \varepsilon_v}{\partial t} = 0
上式にGalerkin法（重み付き残差法の一種）を適用します。まず，上式に重み関数$w$を乗じ，領域全体で積分します。
\int_V w \left[ \frac{\partial}{\partial x} \left( \frac{k}{\gamma_w} \frac{\partial p_w}{\partial x} \right) + \frac{\partial}{\partial y} \left\{ \frac{k}{\gamma_w} \left(\frac{\partial p_w}{\partial y} + \gamma_w \right) \right\} - \frac{\partial \varepsilon_v}{\partial t} \right] dV = 0
Galerkin法では，この重み関数を形状関数と一致させます。ガウスの発散定理を適用し，境界条件を考慮すると次式が得られます。
\int_V N_i \left[ \frac{\partial}{\partial x} \left( \frac{k}{\gamma_w} \frac{\partial p_w}{\partial x} \right) + \frac{\partial}{\partial y} \left\{ \frac{k}{\gamma_w} \left(\frac{\partial p_w}{\partial y} + \gamma_w \right) \right\} - \frac{\partial \varepsilon_v}{\partial t} \right] dV = 0
\sum_{m=1}^{M} \left( \int_{V_m} \left[ \frac{\partial N_i}{\partial x} \frac{k}{\gamma_w} \frac{\partial p_w}{\partial x} + \frac{\partial N_i}{\partial y} \frac{k}{\gamma_w} \frac{\partial p_w}{\partial y} + k \frac{\partial N_i}{\partial t} - N_i \frac{\partial \varepsilon_v}{\partial t} \right] dV - \int_{S_{im}} N_i \bar{q} dS \right) = 0
ここで，上式左辺第1項を整理すると次式のようになります。
\begin{split}
& \sum_{m=1}^{M} \int_{V_m} \left( \frac{\partial N_i}{\partial x} \frac{k}{\gamma_w} \frac{\partial p_w}{\partial x} + \frac{\partial N_i}{\partial y} \frac{k}{\gamma_w} \frac{\partial p_w}{\partial y} \right) dV \\
& = \sum_{m=1}^{M} \int_{V_m} \left[\frac{\partial N_i}{\partial x} \quad \frac{\partial N_i}{\partial y} \right]
\begin{bmatrix}
\displaystyle \frac{k}{\gamma_w} & 0 \\
& \\
0 & \displaystyle \frac{k}{\gamma_w}
\end{bmatrix}
\begin{pmatrix}
\displaystyle \frac{\partial p_w}{\partial x} \\
& \\
\displaystyle \frac{\partial p_w}{\partial y}
\end{pmatrix} dV
\end{split}
全ての$i$について書き下すと次式が得られます。
\sum_{m=1}^{M} \int_{V_m} [B_q]^T [k] [B_q] \{p_w\} dV
[k] = \frac{1}{\gamma_w}
\begin{bmatrix}
k & 0 \\
0 & k
\end{bmatrix}
同様に第2項，第3項，第4項の積分を全ての$i$について書き下すと次式が得られます。
\sum_{m=1}^{M} \int_{V_m} [B_q]^T [k] \{X_w\} dV
- \sum_{m=1}^{M} \int_{V_m} ([B_{\Delta}][N])^T \frac{d}{dt} \{u\} dV
\sum_{m=1}^{M} \int_{S_{2m}} [N]^T [N] \{\bar{q}'\} dS
以上を解析領域全体でまとめると次式が得られます。
[C] \frac{d}{dt} \{u\} - [\bar{K}] \{p_w\} = \{M_3\} - \{P_2\}
[\bar{K}] = \sum_{m=1}^{M} \int_{V_m} [B_q]^T [k] [B_q] dV
\{M_3\} = \sum_{m=1}^{M} \int_{V_m} [B_q]^T [k] \{X_w\} dV
\{P_2\} = \sum_{m=1}^{M} \int_{S_{2m}} [N]^T [N] \{\bar{q}'\} dS
後退差分を適用すると，最終的に解くべき多元連立方程式が得られます。
\begin{bmatrix}
[K] & [C] \\
[C]^T & - \Delta t [\bar{K}]
\end{bmatrix}
\begin{pmatrix}
u(t + \Delta t) \\
p_w(t + \Delta t)
\end{pmatrix}
=
\begin{bmatrix}
0 \\
[C]^T
\end{bmatrix}
\begin{pmatrix}
u(t) \\
p_w(t)
\end{pmatrix}
+
\begin{pmatrix}
F \\
R
\end{pmatrix}
F = - \{M_1(t + \Delta t)\} + \{M_2(t + \Delta t)\} + \{P_1(t + \Delta t) \}
R = \Delta t \{M_3(t + \Delta t)\} - \Delta t \{P_2(t + \Delta t)\}
Pythonによる実装
有限要素解析を実施するため，以下のクラスを実装します。
・線形弾性体材料のクラス（elastic_material）
・節点情報を格納するクラス（node）
・要素剛性マトリクスを計算するクラス（element）
・境界条件を考慮するクラス（boundary）
・FEMのメインルーチン（FEM）
線形弾性体材料のクラス
import numpy as np
class ElasticMaterial:
def __init__(self, young_modulus, poissons_ratio, unit_volume_weight, permeability_coefficient):
self.young_modulus = young_modulus # ヤング係数
self.poissons_ratio = poissons_ratio # ポアソン比
self.unit_volume_weight = unit_volume_weight # 単位体積重量（非使用）
self.permeability_coefficient = permeability_coefficient # 透水係数
# Dマトリクスの作成
def makeDmatrix(self):
tmp = self.young_modulus / \
((1.0 - 2.0 * self.poissons_ratio) * (1.0 + self.poissons_ratio))
Dmatrix = np.matrix([[1.0 - self.poissons_ratio, self.poissons_ratio, 0.0], [self.poissons_ratio,
1.0 - self.poissons_ratio, 0.0], [0.0, 0.0, 0.5 * (1.0 - 2.0 * self.poissons_ratio)]])
Dmatrix *= tmp
return Dmatrix
節点情報を格納するクラス
class Node:
def __init__(self, number, x, y):
self.number = number
self.x = x
self.y = y
要素剛性マトリクスを作成するクラス
from elastic_material import ElasticMaterial
import numpy as np
import numpy.linalg as LA
NODE_DOF = 3
NODE_NUM = 4
class Element:
def __init__(self, ele_id, nodes, material_id):
self.ele_id = ele_id
self.nodes = nodes
self.material_id = material_id
self.ai = [-np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0),
np.sqrt(1.0 / 3.0), -np.sqrt(1.0 / 3.0)]
self.bi = [-np.sqrt(1.0 / 3.0), -np.sqrt(1.0 / 3.0),
np.sqrt(1.0 / 3.0), np.sqrt(1.0 / 3.0)]
self.Jmatrix = []
for i in range(0, NODE_NUM):
self.Jmatrix.append(self.makeJmatrix(self.ai[i], self.bi[i]))
self.Bmatrix = []
for i in range(0, NODE_NUM):
self.Bmatrix.append(self.makeBmatrix(
self.ai[i], self.bi[i], self.Jmatrix[i]))
self.Brematrix = []
for i in range(0, NODE_NUM):
self.Brematrix.append(self.makeBrematrix(
self.ai[i], self.bi[i], self.Jmatrix[i]))
self.Bbarmatrix = []
for i in range(0, NODE_NUM):
self.Bbarmatrix.append(self.makeBbarmatrix(
self.ai[i], self.bi[i], self.Jmatrix[i]))
self.Nvector = []
for i in range(0, NODE_NUM):
self.Nvector.append(self.makeNvector(self.ai[i], self.bi[i]))
# 要素剛性マトリクスの作成
def makeKematrix(self, young_modulus, poissons_ratio, unit_volume_weight, permeability_coefficient, delta_time):
Kematrix = np.zeros(
(NODE_DOF * NODE_NUM, NODE_DOF * NODE_NUM), dtype=np.float64)
Kematrix_c = np.zeros(
(NODE_DOF * NODE_NUM, NODE_DOF * NODE_NUM), dtype=np.float64)
Cmatrix = self.makeCmatrix()
Kbarmatrix = self.makeKbarmatrix(permeability_coefficient)
elastic = ElasticMaterial(
young_modulus, poissons_ratio, unit_volume_weight, permeability_coefficient)
Dmatrix = elastic.makeDmatrix()
for i in range(0, NODE_NUM):
Kematrix += self.Bmatrix[i].T * Dmatrix * \
self.Bmatrix[i] * LA.det(self.Jmatrix[i])
for i in range(0, 8):
for j in range(0, 4):
Kematrix[i, j + 8] += Cmatrix[i, j]
Kematrix[j + 8, i] += Cmatrix[i, j]
for i in range(0, 4):
for j in range(0, 4):
Kematrix[i + 8, j + 8] -= delta_time * Kbarmatrix[i, j]
for i in range(0, NODE_NUM):
for j in range(0, NODE_DOF * NODE_NUM):
Kematrix_c[NODE_DOF * i, j] = Kematrix[2 * i, j]
Kematrix_c[NODE_DOF * i + 1, j] = Kematrix[2 * i + 1, j]
Kematrix_c[NODE_DOF * i + 2, j] = Kematrix[i + 8, j]
for i in range(0, NODE_DOF * NODE_NUM):
for j in range(0, NODE_NUM):
Kematrix[i, NODE_DOF * j] = Kematrix_c[i, 2 * j]
Kematrix[i, NODE_DOF * j + 1] = Kematrix_c[i, 2 * j + 1]
Kematrix[i, NODE_DOF * j + 2] = Kematrix_c[i, j + 8]
return Kematrix
# Cマトリクスの作成
def makeCmatrix(self):
Cmatrix = np.zeros((8, 4), dtype=np.float64)
identify_tensor = np.matrix([[1.0], [1.0], [0.0]])
for i in range(0, NODE_NUM):
Cmatrix -= self.Brematrix[i].T * identify_tensor * \
self.Nvector[i] * LA.det(self.Jmatrix[i])
return Cmatrix
# Kbarマトリクスの作成
def makeKbarmatrix(self, permeability_coefficient):
Kbarmatrix = np.zeros((4, 4), dtype=np.float64)
tmp = np.matrix([[permeability_coefficient, 0.0],
[0.0, permeability_coefficient]])
for i in range(0, NODE_NUM):
Kbarmatrix += self.Bbarmatrix[i].T * \
(tmp / 9.8) * self.Bbarmatrix[i] * LA.det(self.Jmatrix[i])
return Kbarmatrix
# Bマトリクスの作成
def makeBmatrix(self, ai, bi, Jmatrix):
dNdab = self.makedNdab(ai, bi)
dNdxy = LA.solve(Jmatrix, dNdab)
Bmatrix = np.matrix([[dNdxy[0, 0], 0.0, dNdxy[0, 1], 0.0, dNdxy[0, 2], 0.0, dNdxy[0, 3], 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, dNdxy[1, 0], 0.0, dNdxy[1, 1], 0.0,
dNdxy[1, 2], 0.0, dNdxy[1, 3], 0.0, 0.0, 0.0, 0.0], [dNdxy[1, 0], dNdxy[0, 0], dNdxy[1, 1], dNdxy[0, 1], dNdxy[1, 2], dNdxy[0, 2], dNdxy[1, 3], dNdxy[0, 3], 0.0, 0.0, 0.0, 0.0]])
return Bmatrix
# ひずみ算出用Bマトリクスの作成
def makeBrematrix(self, ai, bi, Jmatrix):
dNdab = self.makedNdab(ai, bi)
dNdxy = LA.solve(Jmatrix, dNdab)
Brematrix = np.matrix([[dNdxy[0, 0], 0.0, dNdxy[0, 1], 0.0, dNdxy[0, 2], 0.0, dNdxy[0, 3], 0.0], [0.0, dNdxy[1, 0], 0.0, dNdxy[1, 1], 0.0, dNdxy[1, 2], 0.0, dNdxy[1, 3]], [
dNdxy[1, 0], dNdxy[0, 0], dNdxy[1, 1], dNdxy[0, 1], dNdxy[1, 2], dNdxy[0, 2], dNdxy[1, 3], dNdxy[0, 3]]])
return Brematrix
# Kbarマトリクス用Bマトリクスの作成
def makeBbarmatrix(self, ai, bi, Jmatrix):
dNdab = self.makedNdab(ai, bi)
dNdxy = LA.solve(Jmatrix, dNdab)
Bbarmatrix = np.matrix([[dNdxy[0, 0], dNdxy[0, 1], dNdxy[0, 2], dNdxy[0, 3]], [
dNdxy[1, 0], dNdxy[1, 1], dNdxy[1, 2], dNdxy[1, 3]]])
return Bbarmatrix
# ヤコビマトリクスの作成
def makeJmatrix(self, ai, bi):
dNdab = self.makedNdab(ai, bi)
matxiyi = np.matrix([[self.nodes[0].x, self.nodes[0].y], [self.nodes[1].x, self.nodes[1].y], [
self.nodes[2].x, self.nodes[2].y], [self.nodes[3].x, self.nodes[3].y]])
Jmatrix = dNdab * matxiyi
if LA.det(Jmatrix) < 0.0:
raise ValueError('要素の計算に失敗しました')
return Jmatrix
# dNdab行列の計算
def makedNdab(self, ai, bi):
dN1da = -0.25 * (1.0 - bi)
dN2da = 0.25 * (1.0 - bi)
dN3da = 0.25 * (1.0 + bi)
dN4da = -0.25 * (1.0 + bi)
dN1db = -0.25 * (1.0 - ai)
dN2db = -0.25 * (1.0 + ai)
dN3db = 0.25 * (1.0 + ai)
dN4db = 0.25 * (1.0 - ai)
dNdab = np.matrix([[dN1da, dN2da, dN3da, dN4da],
[dN1db, dN2db, dN3db, dN4db]])
return dNdab
# Nベクトルの計算
def makeNvector(self, ai, bi):
n1 = 0.25 * (1.0 - ai) * (1.0 - bi)
n2 = 0.25 * (1.0 + ai) * (1.0 - bi)
n3 = 0.25 * (1.0 + ai) * (1.0 + bi)
n4 = 0.25 * (1.0 - ai) * (1.0 + bi)
Nvector = np.matrix([n1, n2, n3, n4])
return Nvector
# ひずみ、応力の算出
def calcurate_stress(self, disp_tmp, pwp_tmp, young_modulus, poissons_ratio, unit_volume_weight, permeability_coefficient):
strain_x = np.zeros(4, dtype=np.float64)
strain_y = np.zeros(4, dtype=np.float64)
strain_xy = np.zeros(4, dtype=np.float64)
stress_x = np.zeros(4, dtype=np.float64)
stress_y = np.zeros(4, dtype=np.float64)
stress_xy = np.zeros(4, dtype=np.float64)
pore_water_pressure = np.zeros(4, dtype=np.float64)
elastic = ElasticMaterial(
young_modulus, poissons_ratio, unit_volume_weight, permeability_coefficient)
Dmatrix = elastic.makeDmatrix()
for ii in range(0, NODE_NUM):
strain = np.array(self.Brematrix[ii] @ disp_tmp).flatten()
strain_x[ii] = strain[0]
strain_y[ii] = strain[1]
strain_xy[ii] = strain[2]
stress = np.array(Dmatrix @ strain).flatten()
stress_x[ii] = stress[0]
stress_y[ii] = stress[1]
stress_xy[ii] = stress[2]
pwp = np.array(self.Nvector[ii] @ pwp_tmp).flatten()
pore_water_pressure[ii] = pwp[0]
return strain_x, strain_y, strain_xy, stress_x, stress_y, stress_xy, pore_water_pressure
境界条件を考慮するクラス
import numpy as np
NODE_DOF = 3
NODE_NUM = 4
class Boundary:
def __init__(self, disp_node_no, disp_node_xy, force_node_xy):
self.disp_node_no = disp_node_no
self.disp_node_xy = disp_node_xy
self.force_node_xy = force_node_xy
# 荷重ベクトルの作成
def makeForcevector(self, step, nodes, elements, materials, node_pore_water_pressure, time_step, delta_time):
Fvector = np.zeros(NODE_DOF * len(nodes), dtype=np.float64)
# 上載荷重の入力
for ii in range(0, NODE_DOF * len(nodes)):
Fvector[ii] = self.force_node_xy[step][ii] / time_step
# 過剰間隙水圧の入力
for ele in elements:
m = ele.material_id - 1
permeability_coefficient = materials[m].permeability_coefficient
Kbarmatrix = ele.makeKbarmatrix(permeability_coefficient)
i = ele.nodes[0].number - 1
j = ele.nodes[1].number - 1
k = ele.nodes[2].number - 1
l = ele.nodes[3].number - 1
tmp = np.array([node_pore_water_pressure[i], node_pore_water_pressure[j],
node_pore_water_pressure[k], node_pore_water_pressure[l]])
pwp = delta_time * np.array(Kbarmatrix @ tmp).flatten()
Fvector[i * NODE_DOF + 2] += pwp[0]
Fvector[j * NODE_DOF + 2] += pwp[1]
Fvector[k * NODE_DOF + 2] += pwp[2]
Fvector[l * NODE_DOF + 2] += pwp[3]
return Fvector
# 剛性マトリクス、荷重ベクトルへ境界条件の入力
def setBoundarycondition(self, step, Kmatrix, Fvector, nodes):
Kmatrix_c = np.copy(Kmatrix)
Fvector_c = np.copy(Fvector)
for ii in range(0, len(nodes)):
for jj in range(0, NODE_DOF):
iz = ii * NODE_DOF + jj
if self.disp_node_no[step][iz] == 1:
vecx = np.zeros(NODE_DOF * len(nodes), dtype=np.float64)
vecx = Kmatrix_c[:, iz]
Fvector_c -= self.disp_node_xy[step][iz] * vecx
Kmatrix_c[:, iz] = 0.0
Kmatrix_c[iz, :] = 0.0
Kmatrix_c[iz, iz] = 1.0
return Kmatrix_c, Fvector_c
# 変位ベクトルへ境界条件の入力
def setdispBoundarycondition(self, step, disp, nodes):
disp_c = np.copy(disp)
for ii in range(0, len(nodes)):
for jj in range(0, NODE_DOF):
iz = ii * NODE_DOF + jj
if self.disp_node_no[step][iz] == 1:
disp_c[iz] = self.disp_node_xy[step][iz]
return disp_c
FEMのメインルーチン
from boundary import Boundary
import numpy as np
import numpy.linalg as LA
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import os
NODE_DOF = 3
NODE_NUM = 4
class FEM:
def __init__(self, nodes, elements, materials, bounds):
self.nodes = nodes
self.elements = elements
self.materials = materials
self.bounds = bounds
def analysis(self, number_of_steps, total_time, delta_time):
step = 0
analysis_step = 0
for ii in range(0, number_of_steps):
step += int(total_time[ii] / delta_time[ii])
displacement = np.zeros(
(NODE_DOF - 1) * len(self.nodes), dtype=np.float64)
disp_all = np.zeros(
(step, (NODE_DOF - 1) * len(self.nodes)), dtype=np.float64)
strain_x = np.zeros((4, len(self.elements)), dtype=np.float64)
strain_y = np.zeros((4, len(self.elements)), dtype=np.float64)
strain_xy = np.zeros((4, len(self.elements)), dtype=np.float64)
strain_x_all = np.zeros((step, len(self.elements)), dtype=np.float64)
strain_y_all = np.zeros((step, len(self.elements)), dtype=np.float64)
strain_xy_all = np.zeros((step, len(self.elements)), dtype=np.float64)
stress_x = np.zeros((4, len(self.elements)), dtype=np.float64)
stress_y = np.zeros((4, len(self.elements)), dtype=np.float64)
stress_xy = np.zeros((4, len(self.elements)), dtype=np.float64)
stress_x_all = np.zeros((step, len(self.elements)), dtype=np.float64)
stress_y_all = np.zeros((step, len(self.elements)), dtype=np.float64)
stress_xy_all = np.zeros((step, len(self.elements)), dtype=np.float64)
node_pore_water_pressure = np.zeros(len(self.nodes), dtype=np.float64)
pore_water_pressure = np.zeros(
(4, len(self.elements)), dtype=np.float64)
pore_water_pressure_all = np.zeros(
(step, len(self.elements)), dtype=np.float64)
analysis_time = np.zeros(step, dtype=np.float64)
# 時刻歴計算開始
for ii in range(0, number_of_steps):
time_step = int(total_time[ii] / delta_time[ii])
for _jj in range(0, time_step):
analysis_time[analysis_step] += delta_time[ii]
# 全体剛性マトリクスの作成
Kmatrix = self.makeKmatrix(delta_time[ii])
# 荷重ベクトルの作成
Fvector = self.bounds[ii].makeForcevector(
ii, self.nodes, self.elements, self.materials, node_pore_water_pressure, time_step, delta_time[ii])
# 境界条件の考慮
Kmatrix_c, Fvector_c = self.bounds[ii].setBoundarycondition(
ii, Kmatrix, Fvector, self.nodes)
# 変位ベクトルの計算
Kmatrix_c = csr_matrix(Kmatrix_c)
disp = spsolve(Kmatrix_c, Fvector_c, use_umfpack=True)
# 境界条件の考慮
disp_c = self.bounds[ii].setdispBoundarycondition(
ii, disp, self.nodes)
# ひずみ、応力の算出
for ele in self.elements:
disp_tmp = np.zeros(8, dtype=np.float64)
pwp_tmp = np.zeros(4, dtype=np.float64)
i = ele.nodes[0].number - 1
j = ele.nodes[1].number - 1
k = ele.nodes[2].number - 1
l = ele.nodes[3].number - 1
disp_tmp[0] = disp_c[NODE_DOF * i]
disp_tmp[1] = disp_c[NODE_DOF * i + 1]
disp_tmp[2] = disp_c[NODE_DOF * j]
disp_tmp[3] = disp_c[NODE_DOF * j + 1]
disp_tmp[4] = disp_c[NODE_DOF * k]
disp_tmp[5] = disp_c[NODE_DOF * k + 1]
disp_tmp[6] = disp_c[NODE_DOF * l]
disp_tmp[7] = disp_c[NODE_DOF * l + 1]
pwp_tmp[0] = disp_c[NODE_DOF * i + 2]
pwp_tmp[1] = disp_c[NODE_DOF * j + 2]
pwp_tmp[2] = disp_c[NODE_DOF * k + 2]
pwp_tmp[3] = disp_c[NODE_DOF * l + 2]
m = ele.material_id - 1
young_modulus = self.materials[m].young_modulus
poissons_ratio = self.materials[m].poissons_ratio
unit_volume_weight = self.materials[m].unit_volume_weight
permeability_coefficient = self.materials[m].permeability_coefficient
strain_x_tri, strain_y_tri, strain_xy_tri, stress_x_tri, stress_y_tri, stress_xy_tri, pwp_tri = ele.calcurate_stress(
disp_tmp, pwp_tmp, young_modulus, poissons_ratio, unit_volume_weight, permeability_coefficient)
for ip in range(0, NODE_NUM):
strain_x[ip, ele.ele_id - 1] += strain_x_tri[ip]
strain_y[ip, ele.ele_id - 1] += strain_y_tri[ip]
strain_xy[ip, ele.ele_id - 1] += strain_xy_tri[ip]
stress_x[ip, ele.ele_id - 1] += stress_x_tri[ip]
stress_y[ip, ele.ele_id - 1] += stress_y_tri[ip]
stress_xy[ip, ele.ele_id - 1] += stress_xy_tri[ip]
pore_water_pressure[ip, ele.ele_id - 1] += pwp_tri[ip]
for node in range(0, len(self.nodes)):
displacement[(NODE_DOF - 1) *
node] += disp_c[NODE_DOF * node]
displacement[(NODE_DOF - 1) * node +
1] += disp_c[NODE_DOF * node + 1]
node_pore_water_pressure[node] += disp_c[NODE_DOF * node + 2]
for node in range(0, len(self.nodes)):
disp_all[(analysis_step, (NODE_DOF - 1) * node)
] = displacement[(NODE_DOF - 1) * node]
disp_all[(analysis_step, (NODE_DOF - 1) * node + 1)
] = displacement[(NODE_DOF - 1) * node + 1]
strain = np.zeros((3, len(self.elements)), dtype=np.float64)
stress = np.zeros((3, len(self.elements)), dtype=np.float64)
element_pore_water_pressure = np.zeros(
len(self.elements), dtype=np.float64)
for ele in self.elements:
for ip in range(0, NODE_NUM):
strain[0, ele.ele_id - 1] += strain_x[ip, ele.ele_id - 1]
strain[1, ele.ele_id - 1] += strain_y[ip, ele.ele_id - 1]
strain[2, ele.ele_id -
1] += strain_xy[ip, ele.ele_id - 1]
stress[0, ele.ele_id - 1] += stress_x[ip, ele.ele_id - 1]
stress[1, ele.ele_id - 1] += stress_y[ip, ele.ele_id - 1]
stress[2, ele.ele_id -
1] += stress_xy[ip, ele.ele_id - 1]
element_pore_water_pressure[ele.ele_id -
1] += pore_water_pressure[ip, ele.ele_id - 1]
strain[0, ele.ele_id - 1] *= 0.25
strain[1, ele.ele_id - 1] *= 0.25
strain[2, ele.ele_id - 1] *= 0.25
stress[0, ele.ele_id - 1] *= 0.25
stress[1, ele.ele_id - 1] *= 0.25
stress[2, ele.ele_id - 1] *= 0.25
element_pore_water_pressure[ele.ele_id - 1] *= 0.25
strain_x_all[analysis_step, ele.ele_id -
1] = strain[0, ele.ele_id - 1]
strain_y_all[analysis_step, ele.ele_id -
1] = strain[1, ele.ele_id - 1]
strain_xy_all[analysis_step, ele.ele_id -
1] = strain[2, ele.ele_id - 1]
stress_x_all[analysis_step, ele.ele_id -
1] = stress[0, ele.ele_id - 1]
stress_y_all[analysis_step, ele.ele_id -
1] = stress[1, ele.ele_id - 1]
stress_xy_all[analysis_step, ele.ele_id -
1] = stress[2, ele.ele_id - 1]
pore_water_pressure_all[analysis_step, ele.ele_id -
1] = element_pore_water_pressure[ele.ele_id - 1]
analysis_step += 1
# 計算結果の出力
self.output(analysis_step, analysis_time, disp_all, strain_x_all, strain_y_all, strain_xy_all,
stress_x_all, stress_y_all, stress_xy_all, pore_water_pressure_all)
# 全体剛性マトリクスの作成
def makeKmatrix(self, delta_time):
Kmatrix = np.zeros((NODE_DOF * len(self.nodes),
NODE_DOF * len(self.nodes)), dtype=np.float64)
for ele in self.elements:
m = ele.material_id - 1
young_modulus = self.materials[m].young_modulus
poissons_ratio = self.materials[m].poissons_ratio
unit_volume_weight = self.materials[m].unit_volume_weight
permeability_coefficient = self.materials[m].permeability_coefficient
Kematrix = ele.makeKematrix(
young_modulus, poissons_ratio, unit_volume_weight, permeability_coefficient, delta_time)
for c in range(len(ele.nodes) * NODE_DOF):
ct = (ele.nodes[c // NODE_DOF].number - 1) * \
NODE_DOF + c % NODE_DOF
for r in range(len(ele.nodes) * NODE_DOF):
rt = (ele.nodes[r // NODE_DOF].number -
1) * NODE_DOF + r % NODE_DOF
Kmatrix[ct, rt] += Kematrix[c, r]
return Kmatrix
# 計算結果の出力
def output(self, analysis_step, time_data, displacement, strain_x, strain_y, strain_xy, stress_x, stress_y, stress_xy, pore_water_pressure):
if not os.path.exists("output"):
# ディレクトリが存在しない場合、ディレクトリを作成する
os.makedirs("output")
node_file = open("output/nodes.txt", "w")
element_file = open("output/elements.txt", "w")
columNum = 12
floatDigits = ".5g"
node_file.write("***** Node Data ******\n")
node_file.write("No".rjust(columNum) + "X".rjust(columNum) +
"Y".rjust(columNum) + "\n")
node_file.write("-" * columNum * 3 + "\n")
for node in self.nodes:
strNo = str(node.number).rjust(columNum)
strX = str(format(node.x, floatDigits).rjust(columNum))
strY = str(format(node.y, floatDigits).rjust(columNum))
node_file.write(strNo + strX + strY + "\n")
node_file.write("\n")
element_file.write("***** Element Data ******\n")
element_file.write("No".rjust(columNum) +
"Node No.1".rjust(columNum) +
"Node No.2".rjust(columNum) +
"Node No.3".rjust(columNum) +
"Node No.4".rjust(columNum) + "\n")
element_file.write("-" * columNum * 6 + "\n")
for elem in self.elements:
strNo = str(elem.ele_id).rjust(columNum)
strNodeNo = ""
for node in elem.nodes:
strNodeNo += str(node.number).rjust(columNum)
strNodeNo = strNodeNo.rjust(columNum)
element_file.write(strNo + strNodeNo + "\n")
element_file.write("\n")
for step in range(analysis_step):
node_file.write(f"***** Step {step} *****\n")
node_file.write(
"Time: " + str(format(time_data[step], floatDigits)) + "\n")
node_file.write("Disp. X, Disp. Y:\n")
for i in range(len(self.nodes)):
strXdisp = str(
format(displacement[step, (NODE_DOF - 1) * i], floatDigits)).rjust(columNum)
strYdisp = str(
format(displacement[step, (NODE_DOF - 1) * i + 1], floatDigits)).rjust(columNum)
node_file.write(strXdisp + strYdisp + "\n")
node_file.write("\n")
element_file.write(f"***** Step {step} *****\n")
element_file.write(
"Time: " + str(format(time_data[step], floatDigits)) + "\n")
element_file.write(
"Stress X, Stress Y, Stress XY, Strain X, Strain Y, Strain XY, PWP:\n")
for i in range(len(self.elements)):
element_file.write(str(format(stress_x[step, i], floatDigits)).rjust(columNum) + str(format(stress_y[step, i], floatDigits)).rjust(columNum) + str(format(stress_xy[step, i], floatDigits)).rjust(columNum) + str(format(
strain_x[step, i], floatDigits)).rjust(columNum) + str(format(strain_y[step, i], floatDigits)).rjust(columNum) + str(format(strain_xy[step, i], floatDigits)).rjust(columNum) + str(format(pore_water_pressure[step, i], floatDigits)).rjust(columNum) + '\n')
element_file.write("\n\n")
node_file.close()
element_file.close()
解析データのインプットの例
from node import Node
from element import Element
from elastic_material import ElasticMaterial
from boundary import Boundary
from fem import FEM
import numpy as np
NODE_DOF = 3
def main():
# 解析情報の読み込み
input_file = 'input/input.txt'
file = open(input_file, 'r')
text = file.readline()
text = text.strip()
text = text.split()
number_of_nodes = int(text[0])
number_of_elements = int(text[1])
number_of_materials = int(text[2])
number_of_steps = int(text[3])
nodes = []
# 節点情報の読み込み
for i in range(0, number_of_nodes):
text = file.readline()
text = text.strip()
text = text.split()
node = Node(int(text[0]), float(text[1]), float(text[2]))
nodes.append(node)
elements = []
# 要素情報の読み込み
for i in range(0, number_of_elements):
text = file.readline()
text = text.strip()
text = text.split()
element_id = int(text[0])
elem = [nodes[int(text[1]) - 1], nodes[int(text[2]) - 1],
nodes[int(text[3]) - 1], nodes[int(text[4]) - 1]]
material_id = int(text[5])
element = Element(element_id, elem, material_id)
elements.append(element)
materials = []
# 材料定数の読み込み
for i in range(0, number_of_materials):
text = file.readline()
text = text.strip()
text = text.split()
mate = ElasticMaterial(float(text[1]), float(
text[2]), float(text[3]), float(text[4]))
materials.append(mate)
total_time = []
delta_time = []
disp_node_no = []
disp_node_xy = []
force_node_xy = []
bounds = []
# 解析ステップの読み込み
for i in range(0, number_of_steps):
text = file.readline()
text = text.strip()
text = text.split()
total_time.append(float(text[1]))
delta_time.append(float(text[2]))
number_of_pfix = int(text[3])
number_of_force = int(text[4])
disp_no = np.zeros(NODE_DOF * number_of_nodes, dtype=np.int64)
disp_xy = np.zeros(NODE_DOF * number_of_nodes, dtype=np.float64)
force_no = np.zeros(NODE_DOF * number_of_nodes, dtype=np.int64)
force_xy = np.zeros(NODE_DOF * number_of_nodes, dtype=np.float64)
# ノイマン境界条件の読み込み
if number_of_pfix > 0:
for _j in range(0, number_of_pfix):
text = file.readline()
text = text.strip()
text = text.split()
tmp = int(text[0])
disp_no[NODE_DOF * tmp - 3] = int(text[1])
disp_no[NODE_DOF * tmp - 2] = int(text[2])
disp_no[NODE_DOF * tmp - 1] = int(text[3])
disp_xy[NODE_DOF * tmp - 3] = float(text[4])
disp_xy[NODE_DOF * tmp - 2] = float(text[5])
disp_xy[NODE_DOF * tmp - 1] = float(text[6])
# ディリクレ境界条件の読み込み
if number_of_force > 0:
for j in range(0, number_of_force):
text = file.readline()
text = text.strip()
text = text.split()
tmp = int(text[0])
force_no[j] = tmp
force_xy[NODE_DOF * tmp - 3] = float(text[1])
force_xy[NODE_DOF * tmp - 2] = float(text[2])
disp_node_no.append(disp_no)
disp_node_xy.append(disp_xy)
force_node_xy.append(force_xy)
bound = Boundary(disp_node_no, disp_node_xy, force_node_xy)
bounds.append(bound)
file.close()
# 解析の実行
fem = FEM(nodes, elements, materials, bounds)
fem.analysis(number_of_steps, total_time, delta_time)
if __name__ == '__main__':
main()
参考文献
Zienkiewicz，O.C. and Bettes，P.：Soils and other saturated media under transient, dynamic conditions; General formulation and the validity of various simplifying assumptions，Soil mechanics – Transient and Cyclic loads，1982.
地盤工学会編：地盤技術者のためのFEMシリーズ3 有限要素法をつかう，丸善出版，2020.