Go to Qiita Advent Calendar 2024 Top
LoginSignup
noname20220504
@noname20220504

Are you sure you want to delete the question?

Leaving a resolved question undeleted may help others!

三角形要素による2次元有限要素法の実装方法（特に拘束条件）について教えてください。

Q&A

Closed

有限要素法、FEM、Fortran、拘束条件

解決したいこと

三角形要素による2次元有限要素法の実装方法（特に拘束条件）について教えてください。

発生している問題・エラー

1m×1m、三角形2要素の片持ち梁の右端に1000Nの荷重を載荷したときに生じる変位を求めたいです。
しかしながら、拘束条件を反映できず、節点変位が過大になってしまいます。
全体剛性行列、外力ベクトルの実装までは出来たのですが、その先が分かりません。

displacement is
  -27.101770043373108        7.7433628030121326       -27.101771906018257       -19.358409568667412        9.4395269378821922E-008  -19.358409568667412       -2.5995358643449151E-008   7.7433632686734200


または、問題・エラーが起きている画像をここにドラッグアンドドロップ

### 該当するソースコード
```言語名
Fortran

program FEM2D
implicit none
integer :: i, j, k, l
real(8) :: x, y, E, nu, t, load, nodes(4,2), D(3,3), A(2), B1(3,6), B2(3,6), K1(6,6), K2(6,6),L1(3), L2(3)
real (8) :: K_global(8,8), F(8), displacement(8)
!integer, parameter :: n = 8
!integer,allocatable,dimension(:) :: tanni
!integer :: i, j, k
real :: temp, tanni(8,8)

! 材料物性値
E = 210.0e+9 !ヤング率(Pa)
nu = 0.20 ! ポアソン比
t = 0.010 ! 部材厚(m)
load = 1000.0 ! 載荷荷重(N)

! 節点情報
nodes(1,1) = 0.0 ! 節点番号1のx座標
nodes(1,2) = 0.0 ! 節点番号1のy座標
nodes(2,1) = 1.0 ! 節点番号2のx座標
nodes(2,2) = 0.0 ! 節点番号2のy座標
nodes(3,1) = 1.0 ! 節点番号3のx座標
nodes(3,2) = 1.0 ! 節点番号3のy座標
nodes(4,1) = 0.0 ! 節点番号4のx座標
nodes(4,2) = 1.0 ! 節点番号4のy座標

! コネクティビティ
!  要素1：節点1→2→3
!  要素2：節点1→3→4

! Dマトリクス
D(1,1) = 1.0*E/(1.0-nu**2.0)
D(1,2) = nu*E/(1.0-nu**2.0)
D(1,3) = 0.0
D(2,1) = nu*E/(1.0-nu**2.0)
D(2,2) = 1.0*E/(1.0-nu**2.0)
D(2,3) = 0.0
D(3,1) = 0.0
D(3,2) = 0.0
D(3,3) = (1.0-nu)/2.0*E/(1.0-nu**2.0)

! 各要素の面積A(i)(i=1,2)
A(1)= 1.0/2.0*((nodes(1,1)-nodes(3,1))*(nodes(2,2)-nodes(3,2))-(nodes(2,1)-nodes(3,1))*(nodes(1,2)-nodes(3,2)))
A(2)= 1.0/2.0*((nodes(1,1)-nodes(4,1))*(nodes(3,2)-nodes(4,2))-(nodes(3,1)-nodes(4,1))*(nodes(1,2)-nodes(4,2)))

! Bマトリクス
B1(1,1) = 1.0/(2.0*A(1))*(nodes(2,2)-nodes(3,2)) 
B1(1,2) = 0.0
B1(1,3) = 1.0/(2.0*A(1))*(nodes(3,2)-nodes(1,1))
B1(1,4) = 0.0
B1(1,5) = 1.0/(2.0*A(1))*(nodes(1,1)-nodes(2,2))
B1(1,6) = 0.0
B1(2,1) = 0.0
B1(2,2) = 1.0/(2.0*A(1))*(nodes(3,1)-nodes(3,2))
B1(2,3) = 0.0
B1(2,4) = 1.0/(2.0*A(1))*(nodes(1,1)-nodes(3,1))
B1(2,5) = 0.0
B1(2,6) = 1.0/(2.0*A(1))*(nodes(2,1)-nodes(1,1))
B1(3,1) = 1.0/(2.0*A(1))*(nodes(3,1)-nodes(2,1))
B1(3,2) = 1.0/(2.0*A(1))*(nodes(2,2)-nodes(3,2))
B1(3,3) = 1.0/(2.0*A(1))*(nodes(1,1)-nodes(3,1))
B1(3,4) = 1.0/(2.0*A(1))*(nodes(3,2)-nodes(1,2))
B1(3,5) = 1.0/(2.0*A(1))*(nodes(2,1)-nodes(1,1))
B1(3,6) = 1.0/(2.0*A(1))*(nodes(1,2)-nodes(2,2))

!-------------------------------------------------------
B2(1,1) = 1.0/(2.0A(2))(nodes(3,2)-nodes(4,2))
B2(1,2) = 0.0
B2(1,3) = 1.0/(2.0A(2))(nodes(4,2)-nodes(1,1))
B2(1,4) = 0.0
B2(1,5) = 1.0/(2.0A(2))(nodes(1,1)-nodes(3,2))
B2(1,6) = 0.0
B2(2,1) = 0.0
B2(2,2) = 1.0/(2.0A(2))(nodes(4,1)-nodes(4,2))
B2(2,3) = 0.0
B2(2,4) = 1.0/(2.0A(2))(nodes(1,1)-nodes(4,1))
B2(2,5) = 0.0
B2(2,6) = 1.0/(2.0A(2))(nodes(3,1)-nodes(1,1))
B2(3,1) = 1.0/(2.0A(2))(nodes(4,1)-nodes(3,1))
B2(3,2) = 1.0/(2.0A(2))(nodes(3,2)-nodes(4,2))
B2(3,3) = 1.0/(2.0A(2))(nodes(1,1)-nodes(4,1))
B2(3,4) = 1.0/(2.0A(2))(nodes(4,2)-nodes(1,2))
B2(3,5) = 1.0/(2.0A(2))(nodes(3,1)-nodes(1,1))
B2(3,6) = 1.0/(2.0A(2))(nodes(1,2)-nodes(3,2))

! 要素剛性行列
K1 = A(1)*t*MATMUL(MATMUL(transpose(B1), D), B1)
K2 = A(2)*t*MATMUL(MATMUL(transpose(B2), D), B2)

! 全体剛性行列
K_global(1,1) = K1(1,1) + K2(1,1)
K_global(1,2) = K1(1,2) + K2(1,2)
K_global(1,3) = K1(1,3) + 0.0
K_global(1,4) = K1(1,4) + 0.0
K_global(1,5) = K1(1,5) + K2(1,3)
K_global(1,6) = K1(1,6) + K2(1,4)
K_global(1,7) = 0.0 + K2(1,5)
K_global(1,8) = 0.0 + K2(1,6)

K_global(2,1) = K1(2,1) + K2(2,1)
K_global(2,2) = K1(2,2) + K2(2,2)
K_global(2,3) = K1(2,3) + 0.0
K_global(2,4) = K1(2,4) + 0.0
K_global(2,5) = K1(2,5) + K2(2,3)
K_global(2,6) = K1(2,6) + K2(2,4)
K_global(2,7) = 0.0 + K2(2,5)
K_global(2,8) = 0.0 + K2(2,6)

K_global(3,1) = K1(3,1) + 0.0
K_global(3,2) = K1(3,2) + 0.0
K_global(3,3) = K1(3,3) + 0.0
K_global(3,4) = K1(3,4) + 0.0
K_global(3,5) = K1(3,5) + 0.0
K_global(3,6) = K1(3,6) + 0.0
K_global(3,7) = 0.0 + 0.0
K_global(3,8) = 0.0 + 0.0

K_global(4,1) = K1(4,1) + 0.0
K_global(4,2) = K1(4,2) + 0.0
K_global(4,3) = K1(4,3) + 0.0
K_global(4,4) = K1(4,4) + 0.0
K_global(4,5) = K1(4,5) + 0.0
K_global(4,6) = K1(4,6) + 0.0
K_global(4,7) = 0.0 + 0.0
K_global(4,8) = 0.0 + 0.0

K_global(5,1) = K1(5,1) + K2(3,1)
K_global(5,2) = K1(5,2) + K2(3,2)
K_global(5,3) = K1(5,3) + 0.0
K_global(5,4) = K1(5,4) + 0.0
K_global(5,5) = K1(5,5) + K2(3,3)
K_global(5,6) = K1(5,6) + K2(3,4)
K_global(5,7) = 0.0 + K2(3,5)
K_global(5,8) = 0.0 + K2(3,6)

K_global(6,1) = K1(6,1) + K2(4,1)
K_global(6,2) = K1(6,2) + K2(4,2)
K_global(6,3) = K1(6,3) + 0.0
K_global(6,4) = K1(6,4) + 0.0
K_global(6,5) = K1(6,5) + K2(4,3)
K_global(6,6) = K1(6,6) + K2(4,4)
K_global(6,7) = 0.0 + K2(4,5)
K_global(6,8) = 0.0 + K2(4,6)

K_global(7,1) = K2(5,1)
K_global(7,2) = K2(5,2)
K_global(7,3) = 0.0
K_global(7,4) = 0.0
K_global(7,5) = K2(5,5)
K_global(7,6) = K2(5,4)
K_global(7,7) = K2(5,5)
K_global(7,8) = K2(5,6)

K_global(8,1) = K2(6,1)
K_global(8,2) = 0.0
K_global(8,3) = 0.0
K_global(8,4) = K2(6,1)
K_global(8,5) = K2(6,3)
K_global(8,6) = K2(6,4)
K_global(8,7) = K2(6,5)
K_global(8,8) = K2(6,6)

print *, 'K_global is'
print *, K_global
	
! 境界条件（拘束・荷重）の適用
F(1) = 0.0
F(2) = 0.0
F(3) = 0.0
F(4) = -1.0 * load
F(5) = 0.0
F(6) = 0.0
F(7) = 0.0
F(8) = 0.0

! 単位行列tanniの初期化
tanni = 0.0
do i = 1, 8
    tanni(i, i) = 1.0
end do

! ガウス・ジョルダン法で逆行列を求める
do i = 1, 8
    ! ピボット選択
    if (K_global(i, i) == 0.0) then
        print *, 'Matrix is singular, inverse cannot be computed.'
        stop
    end if
    temp = K_global(i, i)
    K_global(i, :) = K_global(i, :) / temp
    tanni(i, :) = tanni(i, :) / temp
    
    do j = 1, 8
        if (j /= i) then
            temp = K_global(j, i)
            K_global(j, :) = K_global(j, :) - temp * K_global(i, :)
            tanni(j, :) = tanni(j, :) - temp * tanni(i, :)
        end if
    end do
end do

! 変位を出力
!print *, 'The inverse matrix is:'
displacement = MATMUL(tanni,F)
print *, 'displacement is'
print *, displacement

end program FEM2D
### 自分で試したこと
ここに問題・エラーに対して試したことを記載してください。

0

No Answers yet.

Your answer might help someone💌

LoginSign Up