MODEL POWŁOKOWO-BELKOWY MES ANALIZY STATECZNOŚCI RAM PRZESTRZENNYCH O PRĘTACH CIENKOŚCIENNYCH OTWARTYCH

MODELOWANIE INŻYNIERSKIE ISSN 1896-771X 44, s. 131-138, Gliwice 2012

MODEL POWŁOKOWO-BELKOWY MES ANALIZY STATECZNOŚCI RAM PRZESTRZENNYCH

O PRĘTACH CIENKOŚCIENNYCH OTWARTYCH

SŁAWOMIR KOCZUBIEJ,CZESŁAW CICHOŃ

Katedra Informatyki Stosowanej, Politechnika Świętokrzyska e-mail: s.koczubiej@tu.kielce.pl

Streszczenie. W pracy zaproponowano model skończenie elementowy do analizy stateczności ram o prętach cienkościennych i dowolnym przekroju poprzecznym otwartym. W pracy potraktowano obszar węzła jako obiekt przestrzenny, modelowany elementami powłokowymi i dyskretyzację belkowymi elementami cienkościennymi połączono z przestrzenną dyskretyzacją węzłów. Do konsystentnego połączenia modelu powłokowego z modelem belkowym wykorzystano równania więzów dla sformułowania tzw. elementów przejściowych.

1. WSTĘP

Historia rozwoju metod obliczania belek i ram o prętach cienkościennych, z początku prostych – analitycznych (lecz dobrze fizycznie uzasadnionych), potem bardziej złożonych – wykorzystujących możliwości współczesnych metod obliczeniowych, przedstawiona jest w monografii Kujawy [6].

Spotykane dość powszechnie w praktyce inżynierskiej modelowanie ram przestrzennych za pomocą tylko belkowych elementów skończonych nie pozwala uwzględnić efektów lokalnych, takich jak wyboczenie środników ani opisać rzeczywistych warunków podparcia i rzeczywistego charakteru węzłów łączących pręty cienkościenne. Oznacza to, że ramy złożone z prętów cienkościennych powinny być w zasadzie modelowane z wykorzystaniem elementów skończonych powłokowych, co jednakże znacznie podnosi koszt obliczeń.

Racjonalnym rozwiązaniem przyjętym w pracy jest odróżnienie w modelowanej konstrukcji obszarów o charakterze przestrzennym od części o charakterze geometrycznie liniowym. W konstrukcji cienkościennej wyróżnia się: części (węzły ramy, podpory, miejsca przyłożenia obciążenia, miejsca z dodatkowymi wzmocnieniami konstrukcyjnymi), które potraktowano jako obiekty geometrycznie trójwymiarowe 3D i pozostałe części uważane za geometrycznie liniowe 1D. W modelu skończenie elementowym oznaczać to będzie dyskretyzację obszarów 3D przez elementy skończone powłokowe, a obszarów 1D przez belkowe elementy skończone cienkościenne. Podstawowym problemem, jaki w tym przypadku powstaje, jest złożenie tak różnie zdyskretyzowanych (to znaczy opisywanych różnymi modelami matematycznymi) części konstrukcji w jeden dyskretny system elementów skończonych, opisany układem równań równowagi MES.

Naturalnym sposobem postępowania jest skorzystanie z warunku ciągłości przemieszczeń translacyjnych na ściankach, wspólnych dla powłoki i belek cienkościennych

), ( )

