MATLAB으로 유한요소 해석하기 (보 응력 해석)(2)- 코딩하기

안녕하세요~ 매트랩으로 유한요소 해석하기 보응력 해석편 2탄!!

코딩부분을 다뤄볼텐데요.

여기서는 매트랩의 꽃인 매트릭스가 주가 될 것으로 보입니다 ㅎㅎ

 

먼저

파라미터들을 구성해볼까요??

저번 시간에 말씀드린 예제에서

E1 = 1Mpa   E2=2Mpa   E3=3Mpa

A1=100mm^2   A2=25mm^2   A3=400mm^2

L1=20mm   L2=50mm   L3=10mm

F2=20N   F3=-30

 

를 매트랩에 접목시켜 보겠습니다.

1. 기본 수치 지정

각 탄성계수에 해당하는 값들을 적어주고.. 중요한건 주석에다가 

단위들을 잘 적어주는 것이 정말 정말 중요합니다. 

다들 아시겠지만 단위 하나가 설계나 해석 영역에선 결과가 

어마 어마 어마 하게 달라지기 때문이죠^^

저는 탄성계수의 단위를 기본 단위인 Pa로 정했습니다.

단면적들도 역시나 m^2 단위를 사용하고.. 간혹가다 편의상으로 mm^2 을 사용하는 경우도 있는데 저는 코딩할때는 이게 더 편하더라구요

그다음 길이 L1,L2,L3,... 그리고 힘 F2, F3을 정해줬습니다.

여기까지는 그다지 어렵지 않습니다.

2. 구속조건 지정

그 다음 구속조건을 걸어줬습니다.

이 부분은 전시간 (1) 에 부분을 이해했다면 쉽게 이해할 수 있는 부분이죠? 벡터로 표현했습니다. 물체의 성질과 거동으로 인한 작용 식 기억나시져?

3. 노드 지정

노드를 지정해주기 위해서는 좌표계 설정이 필요한데

저는 카르테시안 좌표계(1-D)를 기준으로 할겁니다.

오른쪽을 + 왼쪽을 -라고 정하면

그래서 첫 번째 노드의 위치는 x=0

두 번째 노드의 위치는 L1 만큼 x=+L1

세 번째 노드의 위치는 L1+L2 만큼 x=L1+L2

네번째 노드의 위치는 L1+L2+L3 즉 x=L1+L2+L3 

4. 각 요소 지정 및 특성 입력

각 요소들을 담아주는 행렬을 설정하고 그에 해당되는 요소의 특성들을 담아내겠습니다. 그 특성들은 탄성계수와 단면적이 되겠네요.

첫 번째 열이 의미하는 것은 각 요소의 번호

두 번째 열은 그 요소에 연결되어 있는 노드의 작은 쪽 번호

세 번째 열은 그 요소에 연결되어 있는 노드의 큰 쪽 번호

네 번째, 다섯 번째 열은 보시는 것 처럼 탄성계수와 단면적의 넓이가 되겠네요

왜 이렇게 했냐면 다음 계산 파트로 들어갈건데

계산파트에서의 범용성을 위해서(?) 라고 할 수 있겠네요. 즉 위 1~4번 과정을 충실히 입력하기만 하면

계산 파트는 따로 수정하지 않아도 된다는 편리함이 있거든요 ^^^^

 

자!

그럼 계산파트로 넘어가겠습니다

 

1. Start der Berechnung( 계산 시작을 알림)

이 부분에서는 요소가 총 몇개가 있는지

노드는 몇개가 있는지를 파악하는 단계입니다.

그래서 size 함수를 통해 요소와 노드를 파악하죠.

dof는 자유도를 의미하는데 노드의 갯수와 동일하죠.

2. Elementsteifigkeitsmatrix( 각 요소별 특성 행렬 설립)

각 요소별 특성 행렬을 만들기 위해서

Ke라는 이름을 붙여주었고 K행렬 element라는 뜻입니다.. 줄여서 하하

저번시간에 언급했었던 K행렬e 를 그대로 적어줬습니다.

+1  -1

-1  +1  이렇게 되고 길이는 요소 행렬(e)에서 노드 번호를 찾아서

노드 행렬(n)로 가서 길이를 찾는 알고리즘 이죠

