Zasada prac wirtualnych, metoda Ritza i metoda

elementów skończonych

Leszek Chodor

[e-mail]

1

Zasada prac wirtualnych i funkcjonał Lagran-
ge’a)

1.1

Zasada prac wirtualnych

Równanie prac wirtualnych wyprowadzimy dla bryły o dowolnym kształcie,
dowolnie obciążonej (z dowolnymi statycznymi warunkami brzegowymi) na
powierzchni A

q

oraz z dowolnymi kinematycznymi warunkami brzegowymi

na A

u

.

Załóżmy , że znamy rozwiązanie zagadnienie brzegowego teorii sprężystości
(ZBTS) dla tej bryły, to znaczy znamy aprężenia σ

ij

, odkształcenia ε

ij

oraz

pole przemieszczeń u

i

(i, j = 1, 2, 3).

Nadajmy przemieszczeniom u

i

małe przyrosty δu

i

, takie , aby przemieszcze-

nie po przyroście wynosiło:

e

u

i

= u

i

+ δu

i

i było przemieszczaniem możliwym, to znaczy spełniający kinematyczne wa-
runki brzegowe i będącym funkcją ciągłą co najmniej klasy C

3

. Spełnienie

przez

e

u

i

kinematycznych warunków brzegowych oznacza, że wariacje δu

i

znikają na powierzchniA

u

. Przyrosty δu

i

jako niewielkie zaburzenia stanu

równowagi w rachunku wariacyjnym nazywane są wariacjami, a w mechani-
ce nazywane są przemieszczeniami wirtualnymi.
Pracę sił zewnętrznych L na przemieszczeniach wirtualnyh nazywamy pracą
wirtualną δL:

δL =

Z Z

A

q

+A

u

q

υi

δu

i

dA +

Z Z Z

V

P

i

δu

i

dV

(1)

Wykażemy, że wirtualne praca sił zewnętrznych jest równa wirtualnej pracy
sił wewnętrznych.δW :

δW =

Z Z

V

σ

ij

δε

ij

dV

(2)

W tym celu wykonamy kilka przekształceń matematycznych wyrażenia (1)

1

Najpierw zauważmy, że w pierwszej całce (1) możemy pominąć całkowanie
po powierzchni A

u

, wobec δu

i

= 0 na tej powierzchni i zapisać, że całkujemy

po całej powierzchni A. Jeśli następnie gęstość sił zewnętrznych q

υi

wyra-

zimy w zależności od sił wewnętrznych na brzegu, zgodnie ze statycznymi
warunkami brzegowymi na A

q

(q

υi

= σ

ij

α

υj

), to otrzymamy

δL =

Z Z

A

σ

ij

α

υj

δu

i

dA +

Z Z Z

V

P

i

δu

i

dV

Po zastosowaniu do pierwszej całki, twierdzenia Greena-Gaussa-Ostrogradzkiego
o zamianie całki powierzchniowej na objętościową , otrzymamy:

δL =

Z Z Z

V

∂

∂x

j

(σ

ij

δu

i

)dV +

Z Z Z

V

P

i

δu

i

dV =

Z Z Z

V

h

∂σ

ij

∂x

j

+ P

i



δu

i

+ σ

ij

∂δu

i

∂x

j

i

dV

Ponieważ σ

ij

z założenia spełniają równania Naviera, więc pierwszy skład-

nik w wyrażeniu podcałkowym jest równy zero. Wykorzystując własność
sumacji, drugi składnik można zapisać w postaci:

δL =

Z Z Z

V

σ

ij

·

1

2



∂δu

i

∂x

j

+

∂δu

j

∂x

i



dV

Zauważmy, że wyrażenie w nawiasie kwadratowym na mocy równań Cau-
che’go jest wariacją odkształcenia. Mamy więc:

δL =

Z Z Z

V

σ

ij

δε

ij

dV = δW,

co kończy dowód zasady prac wirtualnych

δL = δW.

(3)

Zauważmy jeszcze, że równanie prac wirtualnych (3) jest słuszne dla bryły
wykonanej z materiału nie tylko liniowo sprężystego, ale również nieliniowego
i niesprężystego, ponieważ do jego wyprowadzenia nie używaliśmy prawa
Hooke’a.

1.2

Twierdzenia Lagrange’a