( ^(b)

(s) x u x

u  (1)

gdzie u^(s)(x) jest wektorem przemieszczeń węzła elementu powłokowego, a u^(b)(x) jest wektorem przemieszczeń węzła elementu belkowego.

W pracy równanie (1) traktowane jest jako równanie więzów i włączone do równania równowagi MES poprzez odpowiednio zdefiniowaną funkcję kary. Konsekwencją takiego postępowania jest konieczność wprowadzenia tzw. elementów przejściowych pomiędzy węzłami powłoki a węzłem cienkościennego elementu skończonego, przedstawionych na rys. 1. Jest to metoda ogólna, możliwa do zastosowania przy rozwiązaniu zarówno zagadnień liniowych jak i nieliniowych. Oryginalność postępowania będzie w tym przypadku polegać na sformułowaniu równań dla elementu przejściowego w ramach przyjętej teorii i jego implementacji w modelu MES.

Rys. 1. Przejściowy element skończony

Innym sposobem rozważanym w [4, 5] jest wykorzystanie przestrzennego elementu węzłowego, który jest obiektem przestrzennym, zdyskretyzowanym powłokowymi elementami skończonymi, w którym dokonywana jest kondensacja statyczna i transformacja redukująca powłokowe stopnie swobody węzła do belkowych cienkościennych stopni swobody zlokalizowanych na ściankach (przekrojach poprzecznych pręta wspólnych dla modelu belkowego i powłokowego), łączących węzeł z prętami ramy.

Przyjęto następujące podstawowe założenia:

 obciążenie konstrukcji jest jednoparametrowe i stanowią go uogólnione siły statyczne, o ustalonej konfiguracji,

 materiał konstrukcji jest jednorodnym i izotropowym continuum materialnym,

 rozważane continuum materialne ma własności fizyczne liniowo-sprężyste,

 odkształcenia są nieskończenie małe, lecz możliwe są duże przemieszczenia i duże gradienty przemieszczeń,

 w metodzie elementów skończonych konstrukcje dyskretyzowano elementami belkowymi cienkościennymi, wykorzystującymi założenia teorii Własowa [3], oraz elementami płytowymi Reissnera–Mindlina z dodatkowym szóstym stopniem swobody,

 korzystano z metody elementów skończonych w sformułowaniu przemieszczeniowym,

 w konstrukcji wszystkie połączenia są spawane.

2. RÓWNANIE WARIACYJNE RÓWNOWAGI

Rozważano ruch ciała w różnych stanach równowagi (rys. 2). Stan C⁰ jest stanem początkowym równowagi, stan C^t jest aktualnym stanem równowagi (na początku przyrostu), a stan C^tt jest sąsiednim stanem równowagi (na końcu przyrostu).

Rys. 2. Trzy stany równowagi ciała: a) początkowy C⁰, b) aktualny C^t i przyrostowy C^tt V0 i V^t są odpowiednio początkową i aktualną objętością ciała, S jest brzegiem ciała z wyróżnionymi częściami S i _u S , na których są zdefiniowane kinematyczne i statyczne _t warunki brzegowe, g jest wektorem intensywności sił objętościowych, a t jest wektorem intensywności sił powierzchniowych. W końcu wektor x określa położenie punktu P ⁰ w konfiguracji C⁰, wektor X lokalizuje ten punkt (P w konfiguracji ^t) C^t, a wektor

x X

u  jest wektorem przemieszczeń.

Zadaniem autorów pracy jest wyznaczenie stanu ciała (przemieszczeń, odkształceń, naprężeń) w konfiguracji sąsiedniej C^tt, przyjmując, że stan ciała w konfiguracji aktualnej

Ct jest znany.

Ciało w stanie C^tt jest w równowadze, jeżeli spełnione jest równanie wariacyjne

(t) ,

(c)

z w

w W W

W  

   (2)

w którym, mając na względzie przyszłą dyskretyzację skończenie elementową, wyróżniono pracę sił wewnętrznych dla continuum W_w^(c) (dyskretyzowanego elementami skończonymi powłokowymi i/lub belkowymi cienkościennymi) oraz pracę wirtualną sił wewnętrznych dla elementów skończonych przejściowych W_w^(t). Przez W_z oznaczono pracę wirtualną sił zewnętrznych (obciążeń).

2.1. Praca wirtualna W_w^(c)

Pracę wirtualną W_w^(c) w stanie C^tt wyraża wzór , d

0

0 T (c) ^



V

w V

W ε σ

 (3)

w którym ε jest wektorem odkształceń Greena–Lagrange'a i σ jest energetycznie z nim sprzężonym wektorem naprężeń Pioli–Kirchhoffa II rodzaju, a indeks górny T oznacza transpozycję.

Równanie (3) nie może być wprost wykorzystane, ponieważ stan ciała w konfiguracji