마찬가지로 vorfaktor라고 된 부분은 E*A/L 을 나타내는 부분으로 이해하시면 됩니다.

그렇게 요소별 특성 행렬은 Ke{i}=vorfaktor*Ke{i} 로 마무리가 됩니다.

3. Erweiterte Elementsteifigkeitsmatrix( 각 요소별 특성 행렬을 조합하기 위한 기초 단계)

이제 이를 조합하는 법만 남았네요.

즉.. 기본 요소 특성행렬을 이용해서 전체 행렬을 조합할 수 가 있어요!!

i=1,2,3,4가 모여서 전체 행렬 K를 이루게 되는 것이죠 후훗

즉

요렇게 조합이 되는 원리랍니다.^^ 복잡해 보였지만 또 그렇지도 않죠?

그리하여..

저렇게 조합을 할수 있게 각 요소 행렬이 지금 2x2행렬 저 윗 행렬 폼을 갖추기 위해서 4x4로 만들어야겠죠?

그래서 저 dof 라는 부분이 존재 했던 겁니다. 

Ke_erw{i}을 4x4 의 zero 행렬로 만들어 준다음

윗 그림 처럼 추가를 해주는 식으로 나가면 되겠죠

각 요소들의 특성을 담은 행렬(e)에서 보시면 두번째 열이 해당 요소(Element)에서 노드번호가 작은 노드를 뜻하고 이 부분을 j라고 할게요

마찬가지로 세번째 열이 각 요소에서 노드번호가 큰 부분을 뜻하니까 이 부분을 k라고 하죠.

그렇다면 예를 들어서

i=1이 될때 -> j=e(1,2) -> 1 이므로

Ke_erw{1} 의 (1,1) 부분에 Ke{1}의 (1,1)부분을 추가시킨다.

또한 k 는 e(1,3) -> 2이므로

Ke_erw{1} 의 (1,2) 부분에 Ke{1}의 (1,2)부분을 추가시킨다.

계속해서

Ke_erw{1} 의 (2,1) 부분에 Ke{1}의 (2,1)부분을 추가시킨다.

Ke_erw{1} 의 (2,1) 부분에 Ke{1}의 (2,1)부분을 추가시킨다.

 

그러면 최종적으로 Ke_erw{1}이 완성이 되는 거죠

그 다음은 또 Ke_erw{2}=4x4 행렬, j=2 k=3 이렇게 되어서

Ke_erw{2}는 다음과 같이

Ke_erw{3}은 다음과 같이

이렇게 됩니다.

4. Gesamtsteifgkeitsmatrix ( 전체 구조의 특성 행렬)

이를 이렇게 다 더해주었습니다. 그러면 총 전체 구조의 특성행렬이 나오겠져??ㅎㅎㅎㅎ

그럼 드디어!!1

특성 행렬은 완성이 되었습니다.

그 다음 행렬 계산량을 줄여볼까요? 경계값 조건을 사용해서 해당부분의 행렬을 지워봅시다.

저렇게 빨간 줄로 된 부분은 행렬이 삭제된 부분이고 그렇게 되면

[2X4 행렬] x [2X1행렬] = [2X1 행렬]

이렇게 줄일 수 있겠죠

코드를 한번 봅시다

보시는 것처럼 

줄여진 K행렬을 K_red 이라고 하고 여기다 K를 입력합시다.

줄여진 F행렬을 K_red  라고 하고 여기다 F를 입력합시다

그리고 u가 0인 부분의 위치를 찾아서

find(u==0)

해당되는 K의 행과 열을 지워줍시다.

K_red= (lagerung,:)=[]; 행

K_red= (:,lagerung)=[]; 열

F_red도 마찬가지로 []을 이용해서 지워줍니다.

그 다음 역함수를 이용해서 u_red을 구할 수 있구요

F도 구하고

응력 (sigma) 도 구하고

꿩도 먹고

알도 먹고

매트랩 GUI를 통한 보 응력해석을 해보았는데요

다들 이해가 가셨는지 모르겠네요

아마도 MATLAB 조금만 하실 줄 안다면!

가능 할 거라 생각합니다. ㅎㅎ

 

 

포스팅 마무리 짓겠습니다^^

여러분 안뇨오오오옹