Twierdzenie Lagrange’a o minimum energii potencjalnej układu, można sfor-
mułować następująco:
Spośród wszystkich możliwych przemieszczeń, te są rzeczywiste, które reali-
zują minimum funkcjonału energii potencjalnej układu.
Twierdzenie to w języku rachunku wariacyjnego zapisuje się jako wymóg
zerowania wariacji δ funkcjonału energii potencjalnej układu Π. , czyli ja-
ko warunek konieczny istnienia ekstremum funkcjonału (analog kryterium
istnienia ekstremum funkcji).

δΠ = 0

(4)

2

W warunku stacjonarności (4) funkcjonał Π jest liczbą przyporządkowaną
trójce funkcji przemieszczeń u

i

. Przyrost tej liczby nazywany jest wariacją

funkcjonału δ i jest analogiem różniczki zmiennej w analizie funkcji (nie-
skończenie małej zmiany tej zmiennej
Twierdzenie Lagrange’a (4) uzyskaliśmy z (3) po uwzględnieniu, że w polu
potencjalnym naprężeń pochodna energii potencjalnej Φ względem odkształ-
cenia jest równa naprężeniom:

∂Φ

∂ε

ij

= σ

ij

.

(5)

Formułę (5) przyjmujemy jako definicję energii potencjalnej Φ. W celu zapi-
sania twierdzenia Lagrange’a w postaci (4) wprowadziliśmy oznaczenie (6),
które traktujemy jako definicję funkcjonału Lagrange’a:

Π =

Z Z Z

V

Φ(u

i

)dV −

Z Z

A

q

q

υi

u

i

dA −

Z Z Z

V

−P

i

u

i

dV.

(6)

Z definicji (5) uzyskamy wyrażenie na energię potencjalną (7) bryły wyko-
nanej z materiału liniowo sprężystego, po przedstawieniu odkształceń przez
naprężenia zgodnie z prawem Hooke’a i wykonaniu przypisanego całkowania
(rozwiązaniu prostego równania różniczkowego):

Φ =

1

2E

[(1 + ν)σ

ij

σ

ij

− νσ

kk

2

].

(7)

W dalszej części często będziemy stosować oznaczenia techniczne, dla któ-
rych (7) przyjmuje postać:

Φ =

1

2E

[σ

x

2

+σ

y

2

+σ

z

2

−2ν(σ

x

σ

y

+σ

y

σ

z

+σ

x

σ

z

)+2(1+ν)(τ

xy

2

+τ

xz

2

+τ

yz

2

)]

(8)

W przypadku brył niesprężystych duże zastosowanie ma twierdzenie Casti-
gliano o minimum energii dopełniającej δΠ

∗

= 0, które będzie przedmiotem

odrębnego opracowania.

2

Przedstawienie równań ZBTS w postaci macie-
rzowej

Przedstawimy macierzową postać podstawowych równań TS w celu wy-
prowadzenia metody elementów skończonych (MES) w zapisie macierzo-
wym.Przyjmujemy następujące oznaczenia: [macierz]ujmujemy w nawiasy
kawadratowe , |wektor| (macierz kolumnową lub jednowierszową) ujmujemy
w kreski pionowe. Zwracamy uwagę, że na potrzeby niniejszego przedstawie-
nia macierzowego wektory i macierze ZBTS, a także indeksowanie wyrażeń
ma inne znaczenie od tego, które stosowano w zapisie tensorowym.

3

2.1

Stan odkształcenia

Stan odkształcenia bryły odkształcalnej opisują związki geometryczne (Cau-
chy’ego), które w zapisie macierzowym przyjmą postać:

|ε| = [∂]|u|

(9)

gdzie wektor odkształceń |ε|, macierz operatorów różniczkowania [∂] oraz
wektor przemieszczeń |u| można zapisać w postaci:

|ε| =















ε

11

ε

22

ε

33

2ε

12

2ε

23

2ε

31















,

[∂] =












∂

1

0

0

0

∂

2

0

0

0

∂

3

∂

2

∂

1

0

0

∂

3

∂

2

∂

3

0

∂

1












,

|u| =









u

1

u

2

u

3









(10)

Symbolem ∂

i

oznaczono operator różniczkowania po współrzędnej x

i

: ∂

i

≡

∂/∂x

i

2.2

Stan naprężenia

Stan naprężenia opisują równanie równowagi (Naviera) , które w zapisie
macierzowym przyjmą postać:

[∂]

T

|σ| + |P | = |0|

(11)

gdzie T jest operatorem transpozycji macierzy ( w tym przypadku macierzy
operatorów różniczkowania (10)), a wektor naprężeń |σ| i wektor sił maso-
wych |P | można zapisać w postaci:

|σ| =















σ

11

σ

22

σ

33

σ

12

σ

23

σ

31















, |P | =









P

1

P

2

P

3









(12)

Statyczne warunki brzegowe na części A

q

powierzchni ciała A można zapisać

w postaci:

[S]|σ| − |q

υ

| = |0|,

(13)

gdzie macierz kosinusów kierunkowych [S] zapiszmy w postaci:

[S] =






α

υ1

0

0

α

υ2

0

α

υ3

0

α

υ2

0

α

υ1

α

υ3

0

0

0

α

υ3

0

α

υ2

α

υ1






,

(14)

4

a wektor obciążenia |q

υ

| oraz wektor normalny do powierzchni ciała |α

υ

| w

postaci:

|q

υ

| =









q

υ1

q

υ2

q

υ3









, |α

υ

| =









α

υ1

α

υ2

α

υ3









(15)

2.3

Związki konstytutywne

Równania fizyczne (konstytutywne)dla ciała liniowo-sprężystego opisują rów-
nania Hooke’a, które w zapisie macierzowym przyjmą postać:

|σ| = [E]|ε|

(16)

Macierz stałych materiałowych [E] można zapisać w postaci

[E] =












λ + 2µ

λ

λ

0

0

0

λ

λ + 2µ

λ

0

0

0

λ

λ

λ + 2µ

0

0

0

0

0

0

µ

0

0

0

0

0

0

µ

0

0

0

0

0

0

µ












,

(17)

gdzie współczynniki materiałowe Lamego można przeliczyć z modułu Youn-
ga E oraz współczynnika Poissone’a ν:
λ=νE/(1+ν)(1-2ν),
µ=G=E/2(1+ν).
Zwracamy uwagę, że tylko dwie z tych stałych materiałowych są niezależne,
a najczęściej stosowane są pary (E,ν) lub (λ,µ).
W szczególnych przypadkach płaskich stanów), mamy:
dla płaskiego stanu odkształcenia

