뭐 이렇게 장간조립교같이 생긴 트러스 다리의 3D 유한요소해석을 Ansys로 진행했습니다. (스크린샷이 따로없고 그냥 컴퓨터 화면 촬영한거라 화질구지네요;; 암튼 오늘 메인은 이게 아니니 대충 넘어갑시다.)
이제 3D 해석결과랑 2D 해석결과를 비교해봅시다. 2D유한요소해석 과정에선 구조물의 강성행렬, 변위행렬을 구하는데요,
보통은 Matlab을 쓰지만 저는 디게 특이하게 Python을 써보겠습니다. (이건 제가 Python을 연습할겸 해보는거고, 보통의 상황에선 Matlab을 쓰시는걸 추천드립니다. 행렬의 가독성이 훨씬 좋습니다.)
일단 강성행렬을 생성하기전 기본적인 2D해석 이론에 대해 간단하게 읊고 넘어갑시다.
이론을 쓱- 보면 코드로 구현하는건 어렵지않습니다....
오늘 메인은 강성행렬 프로그래밍 구현이니 유도과정은 생략하고, 결과와 구현에 중요한 개념을 위주로 다뤄봅시다
위는 트러스 요소의 강성행렬입니다. k는 요소의 등가강성으로, 요소의 단면적과 탄성계수의 곱을 요소의 길이로 나눈값과 같습니다.
유한요소 2D해석에서 중요한 개념은, 요소의 강성행렬을 모두 더해 전체 강성행렬을 조합하는데 있습니다.
요소마다 국부좌표계가 다르기때문에 조합할때 국부좌표계를 전체좌표계로 변환해주는 변환행렬이 필요합니다.
프레임 요소 해석에서는 변환행렬을 따로 곱해주는 과정이 수반되지만 트러스 요소는 강성행렬 자체에 각도정보가 포함 되어있기 때문에 따로 변환행렬을 곱해주는 과정없이 요소별로 각도(theta)값만 계산해주며 다 더해주면됩니다.
이때, 각 요소강성행렬의 전체강성행렬에서의 위치가 중요한데요, 이 위치는 요소가 걸쳐있는 절점번호에 따라 결정할수있습니다..
급하게 그리느라 그림이 약간 거렁뱅이처럼 그려졌는데 의미만 파악하셨다면 ok입니다.
위 그림에서 3,5번 절점에 걸린 트러스 요소를 생각해봅시다. 행렬 위와 우측에 표시된 변위 u는 해당 강성행렬 요소가 기여하는(?) 변위값입니다. 3,5번 절점에 걸려있다면 강성행렬이 기여하는 변위는 위 그림과 같이 되겠지요. 중요한건 이 요소 강성행렬을 전체 강성행렬에 합칠때 이 절점번호 그대로 들어간다는겁니다. 실제 전체 강성행렬은, 절점번호 3,5가 바로 붙어있지 않죠. 그렇다면 4x4강성행렬을 2x2강성행렬 4개로 쪼개어 맞는 위치에 그대로 넣어주시면 되는겁니다.
오른쪽 전체 강성행렬 그림에서 빨간색으로 색칠된 영역에 나눠져 들어간다는겁니다. 결국 각 요소를 전체 강성행렬에서 정확히 맞는 절점의 변위에 매칭시켜주는게 중요합니다. (이 과정이 코딩에서 핵심적인 부분입니다.)
이제 코딩을 짜기전 논리구조를 생각해봅시다.
위의 설계한 장간조립교의 경우 물성치와 단면적이 모두 같습니다. 게다가 길이와 각도가 같은 트러스들이 반복되는 구조이죠. 이렇다면 이 트러스의 타입을 나눠준다음 각 타입별 요소 강성행렬을 구해줍니다.
이제 타입별로 트러스가 걸치는 절점번호를 부여해주어야 합니다.
생성된 4x4크기의 요소 강성행렬은 모두 2x2강성행렬 4개로 쪼개어준다음 각 쪼개준 행렬을 '차례로' 지정된 절점번호에 해당하는 전체행렬에서의 행/열에 추가시킵니다. (예를들어 절점번호가 3인 경우에는 전체 행렬의 (4,4)-(5,5)위치에 2x2 행렬을 넣어줍니다.)
모든 요소에 대해 이 과정을 반복합니다. 그럼 끝~~
트러스 구조물의 2D모습입니다. 위에서 말했듯, 같은 요소 강성행렬을 갖는 트러스 타입이 반복되는 구조입니다.
타입을 몇개 설정해준후 타입의 절점에 따라 전체 강성행렬에 조합시키는 프로그램을 짜봅시다.
import numpy as np
import math
#총 5가지 타입의 프레임 존재
stiffness_matrix_real = np.zeros((54,54))
#type1-각도 180도 요소
A_1=100 #요소 단면적
E_1=200000 #탄성계수
L_1=2500 #요소 길이
k1 = A_1*E_1/L_1 #type1의 등가강성
deg1 = 180
angle1 =math.radians(deg1)
c1 = math.cos(angle1)
s1 = math.sin(angle1)
t1_1 = k1*pow(c1,2)
t1_2 = k1*s1*c1
t1_3 = k1*pow(s1,2)
matrix_type1 = np.array([[t1_1,t1_2,-t1_1,-t1_2],
[t1_2,t1_3,-t1_2,-t1_3],
[-t1_1,-t1_2,t1_1,t1_2],
[-t1_2,-t1_1,t1_2,t1_1]])
#type2-각도 90도 요소
A_2=100 #요소 단면적
E_2=200000 #탄성계수
L_2=5000 #요소 길이
k2 = A_2*E_2/L_2 #type1의 등가강성
deg2 = 90
angle2=math.radians(deg2)
c2 = math.cos(angle2)
s2 = math.sin(angle2)
t2_1 = k2*pow(c2,2)
t2_2 = k2*s2*c2
t2_3 = k2*pow(s2,2)
matrix_type2 = np.array([[t2_1,t2_2,-t2_1,-t2_2],
[t2_2,t2_3,-t2_2,-t2_3],
[-t2_1,-t2_2,t2_1,t2_2],
[-t2_2,-t2_1,t2_2,t2_1]])
#type3-각도 63도 요소
A_3=100 #요소 단면적
E_3=200000 #탄성계수
L_3=math.sqrt(2500*2500+5000*5000) #요소 길이
k3 = A_3*E_3/L_3 #type1의 등가강성
deg3 = 63.435
angle3=math.radians(deg3)
c3 = math.cos(angle3)
s3 = math.sin(angle3)
t3_1 = k3*pow(c3,2)
t3_2 = k3*s3*c3
t3_3 = k3*pow(s3,2)
matrix_type3 = np.array([[t3_1,t3_2,-t3_1,-t3_2],
[t3_2,t3_3,-t3_2,-t3_3],
[-t3_1,-t3_2,t3_1,t3_2],
[-t3_2,-t3_1,t3_2,t3_1]])
#type4-각도 116도 요소
A_4=100 #요소 단면적
E_4=200000 #탄성계수
L_4=math.sqrt(2500*2500+5000*5000) #요소 길이
k4 = A_4*E_4/L_4 #type1의 등가강성
deg4 = 116.565
angle4=math.radians(deg4)
c4 = math.cos(angle4)
s4 = math.sin(angle4)
t4_1 = k4*pow(c4,2)
t4_2 = k4*s4*c4
t4_3 = k4*pow(s4,2)
matrix_type4 = np.array([[t4_1,t4_2,-t4_1,-t4_2],
[t4_2,t4_3,-t4_2,-t4_3],
[-t4_1,-t4_2,t4_1,t4_2],
[-t4_2,-t4_1,t4_2,t4_1]])
sub180_1 = matrix_type1[0:2,0:2]
sub180_2 = matrix_type1[0:2,2:4]
sub180_3 = matrix_type1[2:4,0:2]
sub180_4 = matrix_type1[2:4,2:4]
sub90_1 = matrix_type2[0:2,0:2]
sub90_2 = matrix_type2[0:2,2:4]
sub90_3 = matrix_type2[2:4,0:2]
sub90_4 = matrix_type2[2:4,2:4]
sub63_1 = matrix_type3[0:2,0:2]
sub63_2 = matrix_type3[0:2,2:4]
sub63_3 = matrix_type3[2:4,0:2]
sub63_4 = matrix_type3[2:4,2:4]
sub116_1 = matrix_type4[0:2,0:2]
sub116_2 = matrix_type4[0:2,2:4]
sub116_3 = matrix_type4[2:4,0:2]
sub116_4 = matrix_type4[2:4,2:4]
num = int(input("전체 요소의 갯수를 입력하시오"))
for i in range(num):
ty_select = int(input("요소의 각도 타입을 입력하시오(1-180deg, 2-90deg, 3-63deg, 4-116deg"))
a, b = map(int, input("해당 요소의 절점번호를 작은 순서대로 입력하시오(띄어쓰기로 구분): ").split())
if ty_select == 1:
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(a-1):2*(a-1)+2] += sub180_1
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(b-1):2*(b-1)+2] += sub180_2
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(a-1):2*(a-1)+2] += sub180_3
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(b-1):2*(b-1)+2] += sub180_4
elif ty_select == 2:
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(a-1):2*(a-1)+2] += sub90_1
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(b-1):2*(b-1)+2] += sub90_2
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(a-1):2*(a-1)+2] += sub90_3
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(b-1):2*(b-1)+2] += sub90_4
elif ty_select == 3:
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(a-1):2*(a-1)+2] += sub63_1
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(b-1):2*(b-1)+2] += sub63_2
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(a-1):2*(a-1)+2] += sub63_3
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(b-1):2*(b-1)+2] += sub63_4
elif ty_select == 4:
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(a-1):2*(a-1)+2] += sub116_1
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(b-1):2*(b-1)+2] += sub116_2
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(a-1):2*(a-1)+2] += sub116_3
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(b-1):2*(b-1)+2] += sub116_4
else:
None
np.set_printoptions(threshold=np.inf) #강성행렬을 생략없이 전체 출력
print(stiffness_matrix_real)
전체 짜본 코딩이구요 처음부터 침착하게 설명해보겠읍니다. 주어진 상황에 대해서만 적용시킬수있으며, 트러스구조의 형태가 바뀐다면, 바뀐 물성치와 각도값을 수정해서 입력해주시면 됩니다. (아마 컴공 전공자분한텐 턱한대 맞을수도있는 코딩이네요..;; 기계공 전공인데 미안합니다.)
import numpy as np
import math
#총 5가지 타입의 프레임 존재
stiffness_matrix_real = np.zeros((54,54))
#type1-각도 180도 요소
A_1=100 #요소 단면적
E_1=200000 #탄성계수
L_1=2500 #요소 길이
k1 = A_1*E_1/L_1 #type1의 등가강성
deg1 = 180
angle1 =math.radians(deg1)
c1 = math.cos(angle1)
s1 = math.sin(angle1)
t1_1 = k1*pow(c1,2)
t1_2 = k1*s1*c1
t1_3 = k1*pow(s1,2)
matrix_type1 = np.array([[t1_1,t1_2,-t1_1,-t1_2],
[t1_2,t1_3,-t1_2,-t1_3],
[-t1_1,-t1_2,t1_1,t1_2],
[-t1_2,-t1_1,t1_2,t1_1]])여기서 stiffness_matrix는 전체 강성행렬입니다. ((54,54))는 54x54인 2차원 행렬을 의미합니다. 위의 구조물 절점이 27개이기 때문에 전체 구조물의 자유도는 54개입니다. stiffness_matrix를 54x54행렬로 생성해주고 0으로 채워줍시다. 차후에 계산하며 그 속을 알차게 채워줄겁니다.
밑에는 type1에 대한 요소 강성행렬을 구하는 과정입니다. 이 트러스 구조에서 요소타입은 각도와 길이에 따라 총 5가지 타입으로 분류되는데요, 그중 각도가 180도이고 길이가 2500mm인 타입에 대한 강성행렬을 생성합니다.
나머지 타입에 대해서도 동일한 과정을 거칩니다. 설명 생략하겠읍니다.
sub180_1 = matrix_type1[0:2,0:2]
sub180_2 = matrix_type1[0:2,2:4]
sub180_3 = matrix_type1[2:4,0:2]
sub180_4 = matrix_type1[2:4,2:4]
sub90_1 = matrix_type2[0:2,0:2]
sub90_2 = matrix_type2[0:2,2:4]
sub90_3 = matrix_type2[2:4,0:2]
sub90_4 = matrix_type2[2:4,2:4]
sub63_1 = matrix_type3[0:2,0:2]
sub63_2 = matrix_type3[0:2,2:4]
sub63_3 = matrix_type3[2:4,0:2]
sub63_4 = matrix_type3[2:4,2:4]
sub116_1 = matrix_type4[0:2,0:2]
sub116_2 = matrix_type4[0:2,2:4]
sub116_3 = matrix_type4[2:4,0:2]
sub116_4 = matrix_type4[2:4,2:4]이 과정에 제가 위에서 설명드린, 4x4행렬을 2x2행렬로 나누는 과정입니다. 이때 sub1~4가 전체 강성행렬에서 절점번호 작은 순서대로 들어가므로 그점만 유의하신다면 좋은 결과를 기대해볼수 있겠읍니다..
num = int(input("전체 요소의 갯수를 입력하시오"))
for i in range(num):
ty_select = int(input("요소의 각도 타입을 입력하시오(1-180deg, 2-90deg, 3-63deg, 4-116deg"))
a, b = map(int, input("해당 요소의 절점번호를 작은 순서대로 입력하시오(띄어쓰기로 구분): ").split())
if ty_select == 1:
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(a-1):2*(a-1)+2] += sub180_1
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(b-1):2*(b-1)+2] += sub180_2
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(a-1):2*(a-1)+2] += sub180_3
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(b-1):2*(b-1)+2] += sub180_4
elif ty_select == 2:
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(a-1):2*(a-1)+2] += sub90_1
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(b-1):2*(b-1)+2] += sub90_2
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(a-1):2*(a-1)+2] += sub90_3
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(b-1):2*(b-1)+2] += sub90_4
elif ty_select == 3:
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(a-1):2*(a-1)+2] += sub63_1
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(b-1):2*(b-1)+2] += sub63_2
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(a-1):2*(a-1)+2] += sub63_3
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(b-1):2*(b-1)+2] += sub63_4
elif ty_select == 4:
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(a-1):2*(a-1)+2] += sub116_1
stiffness_matrix_real[2*(a-1):2*(a-1)+2,2*(b-1):2*(b-1)+2] += sub116_2
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(a-1):2*(a-1)+2] += sub116_3
stiffness_matrix_real[2*(b-1):2*(b-1)+2,2*(b-1):2*(b-1)+2] += sub116_4
else:
None
np.set_printoptions(threshold=np.inf) #강성행렬을 생략없이 전체 출력
print(stiffness_matrix_real)마지막으로 절점의 갯수와 요소타입, 절점번호를 입력받아 전체 강성행렬을 조합해 출력하는 단계입니다.
이 코딩의 좋은점은, 절점번호를 입력값으로 받기때문에, 절점번호 순서가 달라져도 유동적인 대처가 가능하다는겁니다. (다만 입력과정이 빡쎄긴합니다. 그러나 Matlab은 절점번호 순서가 바뀌면 코딩을 다 갈아엎어야하지만 얘는 그냥 입력값만 다르게 주면 되는 이점이 있습니다.)
요소의 각도 타입을 입력받고, 요소가 걸쳐있는 절점을 띄어쓰기로 구분하며 입력받으면,
지정된 타입의 요소의 sub matrix (2x2로 나눠진 부분행렬)가 입력받는 절점번호에 해당하는 전체 강성행렬에서의 위치에 더해지는 방식입니다. 이를통해 전체 강성행렬을 구할수있습니다.
마지막에 np.set_printoptions을통해 전체 강성행렬을 '전체 다' 출력할수있습니다. 54x54행렬은 보통 파이썬에서 생략되어 출력되는데 그러면 안돼잖아!!!!!!!! 그래서 전체 다 출력해줄겁니다. (컴 사양에 따라 조금 아파할수도있습니다.아주 조금.)
출력값을 보다보면 가끔 9.79717439e-13과 같이, 0이 들어갈자리에 왠 값이 출력되는것을 볼 수 있습니다. 이 값은 파이썬에서 deg각도를 rad값으로 변환하여 cos에 대입하면 생기는 미세한 오차때문에 발생하는 값입니다.
우리의 상식으론 cos(90deg) = 0 이지만, 파이썬에서 계산해보면 deg가 rad로 변환되는 과정에서 (파이썬은 deg계산을 못합니다. rad로 변환후 계산해줘야합니다.)발생한 미세한 각도 오차가 cos값을 정확히 0이아닌 0에 아주 가까운 미세한 작은값을 출력하게 합니다. 그 값으로 인해 발생한 오차입니다.
그러므로 행렬에서 e-10이상의 터무니없이 작은값이 출력된다면 그 값은 그냥 0으로 보시면 됩니다.
근데 글을 다 쓰고나니 잘못된게 보이네요. 분명 타입이 5가지라고 했는데, 왜 위 코딩에선 4가지밖에 없지..?
알고보니 270deg 요소를 빼먹었읍니다.. 유한요소 2D해석에서 90deg와 270deg는 매우매우 다른값이므로 구분해주셔야 합니다...저의 불찰이네요...
암튼 이런식으로 Python으로도 2D 유한요소 해석을 할 수 있다!!! 이것만 알아두시면 좋을거같습니다.
위 코딩은 구조물의 상황에 맞게 변형하여 사용하시면 되겠습니다.
뻘글 봐주셔서 감사합니다
'🔧 Mechanical Engineering > dynamics🦝' 카테고리의 다른 글
| (동역학) 일과 에너지 (0) | 2024.02.28 |
|---|---|
| (동역학) 각운동량과 만유인력의 법칙 (0) | 2024.02.27 |
| (동역학) 뉴턴 제2법칙과 선형운동량 (0) | 2024.02.24 |
| (동역학) 질점의 곡선운동 (0) | 2024.02.23 |
| (동역학) 질점의 직선운동 (0) | 2024.02.22 |
