Uniwersytet Zielonogórski Rok akademicki

Wydział Inżynierii Lądowej i Środowiska 2009/2010

Instytut Budownictwa

Zakład Mechaniki Budowli

METODY OBLICZENIOWE

Ćw. nr 2. Metoda elementów skończonych dla zadania jednowymiarowego.

Szymon Gucwa

Grupa 22B

1 . Sformułowanie teoretyczne problemu.

Metoda elementów skończonych (MES) jest obecnie powszechnie stosowanym narzędziem obliczeń inżynierskich. Rozwiązanie problemu za pomocą metody elementów skończonych można podzielić na następujące etapy:

1. Dyskretyzacja – konstrukcja zostaje myślowo podzielona na pewną liczbę elementów skończonych

2. Zakłada się, że te połączone są ze sobą w skończonej liczbie punktów znajdujących się na obwodach. Najczęściej są to punkty narożne. Noszą one nazwę węzłów. Poszukiwane wartości wielkości fizycznych stanowią podstawowy układ niewiadomych

3. Obiera się pewne funkcje jednoznacznie określające stan przemieszczeń analizowanej wielkości fizycznej wewnątrz elementów skończonych, w zależności od wartości tych wielkości fizycznych w węzłach. Funkcje te noszą nazwę funkcji węzłowych lub funkcji kształtu.

4. Przyjęte funkcje przemieszczeń definiują jednoznacznie stan odkształceń wewnętrznych elementu. Odkształcenia te łącznie z odkształceniami początkowymi i własnościami fizycznymi materiału określają stan naprężeń w całym elemencie

5. Zostaje określony układ sił skupionych w węzłach równoważący napięcia na brzegach elementu oraz wszelkie inne siły. Równoważne siły skupione są funkcjami przemieszczeń węzłowych i znanych obciążeń zewnętrznych działających na element

6. Zostaje utworzony podstawowy układ równań MES z przemieszczeniami węzłów dyskretyzowanej konstrukcji jako niewiadomymi

K • u = f

Gdzie:

K – macierz sztywności dla całej konstrukcji

u – macierz niewiadomych przemieszczeń węzłowych

f – macierz obciążeń węzłowych

1. Po rozwiązaniu układu równań obliczane są przemieszczenia, odkształcenia i naprężenia we wszystkich elementach skończonych.

2 . Przykład liczbowy bez użycia komputera (n=3) – zginanie belki (sześcienne funkcje kształtu)

Dane;

• materiał – drewno E = 12 000 MPa

• przekrój

• schemat

1. dyskretyzacja

1. przyjęcie elementu i funkcji kształtu

$$N_{1} = 2\frac{x^{3}}{l^{3}} - \ 3\frac{x^{2}}{l^{2}} + 1$$

$$N_{2} = \frac{x^{3}}{l^{2}} - \ 2\frac{x^{2}}{l} + x$$

$$N_{3} = - 2\frac{x^{3}}{l^{3}} + \ 3\frac{x^{2}}{l^{2}}$$

$$N_{4} = \frac{x^{3}}{l^{2}} - \ \frac{x^{2}}{l}$$

1. macierz sztywności pojedynczego elementu skończonego

$$K^{e} = \left| \begin{matrix} \begin{matrix} \frac{12EI}{l^{3}} & \frac{6EI}{l^{2}} \\ \frac{6EI}{l^{2}} & \frac{4EI}{l} \\ \end{matrix} & \begin{matrix} - \frac{12EI}{l^{3}} & \frac{6EI}{l^{2}} \\ - \frac{6EI}{l^{2}} & \frac{2EI}{l} \\ \end{matrix} \\ \begin{matrix} - \frac{12EI}{l^{3}} & - \frac{6EI}{l^{2}} \\ \frac{6EI}{l^{2}} & \frac{2EI}{l} \\ \end{matrix} & \begin{matrix} \frac{12EI}{l^{3}} & - \frac{6EI}{l^{2}} \\ - \frac{6EI}{l^{2}} & \frac{4EI}{l} \\ \end{matrix} \\ \end{matrix} \right|$$

1. globalna macierz sztywności dla całego elementu skończonego i układ równań MES.