[E] = E/(1 − ν

2

)






1

ν

0

ν

1

0

0

0

(1 − ν)/2






(18)

dla płaskiego stanu naprężenia

[E] = E/(1 + ν)(1 − 2ν)






1 − ν

ν

0

ν

1 − ν

0

0

0

(1 − 2ν)/2






(19)

2.4

Równanie prac wirtualnych

Równanie prac wirtualnych w zapisie macierzowym, otrzymamy po pod-
stawienia do (3) macierzowych równań ZBTS: (9),(11), (16). W rezultacie
uzyskamy(20):

Z

V

|σ|

T

δ|ε|dV =

Z

V

|P |

T

δ|u|dV +

Z

A

|q

υ

|

T

δ|u|dA,

(20)

5

lub po macierzowym przetransponowaniu powyższego równania, postać:

Z

V

δ|ε|

T

|σ|dV =

Z

V

δ|u|

T

|P |dV +

Z

A

δ|u|

T

|q

υ

|dA,

(21)

3

Rozwiązanie zagadnienia brzegowego teorii sprę-
żystości metodą Ritza

Załóżmy, że wektor przemieszczenia dowolnego punktu wewnątrz ciała moż-
na wyrazić przez liniową kombinację funkcji dopuszczalnych ϕ

i

(|x|) (funkcji

Ritza):

u(|x|) = ϕ

0

(|x|) + Σ

(i)

a

i

· ϕ

i

(|x|),

(22)

gdzie a

i

są stałymi Ritza, których bedziemy poszukiwać. Współczynniki

a

i

są stałymi w rozumieniu rozpatrywanej przestrzeni fizycznej |x

i

|, ale są

na razie nieznanymi zmiennymi aproksymacji Rizta. Często mówimy w ta-
kim przypadku, że ”stałe” są ”uzmiennione”. W rozumieniu aproksymacji
stałymi są natomiast pozostałe wielkości, w tym ”funkcje” przemieszczeń,
naprężeń i odkształceń. (22) można zapisać w postaci macierzowej

|u| = [N ]|a|,

