LoginSignup
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@GANTZ(GANTZ G)

有限要素法：トラス構造の計算プログラム（Python)

Posted at

目的

トラス構造に対して応力・変形を計算します。

バネ

1.png

力と変位には以下の関係があります。

\begin{bmatrix}
f_{x1}\\
f_{y1}\\
f_{x2}\\
f_{y2}\\
\end{bmatrix}
=
\begin{bmatrix}
-k(\delta_{x1}-\delta_{x2})\\
0\\
k(\delta_{x1}-\delta_{x2})\\
0\\
\end{bmatrix}
=
\begin{bmatrix}
-k & 0 & k & 0 \\
0 & 0 & 0 & 0\\
k & 0 & -k & 0 \\
0 & 0 & 0 & 0\\
\end{bmatrix}

\begin{bmatrix}
\delta_{x1}\\
\delta_{y1}\\
\delta_{x2}\\
\delta_{y2}\\
\end{bmatrix}

＜補足＞
ここで、一点疑問がある$f_{y1} = f_{y2}$ではないのかということである。
しかし、$f_{y1} = 0$で十分である。なぜならばこのトラスに横方向の力がかかっていると考えることに意味がないからである。また、このトラスを通じて横方向の力が伝わることはないからである。

回転マトリクス

座標系を$\theta$回転させたときに同ベクトルを別座標系で表現する。
2.png

\vec{A}
=
\begin{bmatrix}
x_1\\
0\\
\end{bmatrix}_1
=
\begin{bmatrix}
x_1 \cos \theta\\
x_1 \sin \theta\\
\end{bmatrix}_2

といえる。座標系が違うため行列内の数値が異なるがイコールでつなげることができる。

$y$座標も入れて考えると

\begin{eqnarray}
\vec{A}
&=&
\begin{bmatrix}
x_1\\
y_1\\
\end{bmatrix}_1
=
\begin{bmatrix}
x_1 \cos \theta\\
x_1 \sin \theta\\
\end{bmatrix}_2
+
\begin{bmatrix}
-y_1 \sin \theta\\
y_1 \cos \theta\\
\end{bmatrix}_2 \\

&=&
\begin{bmatrix}
\cos \theta & -\sin \theta\\
\sin \theta & \cos \theta\\
\end{bmatrix}

\begin{bmatrix}
x_1\\
y_1\\
\end{bmatrix}_2
\end{eqnarray}

具体例:剛性マトリクス

簡易な例に対して手計算による、応力解析をしてみましょう。

3.png

部材Aについて式を立てます。

\begin{eqnarray}
\begin{bmatrix}
\cos 60 &  \sin 60 \\
-\sin 60 & \cos 60 \\
\end{bmatrix}
\begin{bmatrix}
f_{A1x}\\
f_{A1y}\\
\end{bmatrix}
&=&
\begin{bmatrix}
k&0\\
0&0\\
\end{bmatrix}

\begin{bmatrix}
\delta_1 - \delta_3\\
0\\
\end{bmatrix} \\

&=&
\begin{bmatrix}
k&0&-k&0\\
0&0&0&0\\
\end{bmatrix}

\begin{bmatrix}
\cos 60 &  \sin 60  & 0 & 0\\
-\sin 60 & \cos 60  & 0 & 0\\
0 & 0 & \cos 60 &  \sin 60 \\
0 & 0 & -\sin 60 & \cos 60 \\
\end{bmatrix}

\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{3x}\\
\delta_{3y}\\
\end{bmatrix}

\end{eqnarray}

頂点3側でも同様の式が立ちます。

\begin{eqnarray}
\begin{bmatrix}
\cos 60 &  \sin 60 \\
-\sin 60 & \cos 60 \\
\end{bmatrix}
\begin{bmatrix}
f_{A3x}\\
f_{A3y}\\
\end{bmatrix}
&=&
\begin{bmatrix}
-k&0\\
0&0\\
\end{bmatrix}

\begin{bmatrix}
\delta_1 - \delta_3\\
0\\
\end{bmatrix} \\

&=&
\begin{bmatrix}
-k&0&k&0\\
0&0&0&0\\
\end{bmatrix}

\begin{bmatrix}
\cos 60 &  \sin 60  & 0 & 0\\
-\sin 60 & \cos 60  & 0 & 0\\
0 & 0 & \cos 60 &  \sin 60 \\
0 & 0 & -\sin 60 & \cos 60 \\
\end{bmatrix}

\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{3x}\\
\delta_{3y}\\
\end{bmatrix}

