問題設定

重力を受ける1次元棒の有限変形域における準静的問題を考える。 入力パラメータは以下の通りである。なお、今回はポアソン比 $\nu=0$ として1次元問題を考える。

支配方程式

変形勾配

基準配置を $\bm{X}$、変形後の現在配置を $\bm{x}$ とし、1次元棒の軸方向変位を $\bm{u}(\bm{X}) = \bm{x}(\bm{X}) - \bm{X}$ とする。変形勾配テンソルは $\bm{F} = \partial \bm{x} / \partial \bm{X}$ で定義される。
今、1次元問題を考えるため、$x$方向にのみ変位が生じるとすると、変形勾配テンソルの表現行列は以下のように表される。 $$ [\bm{F}] = \begin{bmatrix} 1 + \frac{\partial u}{\partial X} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$ 変形勾配テンソルを左極分解すると $\bm{F} = \bm{VR}$ となるが、いま回転テンソル $\bm{R}$ は単位テンソル $\bm{I}$ となるため、左ストレッチテンソルは変形勾配テンソルと一致する（ $\bm{V} = \bm{F}$ ）。
つまり、1次元問題における左ストレッチテンソルの表現行列は以下のように表される。 $$ [\bm{V}] = \begin{bmatrix} 1 + \frac{\partial u}{\partial X} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} \lambda & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} $$ ここで、$\lambda = 1 + \frac{\partial u}{\partial X} = \frac{\partial x}{\partial X}$ は軸方向の伸び（主ストレッチ）を表す。

構成則

応力ひずみ関係として Hencky 超弾性構成則を用いる。 ひずみ測度として、左ストレッチテンソルの対数をとった $\bm{\varepsilon} = \ln\bm{V}$ を用いる。 対数ひずみテンソルの表現行列は以下のように表される。 $$ [\bm{\varepsilon}] = \begin{bmatrix} \ln \lambda & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$ Hencky 超弾性構成則に基づく Kirchhoff 応力テンソル $\bm{\tau}$ は以下のように表される。 $$ \bm{\tau} = \left(K+\frac23 \mu \right) \text{tr}(\bm{\varepsilon}) \bm{I} + 2 \mu \bm{\varepsilon} $$ ここで、$K$ は体積弾性係数、$G$ はせん断剛性係数である。ポアソン比 $\nu = 0$ の場合、これらの弾性係数は以下のように表される。 $$ \lambda = \frac{E\nu}{(1+\nu)(1-2\nu)} = 0, \quad \mu = \frac{E}{2(1+\nu)} = \frac{E}{2} $$ よって、Kirchhoff 応力テンソルは $\bm{\tau} = E \bm{\varepsilon}$ と表され、表現行列は次式で与えられる。 $$ [\bm{\tau}] = \begin{bmatrix} E \ln \lambda & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$

強形式

基準配置における準静的問題の力のつり合い式は次式で与えられる。 $$ \frac{\partial}{\partial X_{J}} \left( P_{iJ} \right) + \rho_{0} g_{i} = 0 $$ ここで、$P_{iJ}$ は第一 Piola–Kirchhoff 応力テンソルの成分、$\rho_{0}$ は基準配置における密度、$g_{i}$ は重力加速度ベクトルの成分である。 第一 Piola–Kirchhoff 応力テンソルと Kirchhoff 応力テンソルとの間には $\bm{P} = \bm{\tau} \bm{F}^{-\rm T}$ が成立するので、第一 Piola–Kirchhoff 応力テンソル の表現行列は次式で与えられる。 $$ [P] = [\bm{\tau}] [\bm{F}]^{-1} = \begin{bmatrix} \frac{E \ln \lambda}{\lambda} & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} $$ したがって、$i=1$（$x$方向）のみを考慮する1次元問題に帰着できるので、力のつり合い式は以下のように簡略化される。 $$ \frac{d}{d X} \left( \frac{E \ln \lambda}{\lambda} \right) + \rho_{0} g = 0 $$ ここで、$g$ は $x$ 方向の重力加速度成分である。
境界条件は以下の通りである。 $$ \begin{align*} &u(X=0) = 0 & & \text{(左端固定)} \\ &\lambda(X=L) = 1 & & \text{(右端自由)} \end{align*} $$ 右端自由の条件は、第一 Piola–Kirchhoff 応力テンソルの $x$ 方向成分がゼロであること、すなわち $(E \ln \lambda)/\lambda = 0$ から導かれる。

理論解の導出

力のつり合い式を初期位置 $X$ に関して積分する。 $$ \begin{align*} &\int \frac{d}{d X} \left( \frac{E \ln \lambda}{\lambda} \right) \mathrm{d}X + \int \rho_{0} g \mathrm{d}X = 0 \\ &\rightarrow \frac{E \ln \lambda}{\lambda} + \rho_{0} g X = C \end{align*} $$ 境界条件 $\lambda(X=L) = 1$ を用いると、積分定数 $C$ は $C = \rho_{0} g L$ と求まる。 したがって、以下の非線形方程式を主ストレッチに関して解くことで、任意の位置 $X$ における伸び $\lambda$ を求めることができる。 $$ \begin{align} \frac{\ln \lambda}{\lambda} = \frac{\rho_{0} g}{E} (L - X) \end{align} $$ ちなみにこの問題においては、$(E\ln\lambda)/\lambda$ は 第一 Piola–Kirchhoff 応力テンソルの $P_{11}$ 成分および Cauchy 応力テンソルの $\sigma_{11}$ 成分に相当する。

