E.6 数値解析法

最新版を正確に読む場合には pdf ファイルをどうぞ。これは web 検索のための簡易旧版です。

複雑な系については次節の有限要素法を使えばいいが，1スパンの2点境界値問題であれば，微分方程式を直接数値積分して解を求めることができる。まず左端で与えられた境界条件に加えて未知の量も適当に仮定して右端に向かって数値積分する。次に，その右端の結果と右端で与えられた境界条件の差を解消するように，今度は右から左に数値積分して左端の未知量を更新する。これを繰り返すことによって解を得ることができる。非線形微分方程式であっても数値積分は容易にできる。この節には基礎式となる1階連立常微分方程式を列挙するに留め，数値積分の方法については文献[64]を参照して欲しい。なお，この1階常微分方程式は表-E.3にもまとめた。

長さ $\ell$ の梁に対する無次元化した連立微分方程式にするために以下の量と微分を定義する。

$\displaystyle z_1$	$\textstyle \equiv$	$\displaystyle \dfrac{\left(N\cos\vartheta+V\sin\vartheta\right) \ell^2}{EI},\quad z_2\equiv\dfrac{\left(-N\sin\vartheta+V\cos\vartheta\right) \ell^2}{EI},$
$\displaystyle z_3$	$\textstyle \equiv$	$\displaystyle \dfrac{M \ell}{EI},\quad z_4\equiv\dfrac{u}{\ell},\quad z_5\equi... ...equiv\vartheta, \qquad \dot{( )}\equiv\D*{( )}{\left(\slfrac{x}{\ell}\right)}$	(E.95)

と定義した。ここに $\alpha\subsc{t}'$ は，式(4.85a)で定義した $\alpha\subsc{t}$ を若干定義し直したパラメータである。これに対し，第1次近似理論の場合は，上式のうち

$\displaystyle \dot{z}_3$	$\textstyle =$	$\displaystyle y_2 \dfrac{1+\beta^2 y_1}% {1+\dfrac{\alpha\subsc{t}' \beta^2 y_1}{1+\beta^2 y_1}}$
$\displaystyle \dot{z}_4$	$\textstyle =$	$\displaystyle \left(1+\beta^2 y_1\right)\cos z_6 + \dfrac{\alpha\subsc{t}' \b... ... y_2 \sin z_6}% {1+\dfrac{\alpha\subsc{t}' \beta^2 y_1}{1+\beta^2 y_1}}-1$	(E.98)
$\displaystyle \dot{z}_5$	$\textstyle =$	$\displaystyle -\left(1+\beta^2 y_1\right)\sin z_6 + \dfrac{\alpha\subsc{t}' \... ...^2 y_2 \cos z_6}% {1+\dfrac{\alpha\subsc{t}' \beta^2 y_1}{1+\beta^2 y_1}}$

また伸びを無視した場合，まず第2次近似のつり合い式のうちの式(E.97c)を

$\displaystyle \dot{z}_3$	$\textstyle =$	$\displaystyle \dfrac{y_2}{1+\alpha\subsc{t}' \beta^2 y_1}$
$\displaystyle \dot{z}_4$	$\textstyle =$	$\displaystyle \left(1+\beta^2 y_1\right)\cos z_6 + \dfrac{\alpha\subsc{t}' \beta^2 y_2 \sin z_6}% {1+\alpha\subsc{t}' \beta^2 y_1}-1$	(E.100)
$\displaystyle \dot{z}_5$	$\textstyle =$	$\displaystyle -\left(1+\beta^2 y_1\right)\sin z_6 + \dfrac{\alpha\subsc{t}' \beta^2 y_2 \cos z_6}% {1+\alpha\subsc{t}' \beta^2 y_1}$

となる。数値積分は各種パッケージに含まれているものを利用すればいい。

E.6.2 一つの有限要素解析法

E.6.2.1 微小ひずみの枠組の中での原理

ここでは弾性の有限変位理論を取り扱っていることから，微小ひずみの仮定に基づき，簡便にかつ高精度で棒の有限変位解析を有限要素法を用いて解く手法の一つ[39]を説明する。変形の一般的な記述に極分解の定理 [51]がある。これは簡単に言うと，全体の変形勾配は回転成分とひずみ成分の積で表されるというものである。つまり微小なひずみのまま大きくたわむピアノ線のようなものを対象とする限りは，大きく見える変位成分から回転成分を除いてしまうと，実質的に変形に寄与する変位成分の大きさは非常に小さいことを示唆している。この定理を念頭に置き，次のように考えて剛性方程式を仮定しよう。

ある有限要素を取り出したとき，例えば左節点の断面のたわみ角分だけ回転した要素座標系で見ると，要素長を短くすればするほどその左右端の相対変位（たわみ角も含む）は小さくなり，その相対変位については微小変位理論の剛性方程式が成立すると考えるのである（図は文献[39]を参照のこと）。つまり，ある代表的な構造物の長さを

としたときに，外力ベクトルと変位ベクトルを

で近似する。なお簡単のためにこの章のこの節では，括弧無しの太字で行列を表している。ここで $\fat{T}$ は左の節点の断面の回転角 $\vartheta_1$ で定義される座標変換行列で， $\fat{D}$ は要素の剛体変位成分で

と置いた。ここに $\ell$ は有限要素の長さである。また剛性行列は，Timoshenko梁の幾何剛性を含む式(E.67)を用いて

$\displaystyle \fat{k}_1\supersc{l}$	$\textstyle \equiv$	$\displaystyle \left(\begin{array}{ccc} -\dfrac{\xi}{\beta^2} & 0 & 0 \\ 0 & -\... ...elta_0^2} & \dfrac{-\slfrac{1}{30}-\Delta_2}{\xi\Delta_0^2} \end{array}\right),$	(E.105)
$\displaystyle \fat{k}_2\supersc{l}$	$\textstyle \equiv$	$\displaystyle \left(\begin{array}{ccc} \dfrac{\xi}{\beta^2} & 0 & 0 \\ & \dfra... ...x{Symm.} & & \dfrac{\slfrac{2}{15}+\Delta_2}{\xi\Delta_0^2} \end{array}\right),$