t

Ct^ jest nieznany. Procedura postępowania polega na przyjęciu zlinearyzowanego przyrostowego równania konstytutywnego oraz zachowania pod znakiem wariacji liniowej

części przyrostu odkształcenia ε . Oznaczając przez σ^t naprężenie w chwili t , a przez σ^{t  t} w chwili tt, przyrost naprężenia pomiędzy chwilami t i tt wynosi

, e D ε D σ σ

σ     

 ^{tt t} (4)

gdzie dla materiału liniowo sprężystego D jest macierzą stałych materiałowych, a ε jest wektorem przyrostów odkształceń. Wektor ε możemy z kolei przedstawić w postaci sumy części liniowej e i części kwadratowej η ze względu na przyrosty przemieszczeń

η. e ε 

 (5)

Wykorzystanie w (3) związków (4) i (5) oraz ε^tt ε prowadzi do zależności

 

 

 

0 0

0

T 0 T 0

T 0

(c) ^







V

t

V

t

V

w V V V

W  e D e  η σ  e σ

 (6)

Wektor przyrostów odkształceń liniowych można formalnie zapisać w postaci ,

u L e 

 (7)

w którym wektor przyrostów przemieszczeń wynosi uu^tt u^t, a macierz L jest macierzą operatorów różniczkowych. Oznaczając w dalszym ciągu przez N wektor funkcji kształtu i przez q wektor stopni swobody dla elementu skończonego, interpolację u i u wyrażają wzory

,

, u N q

Nq

u    (8)

gdzie u jest wektorem uogólnionych funkcji przemieszczeń, u jest wektorem przyrostów uogólnionych funkcji przemieszczeń, N jest macierzą funkcji kształtu, q jest wektorem stopni swobody elementu a q jest wektorem przyrostów stopni swobody elementu.

Podstawiając (8) do (7), otrzymano zależność

. q LN q B

e   

 _L (9)

Wykorzystując powyższe zależności we wzorze (6) na pracę wirtualną W_w^(c) dla elementu skończonego otrzymano równania:

1.

 

 

 

0 0

0 T (c)

(c) (c)

T (c) T 0





V

L L

L L

V

V

V q k q k B DB

e D

e 

 (10)

gdzie k jest przyrostową macierzą sztywności; ^(c)_L

2.

 

 

 

0 0

0 T (c)

(c) (c)

T (c) T 0





V

N t N V

t V q k q k B T B V

σ

η  _{ }

 (11)

z macierzą początkowych naprężeń (macierzą geometryczną) k_^(c), w której naprężenia Pioli- Kirchhoffa II rodzaju są przedstawione w formie macierzowej T , a postać macierzy ^t zróżniczkowanych funkcji kształtu B^(c)_N wynika z aproksymacji skończenie elementowej wektora przyrostów odkształceń nieliniowych η;

3.

 

 

 

0 0

T 0 (c) (c)

T (c) T 0





V

t L V

t V q f f B σ V

σ

e 

 (12)

gdzie f^(c) jest wektorem sił wewnętrznych.

2.2. Praca wirtualna W_w^(t)

Równanie więzów (1) zapisane w postaci

(b),

(s) u

u

r  (13)

może być wprowadzone do dyskretyzacji skończenie elementowej na wiele sposobów [1, 7].

W pracy wykorzystano metodę funkcji kary, definiując pracę wirtualną W_w^(t) w stanie C^tt, przedstawiającą wkład równania więzów do opisu ogólnej sztywności konstrukcji, zależnością

  

 



(t)   rr rr   r rr r

W_w k k (14)

gdzie r jest wektorem równań więzów w chwili t , a k jest parametrem kary. Równanie więzów wyraża warunek zgodności z wagą k przemieszczeń uogólnionych dla dowolnego punktu, wspólnego belki i powłoki. Wektor przyrostu funkcji kary r można przedstawić w postaci sumy części liniowej ρ i kwadratowej υ względem przyrostów przemieszczeń uogólnionych

, υ ρ r 

 (15)