現在位置 $x$ は、伸び $\lambda = \partial x / \partial X$ を $X$ で積分することで求められる。 $$ x(X) = \int \lambda(X) \mathrm{d}X $$ 右辺の積分は $X$ を $\lambda$ の関数に変換して計算する。 $$ \begin{align*} &X = L - \frac{E}{\rho_{0} g} \frac{\ln \lambda}{\lambda} \\ &\rightarrow \mathrm{d}X = - \frac{E}{\rho_{0} g} \frac{1 - \ln \lambda}{\lambda^{2}} \mathrm{d}\lambda \end{align*} $$ したがって、現在位置 $x$ は以下のように表される。 $$ \begin{align*} x(\lambda) &= -\frac{E}{\rho_{0} g} \int \frac{1 - \ln \lambda}{\lambda} \mathrm{d}\lambda \\ &= -\frac{E}{\rho_{0} g} \left( \ln \lambda - \frac{(\ln \lambda)^{2}}{2} \right) + C' \end{align*} $$ 境界条件 $x(X=0) = 0$ を用いると、積分定数 $C’$ は $C’ = \frac{E}{\rho_{0} g} \left( \ln \lambda(0) - \frac{(\ln \lambda(0))^{2}}{2} \right)$ と求まる。 ここで、$\lambda(0)$ は $X=0$ における伸びであり、非線形方程式 $\frac{\ln \lambda(0)}{\lambda(0)} = \frac{\rho_{0} g}{E} L$ を解くことで求められる。 以上より、現在位置 $x(X)$ は以下のように表される。 $$ \begin{align*} x(X) &= -\frac{E}{\rho_{0} g} \left( \ln \lambda - \frac{(\ln \lambda)^{2}}{2} \right) + \frac{E}{\rho_{0} g} \left( \ln \lambda_{0} - \frac{(\ln \lambda_{0})^{2}}{2} \right) \\ &= -\frac{E}{\rho_{0} g} \left( \ln \frac{\lambda(X)}{\lambda(0)} - \frac12 (\ln \lambda(X) \lambda(0)) \ln \frac{\lambda(X)}{\lambda(0)} \right) \\ &= -\frac{E}{\rho_{0} g} \left(\ln \frac{\lambda(X)}{\lambda(0)} \right) \left( 1 - \frac12 \ln (\lambda(X) \lambda(0)) \right) \end{align*} $$

主ストレッチの解

まず、重力加速度 $g$ の符号により場合分けする。 $k={\rho_{0} |g| (L-X) }/ E$ とすると、主ストレッチに関する非線形方程式 (1) は以下のように表される。 $$ \begin{align} \frac{\ln \lambda}{\lambda} = \begin{cases} k & (g \ge 0) \ -k & (g < 0) \end{cases} \end{align} $$ はじめに $g>0$ のとき、左辺 $\ln\lambda / \lambda$ は $\lambda = e$ のとき最大値 $1/e$ を取るので、 $$ \begin{align} k = \frac{\rho_0 |g|}{E} (L - X) \le \frac{\rho_0 |g| L}{E} \le \frac{1}{e} \end{align} $$ を満足するパラメータを入力する必要がある。主ストレッチに関する非線形方程式 (2) を変形することで、以下の式が得られる。 $$ \begin{align*} &\lambda = e^{\lambda k} \\ &\rightarrow -k \lambda e^{-k \lambda} = -k \\ &\rightarrow W(z) e^{W(z)} = z \end{align*} $$ ここで、$z = -k$、$W(z) = z \lambda$ とおくと、上式は Lambert の W 関数 の定義式となる。 いま $z = -k \le 0$ であるため、実数解をもつ $z$ の範囲は $-1/e \le z < 0$ となるが、これは式 (3) に整合している。 対応する Lambert の W 関数の値は主枝 $W_{0}(z)$ と分枝 $W_{-1}(z)$ の2通りが存在するが、今回は主枝 $W_{0}(z)$ を対象とする。

一方で、 $g<0$ のときは任意の入力パラメータに対して解を持ち、上と同様の式変形により次式を得る。 $$ \begin{align*} &\lambda = e^{-\lambda k} \\ &\rightarrow k \lambda e^{k \lambda} = k \\ &\rightarrow W(z) e^{W(z)} = z \end{align*} $$

以上を整理すると、主ストレッチは次式で与えられる。 $$ \lambda = \frac{W(z)}{z}, \quad z = \begin{cases} -k &= -\frac{\rho_{0} |g|}{E} (L - X) & (g \ge 0) \\ k &= \frac{\rho_{0} |g|}{E} (L - X) & (g < 0) \end{cases} $$