\end{eqnarray}

2式をまとめると次のようになります。

\begin{bmatrix}
\cos 60 &  \sin 60  & 0 & 0\\
-\sin 60 & \cos 60  & 0 & 0\\
0 & 0 & \cos 60 &  \sin 60 \\
0 & 0 & -\sin 60 & \cos 60 \\
\end{bmatrix}

\begin{bmatrix}
f_{A1x}\\
f_{A1y}\\
f_{A3x}\\
f_{A3y}\\
\end{bmatrix}

 = 

 \begin{bmatrix}
 k&0&-k&0\\
0&0&0&0\\
-k&0&k&0\\
0&0&0&0\\
\end{bmatrix}

\begin{bmatrix}
\cos 60 &  \sin 60  & 0 & 0\\
-\sin 60 & \cos 60  & 0 & 0\\
0 & 0 & \cos 60 &  \sin 60 \\
0 & 0 & -\sin 60 & \cos 60 \\
\end{bmatrix}

\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{3x}\\
\delta_{3y}\\
\end{bmatrix}

ここで次の用に文字を置きます。

T(60) = 
\begin{bmatrix}
\cos 60 &  \sin 60  & 0 & 0\\
-\sin 60 & \cos 60  & 0 & 0\\
0 & 0 & \cos 60 &  \sin 60 \\
0 & 0 & -\sin 60 & \cos 60 \\
\end{bmatrix}

その結果、

T(60)

\begin{bmatrix}
f_{A1x}\\
f_{A1y}\\
f_{A3x}\\
f_{A3y}\\
\end{bmatrix}

 = 

 \begin{bmatrix}
 k&0&-k&0\\
0&0&0&0\\
-k&0&k&0\\
0&0&0&0\\
\end{bmatrix}

T(60)

\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{3x}\\
\delta_{3y}\\
\end{bmatrix}

もう少し整理すると


\begin{bmatrix}
f_{A1x}\\
f_{A1y}\\
f_{A3x}\\
f_{A3y}\\
\end{bmatrix}

 = 
k

T(60)^{-1}

 \begin{bmatrix}
 1&0&-1&0\\
0&0&0&0\\
-1&0&1&0\\
0&0&0&0\\
\end{bmatrix}

T(60)

\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{3x}\\
\delta_{3y}\\
\end{bmatrix}

以上から計算可能な部分を計算すると以下になる。


\begin{bmatrix}
f_{A1x}\\
f_{A1y}\\
f_{A3x}\\
f_{A3y}\\
\end{bmatrix}

 = 
k

 \begin{bmatrix}
0.25&0.43&-0.25&-0.43\\
0.43&0.75&-0.43&-0.75\\
-0.25&-0.43&0.25&0.43\\
-0.43&0.75&0.43&0.75\\
\end{bmatrix}

\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{3x}\\
\delta_{3y}\\
\end{bmatrix}

部材B,部材Cに関しても同様に計算でき次のようになります。


\begin{bmatrix}
f_{C2x}\\
f_{C2y}\\
f_{C3x}\\
f_{C3y}\\
\end{bmatrix}

 = 
k

 \begin{bmatrix}
1&0&-1&0\\
0&0&0&0\\
-1&0&1&0\\
0&0&0&0\\
\end{bmatrix}

\begin{bmatrix}
\delta_{2x}\\
\delta_{2y}\\
\delta_{3x}\\
\delta_{3y}\\
\end{bmatrix}


\begin{bmatrix}
f_{B1x}\\
f_{B1y}\\
f_{B2x}\\
f_{B2y}\\
\end{bmatrix}

 = 
k

 \begin{bmatrix}
0.25& -0.43&-0.25&0.43\\
-0.43&0.75&-0.43&-0.75\\
-0.25&0.43&0.25&-0.43\\
0.43&-0.75&-0.43&0.75\\
\end{bmatrix}

\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{2x}\\
\delta_{2y}\\
\end{bmatrix}

以上3つをまとめる。力は合成することになるため次のようになります


\begin{bmatrix}
f_{1x}\\
f_{1y}\\
f_{2x}\\
f_{2y}\\
f_{3x}\\
f_{3y}\\
\end{bmatrix}

 = 

 \begin{bmatrix}
f_{A1x}\\
f_{A1y}\\
0\\
0\\
f_{A3x}\\
f_{A3y}\\
\end{bmatrix}
+
\begin{bmatrix}
0\\
0\\
f_{C2x}\\
f_{C2y}\\
f_{C3x}\\
f_{C3y}\\
\end{bmatrix}
+