$$\frac{12EI}{l_{1}^{3}}$$	$$\frac{6EI}{l_{1}^{2}}$$	$$- \frac{12EI}{l_{1}^{3}}$$	$$\frac{6EI}{l_{1}^{2}}$$	0	0	0	0	x	w1	=	0
$$\frac{6EI}{l_{1}^{2}}$$	$$\frac{4EI}{l_{1}}$$	$$- \frac{6EI}{l_{1}^{2}}$$	$$\frac{2EI}{l_{1}}$$	0	0	0	0		φ2		0
$$- \frac{12EI}{l_{1}^{3}}$$	$$- \frac{6EI}{l_{1}^{2}}$$	$$\frac{12EI}{l_{1}^{3}} + \frac{12EI}{l_{2}^{3}}$$	$$- \frac{6EI}{l_{1}^{2}} + \frac{6EI}{l_{2}^{2}}$$	$$- \frac{12EI}{l_{2}^{3}}$$	$$\frac{6EI}{l_{2}^{2}}$$	0	0		w3		60
$$\frac{6EI}{l_{1}^{2}}$$	$$\frac{2EI}{l_{1}}$$	$$- \frac{6EI}{l_{1}^{2}} + \frac{6EI}{l_{2}^{2}}$$	$$\frac{4EI}{l_{1}} + \frac{4EI}{l_{2}}$$	$$- \frac{6EI}{l_{2}^{2}}$$	$$\frac{2EI}{l_{2}}$$	0	0		φ4		0
0	0	$$- \frac{12EI}{l_{2}^{3}}$$	$$\frac{6EI}{l_{2}^{2}}$$	$$\frac{12EI}{l_{2}^{3}} + \frac{12EI}{l_{3}^{3}}$$	$$- \frac{6EI}{l_{2}^{2}} + \frac{6EI}{l_{3}^{2}}$$	$$- \frac{12EI}{l_{3}^{3}}$$	$$\frac{6EI}{l_{3}^{2}}$$		w5		80
0	0	$$- \frac{6EI}{l_{2}^{2}}$$	$$\frac{2EI}{l_{2}}$$	$$- \frac{6EI}{l_{2}^{2}} + \frac{6EI}{l_{3}^{2}}$$	$$\frac{4EI}{l_{2}} + \frac{4EI}{l_{3}}$$	$$- \frac{6EI}{l_{3}^{2}}$$	$$\frac{2EI}{l_{3}}$$		φ6		0
0	0	0	0	$$- \frac{12EI}{l_{3}^{3}}$$	$$\frac{6EI}{l_{3}^{2}}$$	$$\frac{12EI}{l_{3}^{3}}$$	$$- \frac{6EI}{l_{3}^{2}}$$		w7		0
0	0	0	0	$$- \frac{6EI}{l_{3}^{2}}$$	$$\frac{2EI}{l_{3}}$$	$$- \frac{6EI}{l_{3}^{2}}$$	$$\frac{4EI}{l_{3}}$$		φ8		0

1. uwzględnienie warunków brzegowych

K	*	u	=	f
1	0	0	0	0
0	1	0	0	0
0	0	$$\frac{12EI}{l_{1}^{3}} + \frac{12EI}{l_{2}^{3}}$$	$$- \frac{6EI}{l_{1}^{2}} + \frac{6EI}{l_{2}^{2}}$$	$$- \frac{12EI}{l_{2}^{3}}$$
0	0	$$- \frac{6EI}{l_{1}^{2}} + \frac{6EI}{l_{2}^{2}}$$	$$\frac{4EI}{l_{1}} + \frac{4EI}{l_{2}}$$	$$- \frac{6EI}{l_{2}^{2}}$$
0	0	$$- \frac{12EI}{l_{2}^{3}}$$	$$\frac{6EI}{l_{2}^{2}}$$	$$\frac{12EI}{l_{2}^{3}} + \frac{12EI}{l_{3}^{3}}$$
0	0	$$- \frac{6EI}{l_{2}^{2}}$$	$$\frac{2EI}{l_{2}}$$	$$- \frac{6EI}{l_{2}^{2}} + \frac{6EI}{l_{3}^{2}}$$
0	0	0	0	0
0	0	0	0	0

1	0	0	0	0	0	0	0	*	0	=	0
0	1	0	0	0	0	0	0		0		0
0	0	656,25	-281,25	-150	225	0	0		w3		60
0	0	-281,25	1125	-225	225	0	0		φ4		0
0	0	-150	-225	4200	1800	0	0		w5		80
0	0	225	225	1800	1800	0	0		φ6		0
0	0	0	0	0	0	1	0		0		0
0	0	0	0	0	0	0	1		0		0

1. wyznaczenie nieznanych przemieszczeń węzłów wg wzoru

u =  K⁻¹ * f

u =	0
	0
	0,182898948
	0,08340192
	0,077549154
	-0,110836763
	0
	0

1. wyznaczenie sił wewnętrznych wg wzoru e¹ =  K¹ * u¹, e² =  K² * u², e³ =  K³ * u³

Gdzie:

K¹, K³, K³ – macierz sztywności dla poszczególnego elementu skończonego

u¹, u², u³ – macierz wyznaczonych przemieszczeń węzłowych dla poszczególnego elementu skończonego

K^1	*	u^1	=	M/T
506,25	506,25	-506,25	506,25	*
506,25	675	-506,25	337,5
-506,25	-506,25	506,25	-506,25
506,25	337,5	-506,25	675