co po linearyzacji prowadzi do zależności analogicznej do (6)

 

 

 

(t)   ρ ρ υ r ρ r

W_w k k k (16)

Podobnie jak to miało miejsce przy rozważaniu pracy wirtualnej W_w^(c) również i obecnie można zdefiniować pojedynczy element skończony przejściowy z wektorem stopni swobody q pomiędzy węzłem belki cienkościennej i węzłem powłoki, co pozwala formalnie zapisać (t)

składowe z równania (16) w formie:

1. k









 

2. k









3. ^k









Zależności na r, B i ^(t)_L B^(t)_ są wyprowadzone w [4].

2.3. Praca wirtualna W_z

Pracę wirtualną obciążenia w stanie C^tt wyraża zależność

   







0 0

T 0 T 0

d d

V S

t t t

t z

t

S V

W  u g  u t

 (20)

która po wykorzystaniu (8) przyjmie postać

   

0 0

0 T

0 T T















   





t

t S

t t t

t t S

t t t

t g t t t t t g

z V t S

W  q p p p N g p N

 (21)

gdzie p^t_g^t jest wektorem intensywności węzłowych sił objętościowych, a p^t_t^t jest wektorem intensywności węzłowych sił powierzchniowych.

3. PRZYROSTOWE RÓWNANIE RÓWNOWAGI MES

Przyrostowe równanie równowagi dla elementu skończonego continuum (powłoki, belki cienkościennej) otrzymuje się z równania wariacyjnego (2), korzystając z zależności (10-12) i (20-21)



 











t t g

L 

 _ (22)

Identyczność ta musi być spełniona dla dowolnej wariacji przyrostu przemieszczeń, co prowadzi do równania równowagi dla elementu skończonego w formie





t g

L  (23)

Sumę macierzy sztywności

, k k

k_T  _L  (24)

nazywa się macierzą sztywności stycznej (macierzą styczną).

Agregując elementy skończone w jeden układ dyskretny, otrzymano przyrostowy układ równań równowagi MES dla konstrukcji w postaci

, , _g^{t t} _t^{t t}

T







 







Q P F P P P

K (25)

w którym dużymi literami oznaczono globalną macierz styczną i globalne wektory dla całej konstrukcji.

W przypadku dyskretyzacji konstrukcji elementami skończonymi powłokowymi, belkowymi i przejściowymi udział tych ostatnich w macierzy globalnej sztywności stycznej i w wektorze globalnym sił wewnętrznych w (25) jest dokonywany przez odpowiednią, dodatkową procedurę agregacji.

4. PRZYKŁAD

Opisany model obliczeniowy MES został zaimplementowany do autorskiego programu AmFEM działającego w systemie MATLAB. Nieliniowy układ równań równowagi rozwiązywano iteracyjnie metodą Newtona–Raphsona, sterując procesem obliczeń za pomocą parametru obciążenia lub wybranego przemieszczenia. Przy opracowywaniu procedur do obliczania macierzy sztywności stycznej i wektora sił wewnętrznych dla elementu skończonego cienkościennego, które obliczano dokładnie w ramach przyjętej teorii geometrycznie nieliniowej, korzystano z obliczeń symbolicznych oferowanych przez system MATLAB i język programowania PERL, wspierający operacje z wykorzystaniem wyrażeń regularnych.

Tematem przykładu jest analiza nieliniowa geometrycznie ramy portalowej o przekroju dwuteowym i obciążeniu pokazanym na rys. 2a. W tym przykładzie zmieniano dodatkowo

parametr a , czyli wielkość części powłokowej w granicach a1 ,5 5h, gdzie h jest wysokością dwuteownika.

Obliczenia dla modelu 3D wykonano w systemie ABAQUS z użyciem elementów S8R6.

Model 3D był dyskretyzowany 1304 elementami powłokowymi i miał N 26094 stopni swobody. W przypadku modeli 1D/3D zostało użytych 87 elementów belkowych i od 136 do 376 elementów powłokowych w zależności od wielkości części przestrzennej

) 15426 6066