\begin{bmatrix}
f_{B1x}\\
f_{B1y}\\
f_{B2x}\\
f_{B2y}\\
0\\
0\\
\end{bmatrix}

先程の行列を入れてまとめると

\begin{bmatrix}
f_{1x}\\
f_{1y}\\
f_{2x}\\
f_{2y}\\
f_{3x}\\
f_{3y}\\
\end{bmatrix}

=
k

 \begin{bmatrix}
0.5&0&-0.25&0.433&-0.25&-0.433\\
0&1.5&0.433&-0.75&-0.433&-0.75\\
-0.25&0.433&1.25&-0.433&-1&0\\
0.433&-0.75&-0.433&0.75&0&0\\
-0.25&-0.433&-1&0&1.25&0.433\\
-0.433&-0.75&0&0&0.433&0.75\\
\end{bmatrix}

\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{2x}\\
\delta_{2y}\\
\delta_{3x}\\
\delta_{3y}\\
\end{bmatrix}

となる。

具体例：境界条件と解

3.png

\begin{eqnarray}
\delta_{3x} &=& 0\\
\delta_{3y} &=& 0\\
\delta_{2y} &=& 0\\

f_{1x} &=& 1\\
f_{1y} &=& 0\\
f_{2x} &=& 0\\
\end{eqnarray}

これを代入して変形する。

\begin{bmatrix}
1\\
0\\
0\\
f_{2y}\\
f_{3x}\\
f_{3y}\\
\end{bmatrix}

=
k
 \begin{bmatrix}
0.5&0&-0.25&0.433&-0.25&-0.433\\
0&1.5&0.433&-0.75&-0.433&-0.75\\
-0.25&0.433&1.25&-0.433&-1&0\\
0.433&-0.75&-0.433&0.75&0&0\\
-0.25&-0.433&-1&0&1.25&0.433\\
-0.433&-0.75&0&0&0.433&0.75\\
\end{bmatrix}

\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{2x}\\
0\\
0\\
0\\
\end{bmatrix}

よって、

\begin{eqnarray}
\begin{bmatrix}
\delta_{1x}\\
\delta_{1y}\\
\delta_{2x}\\
\end{bmatrix}
&=&
(k
 \begin{bmatrix}
0.5&0&-0.25\\
0&1.5&0.433\\
-0.25&0.433&1.25\\

\end{bmatrix}
)^{-1}
\begin{bmatrix}
1\\
0\\
0\\
\end{bmatrix}\\
&=&
\frac{1}{k}
\begin{bmatrix}
2.25\\
-0.144\\
0.5\\
\end{bmatrix}\\
\end{eqnarray}

変位が決定する。

\begin{bmatrix}
f_{2y}\\
f_{3x}\\
f_{3y}\\
\end{bmatrix}

=
k
 \begin{bmatrix}
0.433&-0.75&-0.433\\
-0.25&-0.433&-1\\
-0.433&-0.75&0\\
\end{bmatrix}

\frac{1}{k}
\begin{bmatrix}
2.25\\
-0.144\\
0.5\\
\end{bmatrix}\\
=
\begin{bmatrix}
0.866\\
-1.0\\
-0.866\\
\end{bmatrix}\\

外力も決定する。

プログラムで一般化

具体例でわかった通り、次の手順で計算することになる

  1. 部材ごとの剛性方程式を作成
  2. [(頂点数)*(頂点数)]行列に拡張
  3. 部材ごとの結果を合成する
  4. 外力・境界条件を代入
  5. 不要行・列を取り除く
  6. 外力から逆行列で変位が計算できる
  7. 計算される変位から未確定の外力が計算できる

これを実際にプログラム化すると以下になる

import numpy as np
import copy