(23)

gdzie macierz [N] jest macierzą kształtu.
Funkcje dopuszczalne są z góry przyjętymi dowolnymi funkcjami, ale speł-
niającymi warunki brzegowe.
Uwaga:
1) jeśli funkcja aproksymująca nie spełnia warunków brzegowych, to nie jest
dopuszczalna i nie może byc zastosowana.
2) dokładność aproksymacji dla konkretnego zagadnienia zależy od typu
funkcji Ritza. Najlepszym typem funkcji aproksymujących są ścisłe funkcje
kształtu, to znaczy wynikające ze ścisłego rozwiązania zagadnienia. Dla każ-
dego typu funkcji dopuszczalnych uzyskujemy jednak możliwe rozwiązanie
przybliżone, a jego dokładność można zwiększać poprzez zwiększanie liczby
funkcji dopuszczalnych w liniowej kombinacji (22).

3.1

Aproksymacja równania Cauchy’ego

Po podstawieniu do (9) funkcji aproksymacyjnej (23) równania geometrycz-
ne przyjmą postać:

|ε| = [∂][N ]|a| = [B]|a|,

(24)

gdzie wprowadzono macierz zgodności geometrycznej [B]

[B] = [∂][N ]

(25)

6

3.2

Aproksymacja prawa Hooke’a

Po podstawieniu do (16) funkcji aproksymacyjnej (23) związki fizyczne przyj-
mą postać:

|σ| = [E][B]|a|

(26)

3.3

Aproksymacja równania prac wirtualnych

Równanie prac wirtualnych (21) po wprowadzeniu aproksymacji (23) przyj-
mie postać:

Z

V

δ([B]|a|)

T

[E][B]|a|dV =

Z

V

δ([N ]|a|)

T

|P |dV +

Z

A

δ([N ]|a|)

T

|q

υ

|dA (27)

Ponieważ zmiennymi funkcji aproksymacyjnej są współczynniki |a|, nato-
miast [B] i [N] nie zależą od |a|, więc δ|a|

T

można wyłączyć przed znak

całek i w rezulatacie (27) możemy zapisać w postaci:

δ|a|

T

·

Z

V

[B]

T

[E][B]|a|dV = δ|a|

T

·

Z

V

[N ]

T

|P |dV +δ|a|

T

·

Z

A

[N ]

T

|q

υ

|dA (28)

3.4

Kanoniczne równanie Ritza

W celu wyznaczenia stałych |a| funkcji aproksymujących pole przemieszczeń
skorzystamy z warunku stacjonarności podług tych stałych, czyli zerowania
się pochodnych. W tym celu obie strony (29) różniczkujemy po δ|a|

T

.

, a warunek stacjonarności δL − δW możemy zapisać w postaci:

Z

V

[B]

T

[E][B]|a|dV =

Z

V

[N ]

T

|P |dV +

Z

A

[N ]

T

|q

υ

|dA

(29)

Po wprowadzeniu oznaczeń pomocniczych:

[k]

Z

V

[B]

T

[E][B]|a|dV,

(30)

|F

V

| =

Z

V

[N ]

T

|P |dV,

(31)

|F

A

| =

Z

A

[N ]

T

|q

υ

|dA,

(32)

gdzie: [k]- macierz sztywnosci, F

V

- wektor sił masowych, |F

A

|-wektor ob-

ciążeń.
Kanoniczne równanie metody Ritza można zapisać w zwartej formie:

[k]|a| = |F

V

| + |F

A

|.

(33)

7

4

Wprowadzenie do metody elementów skończo-
nych

4.1

Kanoniczne równanie Ritza fundamentalnym równaniem
MES