(N   . Części belkowe i powłokowe były łączone za pomocą 17 elementów przejściowych dla każdej ze ścianek.

Na rys. 3b i 3c zebrano rzuty ścieżki stanu równowagi na płaszczyznę









w węźle w . ₁

Rys. 3. Rama portalowa: a) geometria i dane materiałowe, b) przemieszczenie u w węźle _x w₁ , c) przemieszczenie i u w węźle _z w ₁

W obu przypadkach krzywe ścieżek stanów równowagi obliczone dla modeli 1D/3D i 3D nieznacznie różnią się. Różnice te mogą wynikać z różnego stopnia nieliniowości zastosowanych elementów. Zastosowane w pracy płaskie elementy powłokowe mogą być stosowane do analizy zadań o małym i średnim stopniu nieliniowości. Natomiast belkowe cienkościenne elementy skończone charakteryzują się możliwością rozwiązywania zagadnień silnie nieliniowych. Zatem zwiększenie części ramy modelowanej elementami belkowymi spowodowało prawdopodobnie poprawę wyników. Należy jednakże pamiętać, że części modelowane elementami powłokowymi nie powinny mieć mniejszych rozmiarów niż

h

a1,5 [2]. Oczywiście w przypadku analizy liniowej wszystkie zaprezentowane modele 3D i 1D/3D dają takie same wyniki.

5. PODSUMOWANIE

W pracy zaproponowano mieszany, powłokowo-belkowy model skończenie elementowy do analizy ram o prętach cienkościennych z różnymi geometriami węzłów, podpór i różnym sposobem przyłożenia obciążeń. Ograniczenie dyskretyzacji za pomocą powłokowych elementów skończonych tylko do pewnych charakterystycznych części konstrukcji znacznie

upraszcza obliczenia w porównaniu z pełną analizą trójwymiarową, zapewniając jednocześnie dobrą jakość rozwiązania.

Ponadto użycie nowoczesnych narzędzi programistycznych (przede wszystkim wyrażeń regularnych i języka programowania PERL) pozwoliło na efektywne wyprowadzenie silnej nieliniowo macierzy stycznej i wektora sił wewnętrznych cienkościennego belkowego elementu skończonego.

LITERATURA

1. Chavan K. S., Wriggers P.: Consistent coupling of beam and shell models for thermo- elastic analysis. „International Journal for Numerical Methods In Engineering” 2004, 59, 14, p. 1861–1878.

2. Cichoń C., Koczubiej S.: Konsystentny model MES dla ram przestrzennych o prętach cienkościennych. „Czasopismo Techniczne” 2008, 21, 1-B/2008, s. 3–20.

3. Kim M. Y., Chang S. P., Kim S. B.: Spatial stability analysis of thin-walled space frames.

„International Journal for Numerical Method in Engineering” 1996, 39, s. 499–525.

4. Koczubiej S.: Model powłokowo-belkowy MES w analizie statycznej i stateczności konstrukcji o prętach cienkościennych otwartych. Praca doktorska. Kielce: Politechnika Świętokrzyska, 2011 (www.tu.kielce.pl/~sk/?notka).

5. Koczubiej S., Cichoń C.: Shell-beam model of thin-walled space structures for geometrical nonlinear analysis. In: CMM-2011 Computer Methods in Mechanics.

Warszawa 2011, p. 255-256.

6. Kujawa M.: Statyka i analiza wrażliwości rusztów zbudowanych z prętów cienkościennych : analiza teoretyczna i badania doświadczalne. Gdańsk: Wyd. Pol.

Gdańskiej, 2009.

7. Wagner W., Gruttmann. F.: Modeling of shell-beam transitions in the presence of finite rotations. „Computer Assisted Mechanics and Engineering Sciences” 2002, 9, 3, p. 405–

418.

SHELL-BEAM FE MODEL IN STABILITY ANALYSIS OF THIN-WALLED STRUCTURES

WITH OPEN CROSS SECTION

Summary. In the paper a finite element model in analysis of thin-wall structures is presented. The standard discretization using beam thin-walled elements is connected with the space discretization of the some parts of the frame. Transition elements are used for consistent coupling between shell (3D) and beam (1D) elements