class Element:
    def __init__(self,angle,k, point1, point2, n):
        # point2を原点としてangleを計算する。
        self.angle = angle
        self.k = k
        self.n = n
        self.p1 = point1 - 1
        self.p2 = point2 - 1
        
    
    def culc_k(self):
        """
        剛性方程式の作成
        [f_1x f_p1y f_2x f_2y]
        =
        KA
        [d_1x d_1y d_2x d_2y]
        """        
        theta = self.angle * np.pi/180
        
        rotation_matrix = np.array([
            [np.cos(theta),  np.sin(theta)],
            [-np.sin(theta), np.cos(theta)]
        ])
        
        zero_matrix = np.zeros((2,2))
        
        T = np.block([
            [rotation_matrix, zero_matrix],
            [zero_matrix, rotation_matrix]
            ])
        K = np.array([
           [1,0,-1,0],
           [0,0,0,0],
           [-1,0,1,0],
           [0,0,0,0]
        ])
        
        KA = self.k * np.linalg.inv(T)@ K @ T
        
        return KA
        
    def culc_K(self):
        #　pointに従って行列を拡張する。
        KA = self.culc_k();
        
        #下と右側に拡張
        K = np.vstack((KA,np.zeros((2*(self.n-2),4))))
        K = np.hstack((K,np.zeros((2*self.n,2*(self.n-2)))))
        
        if self.p1 == 1:
            if self.p2 == 0:
                K[:,[0,2*self.p1]] = K[:,[2*self.p1,0]]
                K[:,[1,2*self.p1+1]] = K[:,[2*self.p1+1,1]]
                      
                K[[0,2*self.p1],:] = K[[2*self.p1,0],:]
                K[[1,2*self.p1+1],:] = K[[2*self.p1+1,1],:]
            else:
                K[:,[2,2*self.p2]] = K[:,[2*self.p2,2]]
                K[:,[3,2*self.p2+1]] = K[:,[2*self.p2+1,3]]
                K[:,[0,2*self.p1]] = K[:,[2*self.p1,0]]
                K[:,[1,2*self.p1+1]] = K[:,[2*self.p1+1,1]]
                
                K[[2,2*self.p2],:] = K[[2*self.p2,2],:]
                K[[3,2*self.p2+1],:] = K[[2*self.p2+1,3],:]            
                K[[0,2*self.p1],:] = K[[2*self.p1,0],:]
                K[[1,2*self.p1+1],:] = K[[2*self.p1+1,1],:]
            
        else:
            K[:,[0,2*self.p1]] = K[:,[2*self.p1,0]]
            K[:,[1,2*self.p1+1]] = K[:,[2*self.p1+1,1]]
            K[:,[2,2*self.p2]] = K[:,[2*self.p2,2]]
            K[:,[3,2*self.p2+1]] = K[:,[2*self.p2+1,3]]
            
            K[[0,2*self.p1],:] = K[[2*self.p1,0],:]
            K[[1,2*self.p1+1],:] = K[[2*self.p1+1,1],:]
            K[[2,2*self.p2],:] = K[[2*self.p2,2],:]
            K[[3,2*self.p2+1],:] = K[[2*self.p2+1,3],:]
        
        #print(K)
        
        
        return K
        

#----------------設定-------------------

n = 3

element_A = Element(angle = 60,k = 1,point1 = 1,point2 = 3,n = n)

element_B = Element(angle = 120,k = 1,point1 = 1,point2 = 2,n = n)

element_C = Element(angle = 0,k = 1,point1 = 2, point2 = 3,n = n)

elements = [element_A,element_B,element_C]

#境界条件
fix_point = 3         #完全固定point
fix_point2 = (2,"y")  #x,y方向どちらか拘束

#外力
force = (1,(1,0))     #point1に（1，0)N加える


#-----------------計算--------------------

#剛性行列計算
j = 0
for i in elements:
    if j == 0:
        K = i.culc_K()
        j=1
    else:
        K = K + i.culc_K()
print("K mattrix")
print(K)


#力ベクトル代入
f= np.zeros((n*2))
f[(force[0]-1)*2 ]   = force[1][0]
f[(force[0]-1)*2 +1] = force[1][1]

print("force")
print(f)

#拘束条件に基づいて行列を削除
del_list = []
del_list.append((fix_point-1)*2)
del_list.append((fix_point-1)*2+1)

if fix_point2[1] == "x":
    del_list.append((fix_point2[0]-1)*2)
if fix_point2[1] == "y":
    del_list.append((fix_point2[0]-1)*2+1)

#大きい順にソート
del_list.sort(reverse=True)

#行列から不要な列・行を削除
#K_mof = copy.deepcopy(K)
K_mof = np.delete(K,del_list,axis = 1)
K_mof = np.delete(K_mof,del_list,axis = 0)

print("K modified")
print(K_mof)

#力不要行列削除
f_mof = np.delete(f,del_list,axis=0)

defomation = np.linalg.inv(K_mof) @ f_mof

#境界条件のために省いた要素を戻す
del_list.sort()
for i in del_list:
    defomation = np.insert(defomation,i,0,axis=0)

print("defomatin")
print(defomation)

#外力を計算する
f_cul = K@defomation

print("force")
print(f_cul)

三角形について手計算の方法と一致することが確認できた。

1
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?