K^2	*	u^2	=	M/T
150	225	-150	225	*
225	450	-225	225
-150	-225	150	-225
225	225	-225	450

K^3	*	u^3	=	M/T
4050	2025	-4050	2025	*
2025	1350	-2025	675
-4050	-2025	4050	-2025
2025	675	-2025	1350

1. wykresy sił wewnętrznych

3 . Program komputerowy

1. opis przygotowania danych do programu

Dane do programu są wprowadzane bezpośrednio w kodzie źródłowym programu

• wektor długości elementów skończonych

l = (/2.0, 3.0, 1.0/)

• wektor obciążeń

FG = (/0.0, 0.0, 60.0, 0.0, 80.0, 0.0, 0.0, 0.0/)

• sztywność giętna dla danego przekroju drewnianego

ej = 337.5

• deklarowanie macierzy sztywności

Ke(1, :)=(/(12.0 * ej)/l(e)* * 3, (6.0 * ej)/l(e)* * 2, ( − 12.0 * ej)/l(e)* * 3, (6.0 * ej)/l(e)* * 2/)

Ke(2, :)=(/(6.0 * ej)/l(e)* * 2, (4.0 * ej)/l(e),(−6.0 * ej)/l(e)* * 2, (2.0 * ej)/l(e)/)

Ke(3, :)=(/(−12.0 * ej)/l(e)* * 3, ( − 6.0 * ej)/l(e)* * 2, (12.0 * ej)/l(e)* * 3, ( − 6.0 * ej)/l(e)* * 2/)

Ke(4, :)=(/(6.0 * ej)/l(e)* * 2, (2.0 * ej)/l(e),(−6.0 * ej)/l(e)* * 2, (4.0 * ej)/l(e)/)

• deklarowanie warunków brzegowych

kg(1, :)=0.

kg(:, 1)=0.

kg(1, 1)=1.

kg(2, :)=0.

kg(:, 2)=0.

kg(2, 2)=1.

kg(7, :)=0.

kg(:, 7)=0.

kg(7, 7)=1.

kg(8, :)=0.

kg(:, 8)=0.

kg(8, 8)=1.

fg(1)=0.

fg(2)=0.

fg(7)=0.

fg(8)=0.

• dane wynikowe zapisywane są do pliku wyniki.txt

OPEN (UNIT = 1, FILE = ′WYNIKI.TXT′,STATUS = ′replace′,ACTION = ′WRITE′)

1. wydruk treści programu

MODULE ModMES

contains

SUBROUTINE rozw(A,B)

IMPLICIT NONE

! argumenty

REAL, INTENT(IN OUT) :: A(:,:) ! macierz

REAL, INTENT(IN OUT) :: B(:) ! wek.prawych stron i rozw.

INTEGER:: i,j,n ! n -rzad macierzy A

LOGICAL:: ok

REAL :: element, Wiersz(SIZE(B,1)),pivtol

n=SIZE(B,1)

ok=SIZE(A,1)== n .AND. SIZE(A,2)== n

IF (.NOT.ok) THEN

WRITE(*,*)' ****** Niewlasciwe wymiary macierzy ****** '

STOP

END IF

pivtol=MAXVAL(A)*1e-10

DO j=1,n

DO i=1,j-1

A(i+1:n,j)=A(i+1:n,j)-A(i,j)*A(i+1:n,i)

END DO

i=MAXVAL(MAXLOC(ABS(A(j:n,j))))+j-1

IF (ABS(A(i,j)) < pivtol) THEN

ok=.FALSE.

RETURN

END IF

IF (i.NE.j) THEN

Wiersz=A(j,:)

A(j,:)=A(i,:)

A(i,:)=Wiersz

element=B(j)

B(j)=B(i)

B(i)=element

END IF

A(j+1:n,j)=A(j+1:n,j)/A(j,j)

END DO

DO i=1,n-1

B(i+1:n)=B(i+1:n)-B(i)*A(i+1:n,i)

END DO

DO j=n,1,-1

B(j)=B(j)/A(j,j)

B(1:j-1)=B(1:j-1)-B(j)*A(1:j-1,j)

END DO

RETURN

END SUBROUTINE rozw

!-----------------------------------------------------------

SUBROUTINE agre(w1,w2,Ke,Kg)

! agreg. mac. Ke do Kg blokami odpowiadajacymi wezlom w1 i w2

IMPLICIT NONE

REAL, INTENT(IN) :: Ke(:,:)

REAL, INTENT(IN OUT):: Kg(:,:)

INTEGER, INTENT(IN) :: w1,w2

INTEGER:: i,j,lsw,nwg,nkg,nwe,nke,W(2)

lsw=UBOUND(Ke,1)/2

W(1)=w1

W(2)=w2