である。 $k\subsc{t}$ はせん断変形の取り扱いを補正する係数（節-4.6.2に例がある）で，断面形状とPoisson比に依存する定数[16]である。あとは，式(E.103)を例えばNewton-Raphson法で解けばいい。そのときの接線剛性 $\fat{k}_t$ は式(E.103)を $\fat{d}_i$ で偏微分するだけで

$\displaystyle \fat{H}_i$	$\textstyle \equiv$	$\displaystyle \left\lfloor \fat{H}_{i1} \fat{H}_{i2} \fat{H}_{i3} \... ...{T}\fat{k}_i\fat{Q}\supersc{t}\right)\left(\fat{d}_2- \fat{d}_1-\fat{D}\right),$
$\displaystyle \fat{P}_i$	$\textstyle \equiv$	$\displaystyle \xi\fat{T}\fat{k}_i\supersc{nl}\fat{T}\supersc{t} \left(\fat{d}_2-\fat{d}_1-\fat{D}\right)/\beta^2,$
$\displaystyle \fat{C}$	$\textstyle \equiv$	$\displaystyle \left(\begin{array}{cccccc} -1 & 0 & -\sin\vartheta_1/\xi & 1 & 0... ...a_1 & \multicolumn{1}{r}{\sin\vartheta_1} & 0 \\ 0 & 0 & 0 \end{array}\right),$	(E.108)
$\displaystyle g$	$\textstyle \equiv$	$\displaystyle \left\{\left(u_2-u_1\right)/L-\left(\cos\vartheta_1-1\right)/\xi ... ...eta_1+ \left\{\left(w_2-w_1\right)/L+\sin\vartheta_1/\xi\right\}\cos\vartheta_1$

と置いた。例えば

ステップまでの解が求められているとすると，式(E.103)から不つり合い力を算定し，式(E.108)の接線剛性を用いて，

ステップの修正解が

のようにして求められる。右辺の接線剛性の逆行列がかかっている項が

ステップの不つり合い力である。収束は例えば

**図 E.3:** 片端ヒンジで片端固定の深いアーチ頂点の荷重変位曲線と変形の進展
$\begin{figure}\begin{center} \unitlength=.01mm \begin{picture}(7399,6000)(1500,-... ...1,Legend(Title) %,-1,Graphics End %E,0, % \end{picture}\end{center} \end{figure}$

**図 E.4:** 浅いアーチ頂点の荷重変位曲線
**図 E.5:** 基部に向かって引張られる柱
$\begin{figure}\begin{center} \unitlength=.01mm \begin{picture}(7247,6000)(1250,-... ...水平変位})/L$}} %,-1,Graphics End %E,0, % \end{picture}\end{center} \end{figure}$

E.6.2.2 数値解析例

図-E.3は文献[17]で対象とされた構造である。また図-E.4は文献[135]でも対象とされた浅いアーチで，AA'およびBB'に分岐経路がある。文献ではBB'は示されていないが，ここのアプローチでは見つけられた経路である。ここで定式化された接線剛性は非対称ではあるが，分岐前後で対角要素（固有値）に負の値が現れ，分岐点と荷重方向および変位方向のピークとでその数が変化する。図-E.5には，斜張橋の塔が面外に座屈することを想定したものであるが，基部に向かって引張られて，魚のかかった釣竿のようになる柱の荷重変位関係を示した。ただし，最初に少しだけ曲がった柱に載荷している。図-6.13や図-6.14等も，同じプログラムで数値的に解いて得たものである。

ここで用いたプログラムでは，メモリ節約等のためにスカイライン法[22]を用い，安定な解析のために弧長法[135]を併用した。このプログラムもまえがきに書いた方法で入手可能である。またこのプログラムは，節-E.3.4で誘導した4種類の座屈公式のそれぞれに対応した剛性行列を用いて，それぞれのモデルの解を求めることができる。このような対応を表-E.3に示した。

