連立一次方程式における既知量と未知量の取扱



連立一次方程式における既知量と未知量の取扱

ここで紹介しているFEM プログラムでは,剛性マトリックスの大きさを変えずに境界条件を導入し,連立方程式を解いている.この方法を用いる場合,連立一次方程式において,未知数と既知数の入れ替えが必要であるため,その考え方を以下に示す.


簡単な事例として,以下の連立一次方程式を考える.

\begin{align} & k_{11} x_1 + k_{12} x_2 + k_{13} x_3 = f_1 \\ & k_{21} x_2 + k_{22} x_2 + k_{23} x_3 = f_2 \\ & k_{31} x_3 + k_{32} x_2 + k_{33} x_3 = f_3 \\ \end{align}

ここで,未知数が $x_1, x_3, f_2$,既知数が $f_1, f_3, x_2$だとすれば,行列演算を行うためには,以下のように書き換える必要がある.

\begin{align} & k_{11} x_1 + 0 + k_{13} x_3 = f_1 - k_{12} x_2 \\ & k_{21} x_2 - f_2 + k_{23} x_3 = 0 - k_{22} x_2 \\ & k_{31} x_3 + 0 + k_{33} x_3 = f_3 - k_{32} x_2 \\ \end{align}

上記関係を,拡張して考えると,以下の通りとなる.


元の剛性方程式
\begin{equation} \begin{bmatrix} k_{11} & \ldots & k_{1i} & \ldots & k_{1j} & \ldots & k_{1n} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ k_{i1} & \ldots & k_{ii} & \ldots & k_{ij} & \ldots & k_{in} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ k_{j1} & \ldots & k_{ji} & \ldots & k_{jj} & \ldots & k_{jn} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ k_{n1} & \ldots & k_{ni} & \ldots & k_{nj} & \ldots & k_{nn} \end{bmatrix} \begin{Bmatrix} \delta_1 \\ \vdots \\ \delta_i \\ \vdots \\ \delta_j \\ \vdots \\ \delta_n \end{Bmatrix} =\begin{Bmatrix} f_1 \\ \vdots \\ f_i \\ \vdots \\ f_j \\ \vdots \\ f_n \end{Bmatrix} \end{equation}

境界条件導入後の剛性方程式

$\delta_i$ および $\delta_j$ が既知とすれば,剛性マトリックスの $k_{ii}$ および $k_{jj}$ を $1$ とし,それ以外の $i$ 列および $j$ 列の要素をゼロとする処理を行う.また剛性マトリックスにおける $i$ 列および $j$ 列の効果を右辺に移項する必要がある.

\begin{equation} \begin{bmatrix} k_{11} & \ldots & 0 & \ldots & 0 & \ldots & k_{1n} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ k_{i1} & \ldots & 1 & \ldots & 0 & \ldots & k_{in} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ k_{j1} & \ldots & 0 & \ldots & 1 & \ldots & k_{jn} \\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\ k_{n1} & \ldots & 0 & \ldots & 0 & \ldots & k_{nn} \end{bmatrix} \begin{Bmatrix} \delta_1 \\ \vdots \\ -f_i \\ \vdots \\ -f_j \\ \vdots \\ \delta_n \end{Bmatrix} =\begin{Bmatrix} f_1 \\ \vdots \\ 0 \\ \vdots \\ 0 \\ \vdots \\ f_n \end{Bmatrix} -\delta_i \begin{Bmatrix} k_{1i} \\ \vdots \\ k_{ii} \\ \vdots \\ k_{ji} \\ \vdots \\ k_{ni} \end{Bmatrix} -\delta_j \begin{Bmatrix} k_{1j} \\ \vdots \\ k_{ij} \\ \vdots \\ k_{jj} \\ \vdots \\ k_{nj} \end{Bmatrix} \end{equation}

この操作により,剛性マトリックスは非対称となりますが,数千自由度の連立一次方程式であれば,numpy の $x = np.linalg.solve(A, b)$ で大きなストレスなく解くことができる.


実際のプログラムでの処理

実際のプログラム(構造解析)での処理を以下に示す.連立方程式を解いた後,既知変位に相当する自由度の項に,既知変位を入れ直すことを忘れないように.


# treatment of boundary conditions
for i in range(0,npoin):
    for j in range(0,nfree):
        if mpfix[j,i]==1:
            iz=i*nfree+j
            fp[iz]=0.0
            for k in range(0,n):
                fp[k]=fp[k]-rdis[j,i]*gk[k,iz]
                gk[k,iz]=0.0
            gk[iz,iz]=1.0

# solution of simultaneous linear equations
disg = np.linalg.solve(gk, fp)

# recovery of restricted displacements
for i in range(0,npoin):
    for j in range(0,nfree):
        if mpfix[j,i]==1:
            iz=i*nfree+j
            disg[iz]=rdis[j,i]

npoin:総節点数
nfree:1節点自由度
mpfix:各節点毎の拘束条件を格納した配列
fp :節点外力ベクトル
rdis :変位拘束する節点の変位量
gk :全体剛性マトリックス
disg :全体剛性方程式の解(変位)



inserted by FC2 system