DO i=1,2

DO j=1,2

nwg= (W(i)-1)*lsw+1

nkg= (W(j)-1)*lsw+1

nwe= (i-1)*lsw+1

nke= (j-1)*lsw+1

Kg(nwg:nwg+lsw-1,nkg:nkg+lsw-1)= Kg(nwg:nwg+lsw-1,nkg:nkg+lsw-1)+&

Ke(nwe:nwe+lsw-1,nke:nke+lsw-1)

END DO

END DO

RETURN

END SUBROUTINE agre

!---------------------------------------

END MODULE ModMES

PROGRAM belka

USE ModMes

IMPLICIT NONE

REAL::Ke(4,4),ej,kg(8,8),l(3),Fg(8),ue(4),Se(4)

INTEGER::e

KG=0.

l=(/2.0,3.0,1.0/)

FG=(/0.0,0.0,60.0,0.0,80.0,0.0,0.0,0.0/)

ej=337.5

Do e=1,3

Ke(1,:)=(/(12.0*ej)/l(e)**3,(6.0*ej)/l(e)**2,(-12.0*ej)/l(e)**3,(6.0*ej)/l(e)**2/)

Ke(2,:)=(/(6.0*ej)/l(e)**2,(4.0*ej)/l(e),(-6.0*ej)/l(e)**2,(2.0*ej)/l(e)/)

Ke(3,:)=(/(-12.0*ej)/l(e)**3,(-6.0*ej)/l(e)**2,(12.0*ej)/l(e)**3,(-6.0*ej)/l(e)**2/)

Ke(4,:)=(/(6.0*ej)/l(e)**2,(2.0*ej)/l(e),(-6.0*ej)/l(e)**2,(4.0*ej)/l(e)/)

CALL agre(e,e+1,ke,kg)

END DO

kg(1,:)=0.

kg(:,1)=0.

kg(1,1)=1.

kg(2,:)=0.

kg(:,2)=0.

kg(2,2)=1.

kg(7,:)=0.

kg(:,7)=0.

kg(7,7)=1.

kg(8,:)=0.

kg(:,8)=0.

kg(8,8)=1.

fg(1)=0.

fg(2)=0.

fg(7)=0.

fg(8)=0.

CALL rozw(kg,fg)

OPEN (UNIT=1,FILE='WYNIKI.TXT',STATUS='replace',ACTION='WRITE')

WRITE(1,*)'PRZEMIESZCZENIA I KAT OBROTU/WARTOŚCI SIL WEWNETRZNYCH'

WRITE(1,*)FG

DO E=1,3

Ke(1,:)=(/(12.0*ej)/l(e)**3,(6.0*ej)/l(e)**2,(-12.0*ej)/l(e)**3,(6.0*ej)/l(e)**2/)

Ke(2,:)=(/(6.0*ej)/l(e)**2,(4.0*ej)/l(e),(-6.0*ej)/l(e)**2,(2.0*ej)/l(e)/)

Ke(3,:)=(/(-12.0*ej)/l(e)**3,(-6.0*ej)/l(e)**2,(12.0*ej)/l(e)**3,(-6.0*ej)/l(e)**2/)

Ke(4,:)=(/(6.0*ej)/l(e)**2,(2.0*ej)/l(e),(-6.0*ej)/l(e)**2,(4.0*ej)/l(e)/)

UE=FG((E-1)*2+1:(E-1)*2+4)

SE=MATMUL(KE,UE)

WRITE(1,*)E,SE

END DO

STOP

END PROGRAM BELKA

3 . Rozwiązanie numeryczne zadania

PRZEMIESZCZENIA I KATY OBROTU

0.0000000 0.0000000 0.18289894 8.34019184E-02 7.75491446E-02 -0.11083675 0.0000000 0.0000000

WARTOŚCI SIL WEWNETRZNYCH

1 -50.370369 -64.444443 50.370369 -36.296295

2 9.6296320 36.296291 -9.6296320 -7.4074068

3 89.629616 7.4074030 -89.629616 82.222214

4 . Porównanie wyników

	Bez użycia komputera	Rozwiązanie numeryczne
Przemieszczenia
w1	0	0.0000000
φ2	0	0.0000000
w3	0,182898948	0.18289894
φ4	0,08340192	8.34019184E-02
w5	0,077549154	7.75491446E-02
φ6	-0,110836763	-0.11083675
w7	0	0.0000000
φ8	0	0.0000000
Siły wewnętrzne
e^1	T1	-50,37037037
	M1	-64,44444444
	T2	50,37037037
	M2	-36,2962963
e^2	T2	9,62962963
	M2	36,2962963
	T3	-9,62962963
	M3	-7,407407407
e^3	T3	89,62962963
	M3	7,407407407
	T4	-89,62962963
	M4	82,22222222