**表 E.3:** 面内梁理論の剛性行列と座屈公式の対応
伸び考慮			伸び無視
梁-柱モデル	モデルA	モデルB	梁-柱モデル	モデルA	モデルB
$u'= \left(1+\epsilon\right)\cos\vartheta+\gamma\sin\vartheta-1$ , $w'=-\left(1+\epsilon\right)\sin\vartheta+\gamma\cos\vartheta$
$\left(N\cos\vartheta+V\sin\vartheta\right)'+p=0$ , $\left(-N\sin\vartheta+V\cos\vartheta\right)'+q=0$
$M'-V\left(1+\epsilon\right)=0$	$M'-V\left(1+\epsilon\right)+N\gamma=0$			$M'-V+N\gamma=0$
$N=EA\epsilon$ , $M=EI\vartheta'$
第2次近似		第1次近似	第2次近似		第1次近似
$V=Gk\subsc{t}A\gamma$		$V=\left(Gk\subsc{t}A+\dfrac{N}{1+\epsilon}\right)\gamma$	$V=Gk\subsc{t}A\gamma$		$V=\left(Gk\subsc{t}A+N\right)\gamma$
式(E.32)	式(E.51)	式(E.48)	式(E.33)	式(E.52)	式(E.50)
			Euler式	改訂式	Engesser式
$\fat{K}\supersc{l}$	$\fat{K}\supersc{l}+\fat{K}\supersc{nl}\subsc{a}$	$\fat{K}\supersc{l}+\fat{K}\supersc{nl}\subsc{b}$	$\fat{K}\supersc{l}$	$\fat{K}\supersc{l}+\fat{K}\supersc{nl}\subsc{a}$	$\fat{K}\supersc{l}+\fat{K}\supersc{nl}\subsc{b}$
			$\epsilon=0$	$\epsilon=0$	$\epsilon=0$
つり合いは厳密に成立する			つり合いは近似的にしか成立しない

E.6.2.3 面内梁の最も一般的な剛性行列

誘導の詳細については文献[37]を参照のこと。有限変位問題に対しても，微小変位（線形化された）理論の剛性行列を用いることができることを前節で示したので，ここには線形化された理論，つまり梁-柱理論レベルの剛性行列を示しておく。 $\fat{f}\equiv\left\lfloor\fat{f}_1 \fat{f}_2 \right\rfloor\supersc{t}$ と定義して

という関係で剛性行列を定義する。なお，添え字mのAとBは，Timoshenko梁理論のせん断力とせん断変形の二つのモデルを区別したものである（表-E.3参照）。要素長を

としたときの，有限変位の枠組での平面骨組の剛性行列 は

$\displaystyle \fat{K}_1\supersc{l}$	$\textstyle \equiv$	$\displaystyle \left(\begin{array}{ccc} \dfrac{1}{\beta^2} & 0 & 0 \\ 0 & \dfra... ... & \dfrac{6}{(1+\epsilon)\Delta} & \dfrac{2-12\Psi}{\Delta} \end{array}\right),$
$\displaystyle \fat{K}_3\supersc{l}$	$\textstyle \equiv$	$\displaystyle \left(\begin{array}{ccc} \dfrac{1}{\beta^2} & 0 & 0 \\ & \dfrac{... ...n)\left(\slfrac{2}{15}+\Delta_2\super{m}\right)} {\Delta^2} \end{array}\right),$	(E.112)
$\displaystyle {}\super{m}\fat{K}_2\supersc{nl}$	$\textstyle \equiv$	$\displaystyle \left(\begin{array}{ccc} 0 & 0 & 0 \\ 0 & -\dfrac{6\Delta_1\supe... ...lon)\left(\slfrac{2}{15}+\Delta_2\super{m}\right)}{\Delta^2} \end{array}\right)$

$\displaystyle \Delta$	$\textstyle \equiv$	$\displaystyle 1+12\Psi, \quad \Psi\equiv\dfrac{\alpha\subsc{t}}{(1+\epsilon)^2}... ...}=\mbox{A}) \\ 1+10\Delta_2\super{m} & (\mbox{m}=\mbox{B}) \end{array}\right.,$	(E.113)
$\displaystyle \Delta_2\super{m}$	$\textstyle \equiv$	$\displaystyle \left\{\begin{array}{ll} 2\Psi-24\Psi^2 & (\mbox{m}=\mbox{A}) \ ... ...-720\Psi^2 & (\mbox{m}=\mbox{A}) \\ 1 & (\mbox{m}=\mbox{B}) \end{array}\right.$

で算定できる。Bernoulli-Euler梁の場合は， $\alpha\subsc{t}=0$ とすればいい。「伸びを考慮」した「モデルB」が，梁に対して最も正確な解を与えてくれる。

E.6 数値解析法

E.6.1 2点境界値問題の一解法

E.6.2 一つの有限要素解析法

E.6.2.1 微小ひずみの枠組の中での原理

E.6.2.2 数値解析例

E.6.2.3 面内梁の最も一般的な剛性行列