具体例

具体的な数値解析例として、以下の入力データを用いる。

このとき、主ストレッチに関する非線形方程式 $\frac{\ln \lambda}{\lambda} = \frac{\rho_{0} g}{E} (L - X)$ をもとに、 $X = 0$ における主ストレッチの値を数値的に求めると次の値が得られる。 $$\lambda(X=0) \approx 0.567$$ この値を代入して、現在位置 $x(X)$ と応力分布 $\sigma(X)$ をプロットしたものを以下に示す。 解析モデルとして、長さ $L$ の1次元棒を一辺 0.02 m の立方体でメッシュ分割し、長手方向以外は1要素でメッシュ分割した有限要素法（FEM）により実施した。 FEMQ1は線形基底関数を、FEMQ2は2次B-spline基底関数を使用した。
線形基底関数のFEMでは、基底関数勾配が要素間で不連続になるため、応力値は要素内一定となり、要素内平均値が理論解と一致するような分布となっている。一方、2次B-spline基底関数のFEMでは、基底関数勾配が要素間で連続になるため、応力分布が滑らかに表現されていることがわかる。

理論解の計算方法

理論解は以下のpythonコードにより計算した。

import numpy as np
import mpmath
from scipy.special import lambertw
import matplotlib.pyplot as plt
import pandas as pd

def solve_for_lam(X_array, coeff, L):
    """
    Solves the equation (ln(lam))/lam - coeff * (L - X) = 0 for lam for each value in X_array.
    """
    z = np.real(-coeff * (L - X_array))
    W = np.real(lambertw(z, k=0, tol=1e-8))

    # Create an array of ones as a base
    lam_solutions = np.ones_like(z, dtype=float)

    # Perform division only where C is not zero
    np.divide(W, z, out=lam_solutions, where=(z != 0.0))

    return lam_solutions

def solve_for_x(lam, coeff):
    lam0 = lam[0]
    print("lambda at X=0:", lam0)

    x_sol = -1.0/coeff * np.log(lam / lam0) * (1.0 - 0.5*np.log(lam*lam0))

    return np.array(x_sol)

# Main execution
if __name__ == "__main__":
    # Input data (Poisson's ratio = 0.0)
    E   = 20.0e0   # Young's modulus
    rho = 1.0   # Initial density
    g   = -1.0  # Gravitational acceleration
    L   = 1.0   # Initial length of the bar

    coeff = rho * g / E

    # Consider only the main branch for Lambert W function
    if -np.exp(-1.0) > -coeff*L:
        print("No solution exists for the given parameters.")
        exit()

    # Define array X
    Xinit = np.linspace(0.0, L, 100)
    
    # Calculate lam for each X
    lam_array = solve_for_lam(Xinit, coeff, L)

    # Stress in x direction
    stress = E * np.log(lam_array) / lam_array

    # Calculate current x(X)
    x_sol = solve_for_x(lam_array, coeff)
    displ = x_sol - Xinit

    # Plot lam(X) distribution
    plt.figure(figsize=(10, 6))
    plt.plot(Xinit, lam_array, marker='.', linestyle='-')
    plt.xlabel(r"$X$")
    plt.ylabel(r"$\lambda$")
    plt.title(r"$\lambda$ s.t. $(\ln\lambda)/\lambda + (\rho*g/E)*(L-X) = 0$")
    plt.grid(True)
    plt.savefig("figs/X_vs_lambda.jpg")  

    # Plot displacement(X) distribution
    plt.figure(figsize=(10, 6))
    plt.plot(Xinit, displ, marker='.', linestyle='-')
    plt.xlabel(r"$X$")
    plt.ylabel(r"$u(X)$")
    plt.title(r"Displacement $u(X)$")
    plt.grid(True)
    plt.savefig("figs/X_vs_displacement.jpg")  

    # Plot stress distribution
    plt.figure(figsize=(10, 6))
    plt.plot(Xinit, stress, marker='.', linestyle='-')
    plt.xlabel(r"$X$")
    plt.ylabel(r"$\sigma(X)$")
    plt.title(r"Stress $\sigma(X)$")
    plt.grid(True)
    plt.savefig("figs/X_vs_stress.jpg")  

    # Plot stress-strain curve
    plt.figure(figsize=(10, 6))
    plt.plot(np.log(lam_array)/lam_array, stress, marker='.', linestyle='-')
    plt.xlabel(r"$\ln\lambda/\lambda$")
    plt.ylabel(r"$\sigma(X)$")
    plt.title(r"Stress-strain")
    plt.grid(True)
    plt.savefig("figs/ss_curve.jpg")  

    plt.show()

    df = pd.DataFrame({
        'X': Xinit,
        'displacement': displ,
        'lambda': lam_array,
        'stress': stress
    })  
    df.to_csv("results_1D_bar_finite_deformation.csv", index=False)