Kanoniczne równanie metody Ritza ((33) jest też fundamentalnym równa-
niem metody elementów skończonych, a definicje określone w rozdz. 3 defi-
niują podstawowe pojęcia tej metody.

4.2

Fundamentalne założenie MES

Istota konwencjonalnej metody elementów skończonych polega na specjal-
nym wyborze współczynników |a| w aproksymacji (23).
Mianowicie jako takie stałe wybiera się przemieszczenia węzłów elementu
skończonego u

(e)

. Zależności z poprzedniego rozdziału zapisuje się dla każ-

dego elementu skończonego (e):

|a| = u

(e)

=


















u

1

(1)

u

2

(1)

u

n

(1)

.....

u

1

(N )

u

2

(N )

u

n

(N )


















N xn

,

(34)

gdzie n- liczba składowych przemieszczenia punktu (w przestrzeni n=3); N-
liczba elementów skończonych. Równanie (23) można zapisać więc w postaci:

|u| = [N ]

(e)

|u|

(e)

=

(35)

4.3

Macierz kształtu elementu

Zwykle przyjmuje się, następującą aproksymację przemieszczeń wewnątrz
elementu u = u(x

1

, x

2

, x

3

):

|u| = [ϕ]|α|,

(36)

gdzie ϕ = ϕ((x

1

, x

2

, x

3

) są funkcjami przemieszczeń, a wektor|α| został

wprowadzony jako pomocnicza zmienna, która zostanie zredukowana w spo-
sób pokazany niżej.
Dla węzła (I) mamy:

|u(x

(I)
1

, x

(I)
2

, x

(I)
3

| = |u

(I)

| = [ϕ(x

(I)
1

, x

(I)
2

, x

(I)
3

]|α| = = [ϕ

(I)

] · |α|.

(37)

8

Zapisując podobnie dla pozostałych węzłów,otrzymamy:

u

(e)

=











u

(I)

u

(J )

.....

u

(N )











= [C]|α|,

gdzie

[C] =








ϕ

(I)

ϕ

(J )

.....

ϕ

(N )








(38)

Z powyższej zalezności wyznaczamy pomocniczy wektor |α|

|α| = [C]

−1

|u|

(e)

(39)

Podstawiajac uzyskaną zalezność do (36) otrzymamy

|u| = [ϕ][C]

−1

|u

(e)

| = [N ]|u

(e)

|,

(40)

gdzie [N] nazywa się macierzą kształtu, krtórą zapiszemy w postaci:

[N ] = [N

(1)

, N

(2)

, ..., N

(N )

],

gdzie

[N

(I)

] =






N

11

(I)

N

12

(I)

N

13

(I)

N

21

(I)

N

22

(I)

N

23

(I)

N

31

(I)

N

32

(I)

N

33

(I)






4.4

Dobór funkcji przemieszczeń

Funkcje przemieszczeń powinny być tak dobrane, by:
1) nie pozwalały na wytworzenie się stanu naprężenia w elemencie, jei prze-
mieszczenia węzłów powodują jedynie ruch elementu jako ciała sztywnego
2) w sytuacji zgodności przemieszczeń węzłów z warunkiem stałych odkształ-
ceń można było zrealizować te stałe odkształcenia
3) odkształcenia na granicy elementu były skończone.

4.5

Równania kanoniczne elementu

Kanoniczne równanie metody Ritza (33) w odniesieniu do elementu skoń-
czonego, przyjmie postać:

[k]

(e)

|u|

(e)

= |F

V

|

(e)

+ |F

A

|

(e)

,

(41)

gdzie

[k]

(e)

= [k

ij

]

(e)

=






k

11

(e)

...

k

1N

(e)

...

...

...

k

N 1

(e)

...

k

N N

(e)






(42)

9

|F

V

|

(e)

= |F

(e)

V i

| =









F

V 1

(e)

...

F

V N

(e)









(43)

|F

A

|

(e)

= |F

(e)

ai

| =









F

a1

(e)

...

F

aN

(e)









(44)

gdzie: N-liczba węzłów elementu skończonego.

4.6

Równania kanoniczne całej konstrukcji

Równanie kanoniczne dla całej konstrukcji zapiszemy w postaci:

[K]|u| = |F

V

| + |F

A

|,

(45)

gdzie składanie rozwiązań z elementów skończonych polega na szczególnym
”sumowaniu” macierzy sztywności [k]

(e)

, wektora sił masowych |F

V

(e)

| oraz

wektora sił powierzchniowych |F

A

(e)

|:

[K

ij

] =

(n)

X

(1)

[k

ij

]

(e)

,

(46)

|F

V i

| =

(n)

X

(1)

[F

V i

]

(e)

,

(47)

|F

Ai

| =

(n)

X

(1)

[F

Ai

]

(e)

,

(48)

Sumowanie dotyczy bloków macierzy sztywności i wektorów równoważników
węzłowych ”przyległych” do węzłów globalnych. Sposób składania rozwiązań
pokażemy na przykładach.

10