COSMOS v procesním inženýrství a strojnictví

COSMOS v procesním inenýrství a strojnictví (èást I. metody dimenzování aparátù) revize: 4/1997

R. itný

tel. 311 59 51, E-mail: [email protected] Katedra strojù a zaøízení pro chemický, potravináøský a spotøební prùmysl, Fakulta strojní, ÈVUT, Technická 4, Praha 6

COSMOS v procesním inenýrství, I.MDA 1997

2

Literatura pouitá (doporuèená tuènì): [1]

Bittnar,Z.-ejnoha,J.: Numerické metody mechaniky 1, ÈVUT Praha 1992

[2]

Bittnar,Z.-ejnoha,J.: Numerické metody mechaniky 2, ÈVUT Praha 1992

[3]

Timoenko, S.P.-Goodier,J.N.: Theory of elasticity, 3-rd edt., McGraw-Hill, N.Y., 1970

[4]

Apetaur,M. et al: Výpoètové metody ve stavbì motorových vozidel, ÈVUT Praha, 1986

[5]

Gallagher,R.H.: Finite element analysis, Prentice Hall, N.J. 1975 (ruský pøedklad Mir, Moskva 1984)

[6]

Jao,I.-Lukavský,J.: Konstrukce aparátù, ÈVUT Praha, 1991

[7]

Bathe,K.J.-Wilson,E.L.: Numerical methods in finite element analysis, Prentice Hall, Englewood Cliffs, N.J., 1976

[8]

Bittnar,Z.-Øeøicha,P.: Metoda koneèných prvkù v dynamice konstrukcí, SNTL, Praha 1981

[9]

Seide,P.: Theoretical manual for COSMOS/M, SRAC, Santa Monica, 1992

[10]

Zienkiewicz,O.C.: The finite element method in engineering science, McGraw-Hill, London 1971

[11]

Brebbia.C.A.: The boundary element method for engineers, Pentech Press, London, 1978

[12]

Køupka,V.-Schneider,P.: Stavba chemických zaøízení I., Skoøepiny tlakových nádob a nádrí, VUT Brno, 1986

[13]

Vejvoda,S.-Vlk,M.: Stavba chemických zaøízení II., Pevnost a ivotnost tlakových nádob, VUT Brno,1982

[14]

Schneider,P.: Základy konstruování chemických zaøízení, VUT Brno, 1984

[15]

Rejent,B.: ivotnost potrubních systémù, ÈVUT Praha, 1989

[16]

ejnoha,J.-Kufner,V.: Prunost, pevnost, plasticita III., ÈVUT Praha, 1990

[17]

Valenta,J.-Nìmec.,J.-Ulrych,Z.: Novodobé metody výpoètu tuhosti a pevnosti ve strojírenství, SNTL, Praha


3

1975 [18]

ASME-Boiler and pressure vessel code, sect. VIII, div.2, 1989

[19]

Schneider,P.-Vykutil,J.: Stavba chemických zaøízení II a mikropoèítaèové aplikace MKP ve statice rotaèních skoøepin, VÚT Brno, 1990

[20]

Owen,D.R.J.-Hinton,E.: Finite elements in plasticity, Pineridge Press, Swansea 1980

[21]

Koláø,V.-Leitner,F.-Zlámal,M.-eníek,A.: Výpoèet ploných a prostorových konstrukcí metodou koneèných prvkù, SNTL Praha, 1972

[22]

esták,J.-Rieger,F.: Pøenosové jevy I., hydrodynamika a sdílení tepla, ÈVUT Praha 1972

[23]

Okrouhlík,M.: Aplikovaná mechanika kontinua I. , ÈVUT Praha 1986

[24]

Okrouhlík,M.: Aplikovaná mechanika kontinua II. , ÈVUT Praha 1989

[25]

Okrouhlík,M.: Aplikovaná mechanika kontinua III., ÈVUT Praha 1990

[26]

Zienkiewicz,O.C.: The finite element method, 3rd ed., N.Y. 1977

[27]

Norrie,D.H.-Vries,G.: An introduction to finite element analysis, Academic Press London 1978 (ruský pøeklad Mir, Moskva 1981)

[28]

Oden,J.T. edt.: Computational methods in nonlinear mechanics, Proceedings of the TICOM second int. conf., North Holland Publ., Amsterdam 1980

[29]

Nìmec, J.: Tuhost a pevnost ocelových èástí, ÈSAV Praha 1963

[30]

ubrt, L.: Teorie desek a skoøepin, ÈVUT, FS Praha 1987

[31]

Anonym: Prunost a pevnost I., ÈVUT FS, Praha 1963

[32]

Anonym: Prunost a pevnost II., ÈVUT FS, Praha 1963

[33]

Geostar, COSMOS/M, Vol. I.,II.,III., uivatelská pøíruèka, Structural Research and Analysis Corporation, St. Monica, 1992


4

[34]

Køupka, V.: Výpoèet válcových tenkostìnných kovových nádob a potrubí, SNTL Praha, 1967

[35]

Press,W.H.-Flammery,B.P.-Teukolsky,S.A.-Vetterling,W.T.: Numerical recipes, Cambridge Univ. Press, 1986

[36]

Gere,J.M.-Timoshenko,S.P.: Mechanics of materials, Wadsworth Int., Belmont, 1985

[37]

esták,J. a kol: Transportní a termodynamická data pro výpoèet aparátù a strojního zaøízení, ÈVUT FSI, Praha, 1981

[38]

Heømanský, B.: Termomechanika jaderných reaktorù, Academia Praha 1986

[39]

Pawlak,T.P.-Yunus,S.M.-Cook,R.D.: Solid elements with rotational degrees of freedom, Part II Tetrahedron elements, International Journal for Numerical Methods in Engineering, Vol. 31, 593-610 (1991)

[40]

Tuèek,A.-Sojka,A.: Zvýení sismické odolnosti potrubních systémù jejich ladìním. Èerpadla, potrubí, armatury, è.2-3, str. 41-48 (1989)

[41]

Smith,P.R.-Van Laan,T.J.: Piping and pipe support systems, McGraw Hill, New York, 1965


5

Seznam nejpouívanìjích symbolù: B

matice deformací tvarových funkcí (deformace=B x uzlová posunutí)

[1/m]

C

matice tlumení

[kg/s]

D,d

prùmìr válce (vnìjí, vnitøní)

[m]

D

matice elastických konstant Hookova zákona (41)

[Pa]

t e

tenzor deformace

E

Youngùv modul prunosti

[Pa]

Et Es Ei Ee EL Es

Teèný modul prunosti

[Pa]

Seèný modul prunosti

[Pa]

F

osamìlá síla (nebo zobecnìná osamìlá síla, napø. i momenty)

Fω

Fourierova transformace síly

f

objemová síla (napø. zrychlení x hustota)

G

smykový modul prunosti

gk

váhové koeficienty gaussovy integraèní formule

[-]

Deformaèní energie (i-internal vnitøní)

[J]

Potenciální energie vnìjích sil (e-external vnìjí)

[J]

Celková potenciální energie (L-Lagrange)

[J]

Energie potøebná na poruení soudrnosti materiálu

[J]

I σ ,II σ ,III σ

první, druhý a tøetí invariant tenzoru napìtí

[N] [N.s] [

N.m-3 ] [Pa] [-]

[ Pa, Pa

2

, Pa 3 ]

H

modul plasticity

J

moment setrvaènosti prùøezu

k

Kármánovo èíslo (zvýení poddajnosti kolena)

K

matice tuhosti

[N/m]

Kt Kσ Kc

teèná matice tuhosti (pøírùstek síly / pøírùstek pusunutí)

[N/m]

matice geometrické tuhosti

[N/m]

L

délka

[m]

M

matice hmot

[kg]

souèinitel intenzity napìtí (lomová mechanika)

M, M o , M k

moment, ohybový moment, kroutící moment

[Pa] 4

[m ] [-]

[

N.m-3 / 2 ]

[N.m]

N(x)

matice tvarových funkcí

[-]

p

tlak

[Pa]

p

vektor "modálních" posuvù (vektor koeficientù lineární kombinace vlastních tvarù)

[m]

q

vektor zobecnìných posunutí uzlových bodù (posuvy + ev. natoèení)

[m]

q

hustota tepelného toku (sloka kolmá k rozhraní)

qω

Fourierova transformace posuvù

Q

matice vlastních vektorù

[m]

r

radiální souøadnice, polomìr

[m]

-2

[ W.m ] [m.s]


S

6

2

[m ]

plocha prùøezu

t t S F ,S q

matice spektrálních výkonových hustot prùbìhù sil a posunutí

s

tlouka skoøepiny

t

èas

T

teplota

[K]

Tf

teplota fázového rozhraní

[K]

T

posouvající síly (nosníky, skoøepiny)

[N]

u,v,w

posunutí ve smìrech x,y,z

[m]

W

hustota deformaèní energie

x,y,z

souøadnice kartézského souøadného systému

α souèinitel teplotní roztanosti β korekèní souèinitel (vliv rozloení smykových napìtí po prùøezu) β pøevrácená hodnota délky "vlnových" poruch membránové napjatosti γ smyková deformace ε deformace ε e ,ε p ,ε c elastická, plastická a creepová deformace κ zpevnìní materiálu λ tepelná vodivost λ násobek referenèního zatíení µ Poissonova konstanta (pøíèná kontrakce) ρ hustota σ napìtí napìtí na mezi kluzu σK dovolené napìtí σD Misesova intenzita napìtí σM hlavní napìtí ve smìru i σi ϕ úhlová souøadnice ψ natoèení prùøezu ω úhlová frekvence Zkratky DOF

Degree Of Freedom (stupeò volnosti)

PSD

Power Spectral Density (Výkonová spektrální hustota)

CQC

Complete Quadratic Combination

NRL

Naval Research Laboratory

SRSS

Square Root Sum of Squares

RMS

Root Mean Square

[

N 2 s2 , m2 s2 ] [m] [s]

[

J.m-3 ] [m] [1/K] [-] [1/m] [-] [-] [-] -3

J.m ] -1 [ W.m .K ] [

-1

[-] [-] -3

[ kg.m ] [Pa] [Pa] [Pa] [Pa] [Pa] [rad] [rad] [rad/s]


7

Obsah 1.

Pøehled aplikací MKP: lineární statika, velké deformace, velká pøetvoøení, plasticita, únava, dynamika, stabilita, pøenos tepla, proudìní.

2.

Principy MKP: odvození tyèového prvku a napsání programu pro seriovì øazené pruiny, základní postup øeení (návrhový diagram-COSMOS), geometrické modelování (entity), zpùsoby vytváøení sítì (makroprvky, automatické generátory na oblastech), základní pøíkazy, pøedchozí pøíklad (tyèové prvky) zapsaný v COSMOSu. Napìtí a deformace: zopakování základních pojmù teorie prunosti, tenzory napìtí (hlavní napìtí a ekvivalentní napìtí), deformace

3.

Lineární statika: Hookeùv zákon (v maticové notaci) a vliv teplotního pole. Mìrná energie deformace, variaèní principy (Lagrange, Castigliano) a jejich pouití v MKP (tvarové funkce). 3.1 3.1.1

Lineární prvky (nosníky, potrubní prvky) Analytická øeení: køivé pruty (Castigliano), tah, ohyb, krut pøímých nosníkù a trubek, zvlátní pøípady: koleno, prstence.

3.1.2

Bernoulliho a Timoenkùv nosník, odvození matic tuhosti a vektoru zatíení

3.1.3

TRUSS, BEAM, PIPE, ELBOW

3.2 3.2.1

Dvourozmìrné prvky (skoøepiny, prstencové prvky) Analytická øeení: Laplaceova rovnice pro membrány (odvození), rotaènì symetrické skoøepiny (diferenciální rovnice prùhybu), teplotní napìtí

3.2.2

Tvarové funkce (izoparametrické prvky) a Gaussova integrace matice tuhosti, rotaènì symetrické skoøepiny s nesymetrickým zaíením (Fourierùv rozvoj)

3.2.3 3.3 4.

PLANE2D, SHELLAX, SHELL4 Tøídimenzionální prvky: SOLID

Nelineární statika koncept mezních stavù 4.1

Velké deformace a posunutí (pøíklad víka nádrí), odvození matice tuhosti pro tyèové prvky

4.2

Elasticko-plastické chování materiálu: plocha plasticity, asociativní teorie plasticity, konstitutivní rovnice (elasticko-plastická matice tuhosti), zpevnìní materiálu (kinematické a izotropní), pøizpùsobení konstrukce.

5.

4.3

Teèení materiálu (creep); modely teèení pouívané programem COSMOS

4.4

Pøírùstková metoda øeení (normální a modifikovaná metoda Newtonova Raphsonova)

Vyhodnocení napjatosti 5.1

Kategorizace napìtí dle ASME CODu

5.2

Únava materiálu: lomová mechanika a ivotnost konstrukce dle ASME (Minerovo pravidlo)


6.

8

5.3

Nìkterá speciální kriteria poruení materiálu (Laminátové konstrukce, kriterium Tsai Wu )

5.4

Optimalizace konstrukcí (pøíhradové konstrukce s minimální vahou, modul OPTIMAL)

Stabilita 6.1

Analytická øeení: stabilita tlaèených nosníkù (Euler, Tetmayer, ...), stabilita válcových skoøepin (tlaèený válec, válec s vnìjím pøetlakem)

6.2

Koneènì-prvková formulace pro pøípad lineární stability a metody øeení vlastní úlohy (Jacobiho metoda, inverzní iterace, iterace podprostoru, Lanczosova metoda)

7.

Dynamika 7.1

Vlastní tvary a vlastní frekvence kmitání konstrukcí

7.2

Vynucené kmitání a seismicita (modální superpozice, Duhamelùv integrál, metody pøímé integrace - Newmarkova a Wilsonova, náhodné buzení se zadaným prùbìhem výkonové spektrální hustoty /seismicita/)

8.

Teplotní pole-matice tepelných kapacit a matice vodivostí, odvození pro rovinný izoparametrický prvek. Výpoèty nestacionárních teplotních polí pøímou integrací.

Pøíloha 1.

Základní pøíkazy programu COSMOS/M


1.

9

Pøehled aplikací MKP Pøíklady komerèních programových systémù pro metodu koneèných prvkù:

NASTRAN, ANSYS, SYSTUS, MARC, ABAQUS, COSMOS/M, PATRAN, ADINA, ROCKFIELD, GT STRUDL, SAP IV, TPS-10 MKP zahrnutá pøímo do systémù pro konstruování ACAD/COSMOS, CADSS5, KATIA, INTERGRAPH, I-DEAS Základním úkolem vech tìchto systémù je pøedevím spoèítat pøesný prùbìh napìtí v libovolnì komplikovaných a zatíených strojích èi aparátech. V metodì koneèných prvkù (MKP) se zpravidla nejprve poèítá pøetvoøení zatíené konstrukce a teprve následnì se zjistí stav napjatosti. Pro dimenzování aparátù existuje øada smìrnic a norem (napø. ASME code [18]), které vak nejsou alternativou MKP; tyto pøedpisy vìtinou pøedpokládají, e konstrukce byla analyzována MKP a hodnotí teprve její výsledky (napø. tím, e stanoví souèinitele bezpeènosti).

Lineární statika Do této skupiny patøí cca 60% vech aplikací (pevnostní kontroly pøíhradových konstukcí, potrubních systémù, vìtiny souèástí zpracovatelských strojù, do jisté míry i tlakových nádob a skoøepinových konstrukcí). Princip lineární analýzy konstrukce lze znázornit na jednorozmìrném pøípadu pruinky zatíené silou F. Nejprve se urèí tuhost pruiny K a z ní posunutí koncového bodu (uzlu) u:

Obecný trojrozmìrný pøípad: Konstrukci nahradíme soustavou elastických prvkù (napø.pruinek), propojených v tzv. uzlových bodech, které budou souèasnì pùsobitìm vnìjích sil

r F

(mohou to být osamìlé síly,

zatíení tlakem, odstøedivé síly, zatíení vlastní vahou, ale i síly vyvolané dilatací, tedy napø. rozdílem montáních a provozních teplot, nebo rozdílem teplot na vnitøním a vnìjím povrchu nádoby ap). Vztah mezi tìmito silami a posuny uzlù je dán soustavou lineárních algebraických rovnic (poèet rovnic je roven poètu N posunutí ve vech uzlových bodech, nìkdy se hovoøí o poètu stupòù volnosti konstrukce)

t r r K.u = F

kde

t K

je ètvercová matice tuhosti,

pùsobících v uzlových bodech.

r u -vektor

posunutí uzlových bodù a

r F -vektor

vnìjích sil,


10

t K jsou tuhosti "náhradních" pruinek a pokud bychom jimi propojili kadý uzel s t kadým, byly by vechny prvky K nenulové. V metodì koneèných prvkù (dále jen MKP) se ovem propojují vdy Prvky matice tuhosti

jenom sousední uzly, take napø. i-tý øádek matice tuhosti obsahuje pouze tuhosti "pruinek", které spojují uzel èíslo

i

t K

s jeho bezprostøedními sousedy. Matice

je potom øídká (vìtina jejích prvkù jsou nuly) a symetrická. Kromì

pruinek se v MKP nahrazuje skuteèná konstrukce napø. malými elastickými trojúhelníkovými èi ètyøúhelníkovými destièkami nebo krychlièkami, jejich elastické chování (tuhost) dokáeme dobøe aproximovat (jsou-li tyto elementy dostateènì malé, mùeme napø. pøedpokládat, e v kadém z nich je prùbìh posunutí x,y,z a deformace

ε

u

lineární funkcí souøadnic

je konstantní).

První fáze MKP tedy spoèívá v sestavení matice tuhosti (u lineární statiky závisí pouze na geometrii souèásti a modulech prunosti materiálu) a následném øeení soustavy (1) napø. Gaussovou eliminací. Výsledkem je vektor posunutí uzlových bodù.

x

Posunutí v obecném bodì interpolaèní bázovou funkci

N i (x)

získáme interpolací uzlových posunutí

ui .

Kadému uzlu pøiøadíme

(èasto se jí øíká tvarová funkce), která je rovna 1 v uzlu èíslo

i

a ve vech

ostatních uzlech je nulová. Potom platí: N

u(x) = ∑ N i (x) ui i=1

V MKP se vyetøovaná souèást rozkládá na elementy (napø. ètyøúhelníky u skoøepinových konstrukcí) tak, aby uzly leely na rozhraní elementù. Tvarové funkce

N i (x)

se definují v kadém elementu zvlá napø. jako lineární nebo

kvadratický polynom, plnì urèený polohou uzlových bodù. Pøíklady tvarových funkcí jsou uvedeny na následujícím obrázku Pozn.: Pojem tvarových funkcí je dùleitý, protoe slouí i jako vodítko pro vytvoøení náhradní matice tuhosti.

Kvalita rùzných

systémù MKP (pøesnost nahrazení díla modelem elastických prvkù) závisí pøedevím na tom, jaké tvarové funkce byly pouity.

Ze známého posunutí (2) lze v jednotlivých elementech stanovit deformace derivováním tvarových funkcí

du N d N i ε= = ∑ ui dx i=1 dx a pouitím konstitutivní rovnice (Hookova zákona) urèit napìtí

σ = E.ε .


11

Vypoètená napìtí, deformace a posuvy mohou být vyuity buï pøímo pro hodnocení funkènosti analyzované konstrukce (vèetnì stanovení její ivotnosti) nebo pro její optimalizaci (napø. optimalizaci tlouky skoøepin, prùøezù nosníkù, viz. dále popisovaný modul OPTSTAR).

Nelineární statika V nìkterých pøípadech by lineární analýza dala zcela falené informace o rozloení deformací a napìtí. Jde o situace, kdy tuhost konstrukce nebo pùsobitì sil se bìhem zatìování výraznì mìní: -

Velká posunutí nebo velké deformace (geometrická nelinearita); je to pøípad málo klenutých vík velkých nádrí nebo støech, dlouhých výloníkù ap.

-

Materiálové nelinearity (pøedevím plasticita, kterou musíme uvaovat napø. pøi dimenzování tlakových nádob nebo creep-teèení materiálu k nìmu dochází za vyích teplot).

-

Kontaktní úlohy-napø. pruný válec, zatlaèovaný do pruné podloky, pøípady, kdy se pøi zatìování vymezují vùle (napø. mezi loiskem a høídelem) ap.

Ve vech tìchto pøípadech analyzujeme konstrukci tak, e ji zatìujeme postupnì (pøírùstková metoda) a v kadém zatìovacím kroku poèítáme novou matici tuhosti. Ukáeme to na pøíkladu, kde významnou roli hrají velká posunutí, viz. následující obrázek 4. Poznámka: "Velká" posunutí je pojem relativní, jak z obrázku patrno; výrazná zmìna tuhosti se projeví jen tehdy, kdy úhel rozevøení obou táhel bude dostateènì velký, jinak by pøi zvyování zatíení

F

dolo nejprve ke

zplastizování nebo ke ztrátì stability táhel. V k-tém zatìovacím kroku bude vztah

∆F popisován rovnicí t t (k) t (k) r ( K L - K u - K σ )∆u = ∆F, resp. ( K L - K u - K σ )∆ u (k) = ∆ F (k) ,

mezi pøírùstkem posunutí

t KL

∆u

a zatíení

je matice tuhosti odpovídající výchozímu nezatíenému stavu,

vlivem zmìny geometrie (zvìtení úhlu rozevøení táhel);

t (k) Ku

t t (k) K L - Ku

kde

je matice, vyjadøující zmenení tuhosti

je matice tuhosti vypoètená "klasickým"

zpùsobem, ale pro "deformovanou" geometrii odpovídající k-tému zatìovacímu kroku (matice posunutí),

t (k) Kσ

t Ku

je funkcí

je tzv. matice geometrické tuhosti (té matice poèáteèních napìtí) a popisuje zmenení tuhosti

vlivem zmìny smìru pùsobení vnitøních sil. Jde o tento efekt: Na zaèátku k-tého zatìovacího kroku jsou vnitøní síly v táhlech

R(k)

v rovnováze s vnìjím zatíením

F (k) . Po zatíení pøírùstkem ∆ F (k)

se vak zvìtí úhel rozevøení


12

táhel a rovnováha reakcí s výchozím zatíením se poruí (a proto je matice geometrické tuhosti funkcí pøedpìtí vnitøních sil). Pøi vyetøování materiálových nelinearit (elastoplastického stavu) se postupuje podobnì, tj. po zatìovacích (nebo v pøípadì creepu èasových) krocích, a v kadém se poèítá efektivní matice tuhosti jako rozdíl klasické "lineární" matice tuhosti

t KL

a matice plasticity

t K p ( σ ′ ) , která závisí na okamitém stavu napjatosti.

Dynamika Aplikace dynamické analýzy jsou ménì èetné ne pøedchozí; toto jsou základní úlohy: -

Výpoèet vlastních frekvencí konstrukce. Zpravidla nás zajímá jen nìkolik nejmeních vlastních frekvencí,

kterým odpovídá nejmení útlum, nejvìtí amplitudy a tedy i nejvìtí dynamická napìtí. Obvykle známe budící frekvence (napø. nejnií frekvence, buzená kompresorem pøipojeným do potrubní sítì je dána frekvencí jeho otáèek nebo pro daný prùmìr kolony a rychlost vìtru lze spoèítat frekvenci periodického odtrhávání vírù ap.) a úkolem návrháøe je upravit konstrukci tak aby její vlastní frekvence byly co nejdále od frekvence budící a nedolo k rezonanci (u potrubních sítí toho lze dosáhnout napø. zmìnou ukotvení, tuhosti /kompenzátory/ èi pouitím tlumicích prvkù). -

Stanovení odezvy konstrukce na harmonickou budící sílu (napø. periodické budící síly vyvolané vìtrem pøi

obtékání válcových aparátù, dynamické úèinky rotujících prvkù ap.). Poadovaným výsledkem je obvykle závislost odezvy konstrukce (amplitud prùhybù, deformací a sil) na frekvenci budících sil a posouzení vlivu tlumení. -

Stanovení odezvy na zcela obecný prùbìh sil (napø. na skokovou nebo lineární zmìnu pùsobících sil).

Podobnì lze øeit pøípady, kdy je pøedepsán prùbìh zrychlení v urèitých místech. Analýza je zamìøena na popis vln deformací èi napìtí (napø. v souèástech zpracovatelských strojù). -

Náhodné buzení; pøedevím vliv seismicity. V tìchto pøípadech obvykle není znám konkrétní èasový prùbìh

zrychlení, ale tzv. výkonová hustota PSD (Power Spectral Density) budícího signálu (napø. prùbìh støední amplitudy zrychlení základù kolony v závislosti na frekvenci; zpravidla se uvauje konstantní PSD od nulové a po jistou mezní frekvenci). Výsledkem této analýzy jsou støední hodnoty posunutí a pøedevím napìtí. kolní pøípady volného i vynuceného kmitání (soustava s jedním stupnìm volnosti) jsou uvedeny na obrázku 5


Volné kmitání:

13

d2 u M 2 + Ku = 0 dt u = ua eiωt , (-M ω 2 + K)ua eiωt = 0, ω=

K . M

V obecném pøípadì konstrukce s N stupni volnosti popisuje volné kmitání rovnice

r t d2 u t r M 2 + Ku = 0, dt kde

t K

je matice tuhosti konstrukce a

t M

tzv. matice hmot (nìkdy se pouívá postup, pøi kterém je hmotnost celého

systému rozdìlena do uzlových bodù a rovnice (7) pak popisuje kmitání soustavy hmotných bodù propojených pruinami; matice hmot je v tomto pøípadì diagonální). Øeením (7) jsou harmonické prùbìhy posunutí

r r u k (t) = u k eiω t , k

kde

ωk

je k-tá vlastní frekvence a

r uk

vektor amplitud posunutí vech uzlových bodù

koneènìprvkového modelu (nazývá se k-tý vlastní tvar kmitù). Dosazením (8) do (7) se dostaneme k tzv. vlastní úloze,

t t r (- ω 2k M + K)u k = 0,

jejím øeením je právì N vlastních frekvencí a vlastních tvarù kmitù (modální analýza). Problém vynuceného kmitání (a event. tlumení) je pøirozenì sloitìjí, MKP ho

pøevádí na soustavu obyèejných diferenciálních rovnic, které jsou zobecnìním rovnice (7)

t d 2 ur t dur t r r M 2 + C + Ku = F(t), dt dt

kde kromì døíve uvedených matic tuhosti a hmot se objevuje matice tlumení

t C . I pro øeení této rovnice se pouívají

výsledky modální analýzy, které umoní pøevést (10) na soustavu N diferenciálních rovnic, které jsou nezávislé. Tento fakt jenom podtrhuje dùleitost metod øeení vlastní úlohy (9).


14

Stabilita Klasickým pøípadem je vzpìr nosníku; pokud je dostateènì tíhlý, lze pouít Eulerovu teorii, viz. následující obrázek Mez stability nosníku odpovídá situaci, kdy moment vyvolaný osovou silou

F.u

je

právì roven momentu, který odpovídá elastické deformaci nosníku, t.j. kde J-je moment setrvaènosti prùøezu nosníku. Rovnice (11) je úplnì stejná jako rovnice volného kmitání (6), místo poèáteèních podmínek (výchylky a rychlosti v èase nula) zde vak figurují podmínky okrajové u(0)=u(L)=0:

πx , u(x) = ua sin L

(-EJ

EJ

πx π2 =0 , 2 + F)ua sin L L

F krit =

EJ π 2 . L2

d2u + Fu = 0, dx 2

V obecném pøípadì je lineární stabilita konstrukce popisována rovnicemi, která jsou analogií rovnice (9) pro kmitání

t r t (-λ K σ + K)u a = 0,

kde

t K

je matice tuhosti poèítaná dle lineární teorie a

posunutí), spoèítaná pro zvolené zatíení konstrukce. tuhost

t t t λ K σ ; kdy je výsledná tuhost K - λ K σ

λ

t Kσ

je matice geometrické tuhosti (viz. analýza velkých

je násobek tohoto zatíení jemu odpovídá geometrická

nulová, dojde ke zhroucení konstrukce (tato pøedstava není úplnì

pøesná, protoe prvky matice geometrická tuhosti nejsou násobkem prvkù matice stability pøi takové hodnotì násobku referenèního zatíení

t K , ve skuteènosti dojde ke ztrátì

λ , kterému odpovídá nulový determinant výsledné matice

tuhosti, pøesnìji nejmení vlastní èíslo vlastní úlohy (13)). Take: kritické zatíení se poèítá úplnì stejnì jako nejmení vlastní frekvence kmitání, pouze místo matice hmot se uvauje matice geometrické tuhosti. Metodou vlastních èísel lze øeit stabilitu prutových konstrukcí, skoøepin, nádob zatíených vnìjím pøetlakem (pionýrské práce v této oblasti vykonal von Mises, kdy konstruoval trupy ponorek za první svìtové války). Problém je v tom, e vypoètená kritická zatíení nejsou na stranì bezpeènosti; ke skuteènému zhroucení konstrukce dojde mnohem døíve, zejména pro tyto pøípady (Køupka [12]):


15

- válcová skoøepina, zatíená osovým tlakem nebo ohybovým momentem. Dle lineární teorie je kritické napìtí závislé pouze na relativní tlouce skoøepiny (s/R)

σ krit = K 0 E

s , R

kde

K0

K0

je 2 a 10 krát mení. Tento rozpor lze vysvìtlit pøedevím pouitím teorie velkých deformací (pak vyjde

je pro tlaèenou válcovou skoøepinu teoreticky rovno 0.6, zatímco experimenty ukazují, e skuteèná hodnota

teoretická hodnota

K0

0.234) a dále uvaováním plastického ovlivnìní ztráty stability. Poznamenejme, e zhruba

stejná situace je u kulových skoøepin zatíených vnìjím pøetlakem; i tam je kritické napìtí dáno vztahem (14) s hodnotou

K 0 = 0.6

dle lineární teorie a 0.15 z výsledkù zkouek (kdy dochází jenom k lokální ztrátì stability).

- válcová skoøepina, zatíená vnìjím pøetlakem; pro nekoneènì dlouhou skoøepinu, vyplývá z lineární teorie stability vztah pro kritický pøetlak

p krit =

2E s 3 ( ) , 2 1- µ D

zatímco skuteèný kritický pøetlak je cca o 40% mení.

Teplotní pole Výpoèty teplotních polí patøí také do problematiky dimenzování aparátù, protoe teprve na základì vypoèteného rozloení teplot lze stanovit dilatace a obecnì teplotní napìtí. Teplotní pole je obvykle tøeba uvaovat jako nestacionární (a modelovat tøeba dìj nábìhu chemického reaktoru), protoe pøedem nevíme, zda budou dominovat dilataèní napìtí (odpovídající ustáleným provozním teplotám) nebo teplotní napìtí, zpùsobená velkými gradienty teplot pøi nábìhu aparátu. Jiné aplikace se týkají fázových pøemìn napø. v oblasti svaøování (viz. modul SYSWELD programu SYSTUS, který umí vyhodnocovat i materiálové zmìny /pøechody austenit-martenzit ap./) a øada dalích analýz, napø. stanovení úèinnosti izolací, výpoèet teplotních polí pøi tøískovém obrábìní,.... I kdy je to moná na první pohled zvlátní, tak moduly MKP pro výpoèet teplotních polí se pouívají i pro analýzu proudìní, pøièem se vychází z analogie: Napø. tok pórézním prostøedím popisuje stejná diferenciální rovnice jako stacionární rozloení teplot (místo teplot je tøeba uvaovat hydrostatickou výku a místo souèinitele tepelné vodivosti souèinitel propustnosti materiálu); takto se napø. modelují toky spodních vod (tøeba prosakování vody pod pøehradními tìlesy). Pøi výpoètech teplotních polí MKP se dalím uzlovým parametrem kromì posunutí promìnné teploty

T

a ty popisuje soustava obyèejných diferenciálních rovnic

r t dT t r r + K λ T = Q, Mc dt

u

stávají èasovì


kde

t Mc

16

je matice tepelných kapacit (stejnì jako matici hmot ji mùeme aproximovat diagonální maticí, její prvky

jsou hmotnosti uzlù krát mìrná tepelná kapacita

cp )

a

t Kλ

je matice tepelné vodivosti. Vektor pravé strany

r Q

zahrnuje zdroje tepla a okrajové podmínky: - objemové nebo bodové zdroje tepla (chemická reakce, indukèní nebo odporový ohøev) - konvektivní pøenos tepla povrchem (v tom pøípadì je tøeba u hranièních elementù zadávat souèinitel pøestupu tepla

α

a teplotu vnìjího prostøedí)

- radiaèní pøestup tepla (sáláním). Tento pøípad je nejchoulostivìjí, protoe radiaèní výkon není lineární funkcí teploty a celá úloha se rázem stane nelineární.

Dalí oblasti aplikací MKP Øada komerèních systémù MKP umoòuje i výpoèty proudìní (napø. program COSMOS øeení Navierových Stokesových rovnic pro laminární proudìní newtonských i mocninových kapalin v modulu FLOWSTAR, výpoèet laminárního i turbulentního proudìní v potrubních sítích s uvaováním charakteristik èerpadel ap. v modulu PIPESTAR ap.), výpoèty elekromagnetických polí (solenoidy ap., viz napø. modul ELSTAR programu COSMOS), akustické analýzy (modul SYSNOISE programu SYSTUS, umoòující stanovit napø. prùbìhy akustického tlaku v kabinì automobilu). Tìmito aplikacemi se dále zabývat nebudeme.

Jiné metody numerického øeení Prakticky vechny výe uvedené aplikace lze zvládnout i jinými numerickými metodami, jmenovitì - metodou sítí (koneèných diferencí, kontrolních objemù ...), která se uplatòuje pøedevím pøi výpoètech proudìní (zvlátì stlaèitelných tekutin, napø. v transonické oblasti) - metodou hranièních prvkù, která má hodnì spoleèného s MKP i kdy zatím zdaleka nedosáhla jejího rozíøení (a asi ani nedosáhne, protoe není tak univerzální a mùe být výhodnìjí pouze pøi analýze teplotních polí, nebo lineární elasticity). Její princip je zaloen na skuteènosti, e napø. teplotní pole nebo pole posunutí v homogením tìlese je plnì urèeno rozloením teplot a tepelných tokù, resp. posuvù a napìtí na povrchu tìlesa. To znaèí, e kdybychom znali posunutí i vnìjí zatíení (napø. tlak) na celém povrchu souèásti, mohli bychom pomìrnì snadno vypoèítat posunutí, deformace a napìtí v libovolném bodì uvnitø tìlesa (kadá hodnota by byla vyjádøena plonými integrály). Háèek je v tom, e my nikdy neznáme obì tyto velièiny souèasnì: napø. víme, e na urèité èásti povrchu souèásti je pøedepsaný tlak (a pak neznáme posunutí) zatímco na jiné èásti je zadané posunutí (a pak zase neznáme síly - reakce v uloení). V metodì hranièních prvkù se koneènými prvky a uzly pokrývá pouze povrch analyzované souèásti (napø. trojúhelníkovými elementy u trojrozmìrného objektu) a sestavuje se soustava lineárních algebraických rovnic jen pro neznámé posuvy/napìtí v tìchto hranièních elementech. Z toho je patrná základní výhoda: Pøi stejné pøesnosti nahrazení staèí mnohem mení poèet elementù a uzlových bodù, to znamená e staèí øeit soustavy rovnic s mením poètem neznámých (poznamenejme, e v MKP se i na perzonálních poèítaèích nebo pracovních stanicích bìnì øeí soustavy desítek tisíc lineárních algebraických rovnic) a jednoduí je i modelování objektu a vytváøení sítì. Tyto výhody jsou vak sporné: matice soustavy øeených rovnic je sice mení, ale není ji øídká a dokonce ani


17

symetrická, co znemoòuje pouít nìkteré speciální a vysoce efektivní metody øeení soustav jako je frontální metoda, LDL dekompozice nebo gradientní metody. Pokud v zadání vystupují objemové síly (nebo dokonce nelinearity) je stejnì nutné dìlit na koneèné prvky i vnitøek vyetøované oblasti. Nesporná je asi jen jedna výhoda: Dobøe se modelují dìje v (polo)nekoneèných oblastech, napø. teplotní pole v neomezeném homogenním prostøedí, které obklopuje tìleso ze zadanou teplotou (nebo pøípady vnìjího obtékání tìles), Brebbia [11].


2.

18

Principy MKP Základní principy MKP vysvìtlíme na pøíkladì návrhu programu pro statické øeení lineárních

dvourozmìrných prutových konstrukcí (soustavy kloubovì spojených pruinek). Ukáeme, e takový program lze sestavit i bez znalosti variaèních principù a sloité matematiky. Návrh algoritmu zaèíná vdy u jednoho obecného elementu, v naem pøípadì pruiny (táhla, tyèe) se zadanou tuhostí k. Pro obecnì orientované táhlo se dvìma uzlovými body 1,2, jsou ve høe ètyøi posunutí ( u x1 ,u y1 ,u x2 ,u y2 ) a ètyøi síly ( F x1 ,F y1 ,F x2 ,F y2 ). Pokusíme se odvodit matici tuhosti, která umoní vyjádøit lineární závislost mezi tìmito posuvy a vnitøními silami soustavou ètyø lineárních rovnic. Tato soustava by se dala napsat pøímo, ale ukáeme si postup, který je èasto výhodný ve sloitìjích pøípadech. Spoèívá v tom, e matice tuhosti se stanoví jednodueji v tzv. lokálním souøadném systému ( ζ ,η ), kde osa

ζ

je

totoná s osou pruinky; potom platí, e napø. síla pùsobicí v uzlu 1 ve smìru osy ( F ζ 1 ) je rovna souèinu tuhosti pruiny (k) a rozdílu posunutí v uzlech 1 a 2 ( uζ 1 - uζ 2 ). V maticovém tvaru:

 k  - k

- k   uζ 1  F ζ 1 ,  .  =   k   uζ 2  F ζ 2

E[Pa].S[ m2 ] kde tuhost k[N.m ] = L[m] -1

Sloky posunutí a síly v globálním kartézském souøadném systému (x,y) jsou se slokami v lokálním souøadném systému vázány lineární transformací (násobením transformaèní maticí

t T)

0  u x1  cosϕ        uζ 1 sin u ϕ 0   y1  .  ,   = 0 cosϕ   uζ 2  u x2       u y2  0 sin ϕ  a tatá transformace platí i pro sloky sil. Soustavu rovnic (17) lze nyní transformací (18) pøevést na poadovanou soustavu ètyø rovnic, které popisují vztah mezi posunutím a silami v globálním s.s. co po vynásobení matic na levé stranì dá koneèný tvar


0  cosϕ    sin ϕ 0  k   . 0 cosϕ   - k    0 sin ϕ    cos2 ϕ   cosϕ . sin ϕ k. - cos2 ϕ    - cosϕ . sin ϕ

t K

19

0  F x1  u x1  cosϕ       0 0  u y1  sin ϕ 0  F ζ 1  F y1 - k   cosϕ sin ϕ  =  ,  .  =   .  .  k  0 0 cosϕ sin ϕ   u x2  0 cosϕ   F ζ 2  F x2        F y2  u y2  0 sin ϕ 

cosϕ . sin ϕ

- cos2 ϕ

sin2 ϕ

- cosϕ . sin ϕ

- cosϕ . sin ϕ

cos2 ϕ

- sin2 ϕ

cosϕ . sin ϕ

- cosϕ . sin ϕ  u   x1  F x1 - sin2 ϕ   u y1  F y1  .  =   , cosϕ . sin ϕ   u x2  F x2      sin2 ϕ   u y2  F y2

t r r - - - - - - - - - - - - - - - - - K - - - - - - - - - - - - - - - - - - -u - - -F je matice tuhosti jednoho obecného tyèového prvku; my ovem potøebujeme matici celé soustavy pro libovolný

poèet libovolnì orientovaných prvkù. Pokud má soustava M uzlových bodù, bude rozmìr této globální matice tuhosti N=2.M (protoe kadému uzlu odpovídají 2 stupnì volnosti):

 F x1  K 11 K 12 . . K 1 2M   u x1       F y1 .  u y1  K 21 K 22 . .      . . . . . . . =  . .       . . . . .  . .         K 2M 1 . . . K 2M 2M   u yM   F yM  Prvky globální matice tuhosti získáme sèítáním odpovídajících prvkù matic tuhosti vech tyèových elementù (globální matice tuhosti se nejprve vynuluje a pak probíhá cyklus pøes vechny elementy: v kadém kroku se pøièítá 16 prvkù lokální matice tuhosti (20) k pøísluným prvkùm globální matice (21) /viz. následující program/). Korespondence mezi øádky a sloupci lokální a globální matice tuhosti je dána tzv. maticí konektivity, její øádky obsahují èísla uzlù (1,2,...,M) jednotlivých elementù (matice konektivity má v naem pøípadì 2 sloupce /protoe kadý prvek má právì dva uzlové body/ a poèet øádkù je roven poètu elementù). Vektor pravé strany (zatíení) obsahuje souèty sil pruin, které pùsobí v jednotlivých uzlech a má-li být zachována rovnováha musí být prvky tohoto vektoru pøímo vnìjí síly (v kloubech /-uzlech/, kde ádná vnìjí síla nepùsobí, se zadává 0). Soustava rovnic (21) jetì stále není v této fázi øeitelná, protoe nezahrnuje okrajové podmínky (to se projeví tím, e

u0i , pøedepsaná ve vybraných ui = u0i . Nejjednoduí úprava

matice tuhosti je singulární1). Okrajové podmínky musí být v naem pøípadì alespoò 3 posunutí uzlech. Kadá tato podmínka znamená, e je tøeba nahradit jednu rovnici soustavy (21), rovnicí

1

Chybové hláení SINGULAR STIFFNESS MATRIX, se kterým se budete asi nejèastìji setkávat, je zpravidla zpùsobeno právì opomenutím okrajových podmínek (nebo je jich

málo a konstrukce je v nìkterém smìru volnì pohyblivá).


20

soustavy (21), zajiující alespoò pøibliné splnìní tìchto podmínek je zaloena na následujícím triku: diagonální prvek matice tuhosti, který odpovídá pøedepsanému posunutí

u0i

se vynásobí velmi velkým èíslem (napø.

α = 1010 ) a do pravé strany se

místo zatíení dosadí posunutí, vynásobené upraveným diagonálním prvkem; i-tá rovnice bude mít po této úpravì tvar

α K ii ui (+ K i1 u1 + K i2 u2 +... ) = α K ii u0i , a je snad patrné, e pøi zanedbání relativnì malých èlenù uvedených v závorce, bude splnìní okrajové podmínky zajitìno. Výhoda tohoto postupu je v tom, e staèilo zmìnit jediný prvek matice tuhosti a navíc jeho pùvodní hodnota není ztracena; matici

t K

lze pozdìji rekonstruovat (dìlením pøísluného diagonálního prvku konstantou

α)

a pouít pøi výpoètu reakcí (sil, která pùsobí v tìch uzlech, kde bylo pøedepsáno posunutí). Po tìchto úpravách se ji soustava rovnic dá øeit (napø. Gaussovou eliminací, i kdy v MKP se pouívají mnohem efektivnìjí postupy) a výsledkem je vektor posunutí uzlových bodù. Síly (a napìtí) pùsobící v jednotlivých prutech a reakce v uloení se zjistí a následovnì z (20), resp. (21) (tj. násobením matic tuhosti vektorem vypoètených posunutí). Text programu, napsaného v jazyce Fortran (tento jazyk se v MKP pouívá prakticky vdy), je následující:

C C Øeení dvourozmìrných prutových konstrukcí: lineární elasticita, malé prùhyby i deformace C Vstupní data: soubor "TRUSS.DAT" C DIMENSION A(10000) OPEN(10,file='TRUSS.DAT') READ(10,*)M,NE,NDF,NDU C M-poèet uzlù, NE-poèet prutù, NDF-poèet zadávaných sil, NDU-poèet zadávaných posuvù. C C Dynamické rozdìlení pamìti; slouí k tomu aby operaèní pamì poèítaèe byla maximálnì vyuita. Program bude trochu sloitìjí, ale je to profesionální postup. C Vektor a(10000) /operaèní pamì, kterou pøidìlíme datùm/ rozdìlíme na úseky, odpovídající potøebným pracovním vektorùm, jejich velikost je dána C poadavky právì øeené úlohy (tj. hodnotami M,NE,NDF,NDU). Význam tìchto úsekù je tento: C C X(M),Y(M),U(2.M),F(2.M),K(2.M,2.M),IUE(NE,2),S(NE),DF(NDF),KDF(NDF),DU(NDU),KDU(NDU) C X,Y - vektory souøadnic uzlových bodù C U,F - vektor posunutí uzlù a vektor pravé strany globální soustavy (vektor vnìjích sil) C K - matice tuhosti globální soustavy C IUE - matice konektivity (dvojice èísel uzlù jednotlivých elementù) C S - plochy prùøezù jednotlivých prutù C DF - pøedepsané sloky vnìjích sil; ve kterých uzlech urèuje následující vektor KDF C KDF - indexy uzlù, v nich se zadává pøísluná sloka síly FD; kladná hodnota KDF znaèí, e DF oznaèuje sloku ve smìru X, zatímco záporná hodnota sloku Y. C DU - pøedepsané sloky posunutí; které sloky a ve kterých uzlech urèuje následující vektor KDU C KDU - indexy uzlù v nich se zadává posunutí DU; platí stejná konvence jako u KDF. C C Dle tohoto návrhu rozmístìní zaèíná napø. matice konektivity IUE(NE,2) na pozici C LIUE=6*M+4*M*M+1 vektoru A(10000). C LIUE=6*M+4*M*M+1 LKDF=LIUE+3*NE LKDU=LKDF+2*NDF C Ètení dat a vlastní výpoèty pøenecháme proceduøe COMP(M,NE,NDF,NDU,X,Y,U,F,K,IUE,S,DF,KDF,DU,KDU) CALL COMP / (M,NE,NDF,NDU,A(1),A(M+1),A(2*M+1),A(4*M+1),A(6*M+1),A(LIUE),A(LIUE+2*NE),A(LKDF),A(LKDF+NDF),A(LKDU),A(LKDU+NDU) ) END SUBROUTINE COMP(M,NE,NDF,NDU,X,Y,U,F,K,IUE,S,DF,KDF,DU,KDU) REAL X(M),Y(M),U(2*M),F(2*M),K(2*M,2*M),IUE(NE,2),S(NE),DF(NDF),KDF(NDF),DU(NDU),KDU(NDU),KL(4,4),L INTEGER IP(4) C Ètení dat: modul prunosti, souøadnice uzlù, matice konektivity, prùøezy vech prutù, zatíení, pøedepsaná posunutí READ(10,*)E,X,Y,((IUE(I,J),J=1,2),I=1,NE),S,DF,KDF,DU,KDU C Nulování matice tuhosti i vektoru pravé strany K=0. F=0.


21

C Sestavení globální matice tuhosti v cyklu pøes NE-prutù DO IE=1,NE C I1,I2-indexy uzlù, IP(4)-indexy uzlových parametrù (posunutí uzlù) elementu IE I1=IUE(IE,1) I2=IUE(IE,2) IP(1)=2*I1-1 IP(2)=2*I1 IP(3)=2*I2-1 IP(4)=2*I2 C Výpoèet lokální matice tuhosti prutu èíslo IE dle vztahu (20) (L-délka prutu, CS-cosinus fi,...) L=SQRT((X(I1)-X(I2))**2+(Y(I1)-Y(I2))**2) SS=(Y(I2)-Y(I1))/L CS=(X(I2)-X(I1))/L RK=E*S(IE)/L KL(1,1)=RK*CS**2 KL(1,2)=RK*CS*SS .... KL(4,4)=RK*SS**2 C Pøiètení lokální matice tuhosti elementu IE ke globální matici DO I=1,4 DO J=1,4 K(IP(I),IP(J))=K(IP(I),IP(J))+KL(I,J) ENDDO ENDDO ENDDO C Vnìjí zatíení DO JF=1,NDF IF(KDF(JF).GT.0)THEN F(2*KDF(JF)-1)=DF(JF) ELSE F(-2*KDF(JF))=DF(JF) ENDIF ENDDO C Zadávaná posunutí DO JU=1,NDU IF(KDU(JU).GT.0)THEN I=2*KDU(JU)-1 ELSE I=-2*KDU(JU) ENDIF K(I,I)=K(I,I)*1E10 F(I)=K(I,I)*DU(JU) ENDDO C Øeení soustavy 2.M lineárních algebraických rovnic Gaussovou eliminací (GELG je procedura fy. IBM, která výsledek /POSUNUTÍ/ vrací do vektoru F) CALL GELG(F,K,2*M,1,1E-6,IERROR) WRITE(*,1)F END

Jako ukázku vstupních dat uvedeme pøíklad dvou symetrických táhel se zatíeným styèníkem

E = 2.1011 Pa, S = 10 -4 m2 , F = 1000 N Data pøipravená v souboru TRUSS.DAT: 3 (uzly) 2 (elementy) 1 (síla) 4 (posuvy) 0 1 0.5 (souøadnice x) 001 (souøadnice y) 1 3 2 3 (matice konektivity) 1E-4 1E-4 (prùøezy prutù) -1000 -3 ( F y v uzlu 3)

0 0 0 0

( u x1 ,u y1 ,u x2 ,u y2 )

1 -1 2 -2 Pozn.: Pro praktickou aplikovatelnost tomuto programu nìco chybí; výpoèet napìtí v prutech a reakcí (zatím je výsledkem pouze posunutí uzlù zatíené konstrukce). Dokáete to?


22

Vyuívání komerèních programù analýzy MKP vysvìtlíme na systému COSMOS/M fy. S.R.A.C. St.Monica, který je v ÈR pomìrnì rozíøený (je integrován m.j. do nejnovìjích verzí programù AutoCAD nebo Pro/Engineer). I u jiných systémù je postup analogický, lií se jen názvy pøíkazù, jejich parametry a pøedevím bohatostí knihoven elementù a algoritmù øeení.

Geometrický model konstrukce

Základní entity, z nich je tvoøen model jsou POINT (vztaný bod), CURVE (køivka), SURFACE (povrch, troj nebo ètyøúhelníková plocha), REGION (oblast vymezená køivkami, vícenásobnì souvislá plocha), VOLUME (objem, ètyø nebo estistìn), PH (PolyHedra - uzavøený mnohostìn), PART (souèást, vymezená mnohostìny /vícenásobnì souvislá oblast s dutinami/). Pøedstavte si tøeba tenkostìnnou nádobu (skládanou z entit typu SURFACE) a pøipojené potrubí (modelované køivkami) nebo geometrický model pøíruby (sloený z prstencùobjemù). Postup modelování geometrie je prakticky stejný jako napø. u programu Autocad (i kdy to není obvyklé lze geometrický model vytvoøit pøímo Autocadem a pøenést ho do systému COSMOS jako soubor typu DXF). Definice vztaných bodù (POINTS, pozor, to jetì nejsou uzlové body): PT definice zadáním souøadnic x,y,z PTINTCC,PTINTCS bod jako prùseèík (INTersection) dvou køivek (Curves) nebo køivky a plochy (Surface) Definice prostorových køivek (CURVES): CRLINE, CRPLINE, CRSPOLY spojení dvou bodù úseèkou nebo definice polygonu (Polylines) CRARC, CRPCIRCLE, krunice a kruhové oblouky CR4PT, CRHELIX, CRELLIPSE, kubická køivka procházející 4 body, roubovice, elipsa CRSPLINE, CRBEZCORD, spliny a Bezierovy køivky CRFILLET,CRJOIN,CRBLENDS, spojování dvou køivek (FILLET-zaoblení s daným polomìrem ap.) CRINTSS, køivka definovaná jako prùseènice dvou ploch CRTANPT,CRTANLIN,CRNORMPT, teèna z bodu ke køivce, spoleèná teèna køivek, normála ke køivce Definice prostorových ploch-zakøivené troj nebo ètyøúhelníky (SURFACES): SF3PT,SF4PT, trojúhelníková, resp. ètyøúhelníková ploka definovaná vztanými body (vrcholy) SF2CR,SF4CR,SF3CR,SF4PCR, ètyøúhelníková ploka definovaná dvìma nebo ètyømi køivkami, trojúhelníková plocha, plocha, procházející ètyømi paralelními køivkami SFEXTR,SFSWEEP,SFGLIDE,SFDRAG, vytváøení ploch pohybem tvoøících køivek (EXTRusion-taení ve smìru osy souø. systému, SWEEP-rotace tvoøící køivky, GLIDE-tvoøící køivka "kloue" po vodící køivce, pøièem si zachovává smìr /translaèní pohyb/, DRAG-toté co GLIDE, ale tvoøící køivka zachovává úhel vzhledem k vodící køivce /a ta musí být hladká aby v místech skokové zmìny smìrnice nevznikly ve vytváøené ploe "díry"/). Definice objemù-ètyø nebo estistìnù (VOLUMES): VL8PT,VL4CR,VL2SF,VLCRSF, estistìny definované 8-mi vrcholy, 4-mi paralelními køivkami, dvìma protilehlými povrchy, plochou a køivkou VLEXTR,VLSWEEP,VLGLIDE,VLDRAG, vytváøení objemù pohybem tvoøících ploch (taením, rotací, klouzáním). Definice souèástí (PARTS):


23

PH,PHEXTR,PHSWEEP, povrch mnohostìnu, definovaný výètem ploch nebo oblastí (REGIONS) PART, objemová entita, definovaná výètem mnohostìnù PH, které tvoøí vnìjí i vnitøní povrchy (dutiny). Pøedchozí pøíkazy umoòují definovat vdy jen jednu entitu; existují dalí, které s ní umí manipulovat (posunout, natoèit) èi duplikovat (napø. pøi modelování otvorù pro rouby v pøírubì vytvoøíme pøíkazem CRCIRCLE krunici /obvod otvoru/, duplikujeme ji pøíkazem CRGEN /parametrem je úhel natoèení otvorù a jejich poèet/ a vlastní plochy otvorù vytvoøíme taením vzniklé mnoiny køivek ve smìru osy SFEXTR).

Volba typu elementu a generování sítì

Kadý systém MKP obsahuje knihovnu rùzných typù koneèných prvkù, které lze vìtinou vzájemnì kombinovat; napø. u tlakové nádoby pouijeme pro modelování plátì skoøepinové prvky (pokrývající plochy /SURFACES/ geometrického modelu) a jednorozmìrné prvky typu trubka /PIPE/ (pokrývající køivky modelující pøipojená potrubí). Nìkteré elementy programu COSMOS jsou uvedeny v následující tabulce: MASS

hmotný bod, 1 uzel, 6 stupòù volnosti (DOF /Degree Of Freedom/, 3

TRUSS3D

tyèový prvek, 2 uzly, 6 DOF (posuvy)

BEAM3D

nosníkový prvek, 2 uzly, 12 DOF (posuvy a natoèení)

PIPE

trubka, 2 uzly, 12 DOF (podobnì jako BEAM3D, ale uvauje se vnitøní

SHELLAX

rotaènì sym. skoøepina, 2 uzly, 6 DOF (tenký kuelový prstenec)

PLANE2D

ètyøúhelníkový zakøivený prvek, 4 nebo 8 uzlù (DOF jsou pouze posuvy), umoòuje modelovat rovinnou deformaci, rovinnou napjatost, prstenec se ètyøúhelníkovým prùøezem, nevazkou kapalinu. Pouitelný pro velké deformace i posuvy, plasticitu

TRIANG

zakøivený trojúhelníkový prvek, 3 nebo 6 uzlù, podobné vlastnosti jako

SHELL3

trojúhelníkový nebo ètyøúhelníkový skoøepinový prvek (dle teorie tenkých

SHELL4

skoøepin) s tøemi a devíti uzlovými body (DOF jsou posuvy i natoèení)

SHELL9 SHELL4T

tlustostìnná skoøepina

SHELL4L

laminátová skoøepina (sloená z vrstev rùzných vlastností)

SOLID

estistìn s 8-mi a 20-ti uzlovými body (DOF jsou pouze typu posunutí)

GAP

speciální dvouuzlový prvek-vymezování vùlí, vliv tøení

RLINK

dvouuzlový prvek reprezentující pøenos tepla sáláním

Elementy, které budou pouity k modelování zvolené èásti modelu se volí pøíkazem EGROUP.


24

Generování sítì koneèných elementù (tj. i vytváøení uzlových bodù) se provádí dvìma zpùsoby: - parametrické generování sítì (vpodstatì pravidelné dìlení køivek, ploch nebo objemù na zvolený poèet elementù /s moností zhuování sítì v urèitém smìru/); pøíkazy M_CR, M_SF, M_VL (Mesh CuRve, SurFace, VoLume). - automatické generování sítì (v tom pøípadì musí pøedcházet dalí krok: definice oblastí /REGIONS nebo PARTS/ vymezených hranicemi /CONTOURS nebo PH/, které se skládají z døíve ji definovaných køivek èi ploch; pøíkazy CR, RG, PH, PART). Takto vymezené oblasti jsou potom automaticky pokrývány nepravidelnì uspoøádanými elementy, pøièem se zadává buï jejich poèet nebo prùmìrná velikost; pøíkazy MA_RG, MA_PH (Mesh Automatic ReGion nebo PolyHedra pro ploné elementy SHELL) nebo MA_PART (pro ètyøstìnné elementy TETRA). Pozn.: Nejnovìjí programy MKP umoòují automaticky pøedefinovávat sí elementù dle výsledkù pøedchozí analýzy (napø. zhustit sí v oblasti vysokých gradientù napìtí nebo posunutí). Je to tak zvaná H-metoda, která se dokonce nìkdy kombinuje s P-metodou, která volí automaticky i typ elementù (pøesnìji stupeò polynomù tvarových funkcí prvkù); viz. pøíkaz ADAPTIVE.

Zadání materiálových parametrù elementù Pro statické výpoèty uvaující lineární izotropní materiál, staèí zadat modul prunosti E /oznaèení EX/ a Poissonovu konstantu

µ

/NUXY/ nebo modul prunosti ve smyku G /GXY/. Uvaujeme-li vliv teplotního pole je

tøeba specifikovat i souèinitel teplotní roztanosti vodivosti

λ

α

/ALPX/ a pokud teplotní pole poèítáme i souèinitel tepelné

/KX/ a mìrnou tepelnou kapacitu materiálu

cp

/C/. Pøi analýze kmitání nebo pøi uvaování vlivu

zrychlení (tøeba zatíení vlastní vahou) musí být známa i hustota

ρ

/DENS/. V nìkterých pøípadech nelze

pøedpokládat, e materiál je izotropní (døevo, kompozity) a v tom pøípadì se elastické konstanty i souèinitele tepelné vodivosti a teplotní roztanosti specifikují pro jednotlivé smìry vztaené k lokálnímu souøadnému systému elementu ( Ex ,

E y ,... λ x ,... α x ,... ).

U nelineárních modelù je tøeba pracovat s funkèní závislostí materiálových

parametrù; u nejjednoduích modelù elastoplacity staèí napø. mez kluzu /SIGYLD/ a teèný modul prunosti /ETAN/, pøi uvaování creepu konstanty pøísluného modelu (napø. 7 konstant odpovídajících exponenciálnímu modelu CREPX). Obvyklý je popis funkèní závislosti tabulkou hodnot (v ní pak program interpoluje); takto se zadávají napø. závislosti modulù prunosti, souèinitelù tepelné vodivosti ap. na teplotì (pøíkaz CURDEF), nebo køivky závislosti napìtí - deformace (pøíkazem MPC).

Základní pøíkazy: MPROP - zadání libovolného materiálového parametru (EX,NUXY,...) pro urèitou skupinu elementù RCONST - zadávání dalích parametrù jako napø. tloutìk skoøepin (zde tøeba nahlédnout do pøedpisu pro jednotlivé typy elementù), PIC_MAT-výbìr materiálových parametrù pro zvolený materiál z pøipravené databáze, která zahrnuje napø. tyto


25

materiály (PC_STEEL-uhlíková ocel, CS_STEEL-nerez CF-8M nebo CF-20, MC_IRON-tvárná litina ASTM-A220, D_NICKEL-Duranickel 301, W_COPPER-slitiny mìdi) :

Oznaèení

EX

GXY

NUXY

ALPX

DENS

E [MPa]

G [MPa]

µ [-]

α[K ]

ρ [kg.m ] λ [W.m . c p [J. kg

PC_STEEL

206840

79290

0.28

1,332.10 -5

7834

43

430

CS_STEEL

193050

79290

0.26

1,512.10 -5

7750

16.3

523

MC_IRON

186200

86180

0.27

1,206.10 -5

7336

47

519

D_NICKEL

206800

77220

0.34

1,386.10 -5

8416

14.8

419

W_COPPER

106870

39990

0.33

1,980.10 -5

8442

26

377

-1

KX

C -1

-3

-1

Poznamenejme, e COSMOS poèítá implicitnì v anglosaských jednotkách (které poaduje napø. americká norma pro navrhování tlakových nádob ASME), kde jednotkou tlaku je psi [=6.8948 kPa], délky inch [palec=2.54 cm], hmoty lbm [libra=0.4536 kg], síly lbf [váhy=4.448 N] a teploty F [Fahrenheit]. V tìchto jednotkách by napø. hodnoty materiálových parametrù pro bìnou uhlíkovou ocel byly:

Oznaèení

PC_STEEL

EX

GXY

NUXY

ALPX

DENS

KX

C

E [Mpsi] G [Mpsi] µ [-]

α[F ]

ρ [lb.s .in λ [BTU.in c p [BTU.in

30

7,4.10 -6

0.0007324

11.5

0.28

-1

2

0.000575

40.572

Zatíení a okrajové podmínky (Forces&Displacements) V této fázi formulace úlohy jsou ji vytvoøeny dvì databáze: geometrický model (seznamy vztaných bodù, køivek, ploch, objemù) a prvkový model (seznam uzlù, elementù a materiálových parametrù). Zadávané síly,


26

posunutí, teploty, atd. se vztahují k prvkové databázi (lze napø. zadat sílu pùsobící v uzlu /NODE/, ale ne ve vztaném bodì /POINT/). Systém COSMOS nicménì geometrický model pouívá i v této etapì: pøíkazem P_SF lze napø. zadat zatíení tlakem, který pùsobí na plochu (co je geometrická entita) a program sám vyhledá elementy (tøeba SHELL), které na této ploe leí a kadému z nich pøedepsané zatíení pøiøadí. Posunutí-DISPLACEMENTS DND-pøedepsané posunutí nebo natoèení v uzlu (a v zadaném smìru, uívá se oznaèení UX,UY,UZ,RX,RY,RZ) DCR,DSF,DCT-pøedepsaná posunutí nebo natoèení ve vech uzlech leících na køivce, povrchu nebo hranici oblasti (CuRve, SurFace, ConTour). Osamìlé síly-FORCES FND,FPT,FCR,FSF,FCT-pøedepsané síly (FX,FY,FZ) nebo momenty (MX,MY,MZ) v uzlu nebo ve vech uzlech leících na vyjmenovaných køivkách, povrích èi hranicích. Zatíení tlakem-PRESSURE PEL,PCR,PSF,PRG-pøedepsaný tlak na jednu stranu elementu ev. elementù napojených na køivku, povrch nebo oblast (ReGion). Zrychlení-GRAVITY ACEL,OMEGA,DOMEGA-zrychlení, úhlová rychlost, úhlové zrychlení. Teploty-TEMPERATURES NTND,NTCR,NTSF,NTVL-pøedepsaná teplota v uzlu nebo uzlech napojených na køivku, plochu, objem. Bodové zdroje tepla-NODAL HEAT QND,QCR,QSF,QVL,QRG-tepelný výkon [W] generovaný v uzlu nebo uzlech napojených na uvedené entity Objemový zdroj tepla-ELEMENT HEAT -3

QEL,QECR,QESF,QEVL,QERG-objemový zdroj tepla [ W.m ] v elementu nebo elementech napojených na uvedené entity Tepelný tok-HEAT FLUX HXEL,HXCR,HXSF,HXRG-pøedepsaná hustota tepelného toku stìnou elementu (elementù spojených s entitami) Souèinitel pøestupu tepla-CONVECTION CEL,CECR,CESF,CERG-souèinitel pøestupu tepla a teplota vnìjího prostøedí. Stejnì jako bylo moné zadávat napø. teplotní závislosti materiálových parametrù formou tabulek, lze pøedepsat i èasové prùbìhy sil a okrajových podmínek; pøíklad CURDEF,TIME,1, - definice èasového prùbìhu èíslo 1 - zadáváme dvojice: èas, funkèní hodnota ACTSET,TC,1 PSF,3,1500;

- aktivujeme èasový prùbìh èíslo 1

- zatíení tlakem na ploe 3; tlak ale není 1500 Pa (jak by tomu bylo pøi vynechání pøíkazu ACTSET),


27

nýbr 1500 x hodnota, interpolovaná z èasového prùbìhu èíslo 1 v aktuálním èase (týká se samozøejmì pouze nestacionární analýzy).

Provedení analýzy Odstartování vlastního výpoètu; teprve v této fázi se budou vytváøet matice tuhosti, hmot, vodivostí ap. a øeit soustavy algebraických nebo diferenciálních rovnic pro uzlové parametry (posuvy, teploty) nebo vlastní èísla (frekvence kmitání èi kritická zatíení). Následující zpracování zahrnuje výpoèet deformací, napìtí a event. ivotnosti konstrukce. Základní moduly a volitelné metody øeení (jejich pøesnìjí popis bude uveden pozdìji nebo ho najdete v doporuèené literatuøe), které pouívá program COSMOS jsou uvedeny v následující tabulce

STATIC (lineární elasticita, malé deformace) NONLINEAR

t r r K.u = F t r r K∆u = ∆F

(elastoplasticita, velké deformace, creep) BUCKLING (lineární stabilita) FREQUENCY (vlastní

Eliminaèní metoda zaloená na blokovém rozkladu matic (sky line) Newtonova Raphsonova metoda a její modifikace

t t r (K - λ K σ )u = 0 t t r (K - ω 2 M)u = 0

kmity)

Výpoèet nejmení vlastní hodnoty (násobku zatíení) inverzní iterací Inverzní a Jacobiho iterace, Iterace podprostoru, Lanczosova metoda, Guayanova redukce

DYNAMIC (vynucené kmitání, seismicita)

r r t d 2 u t du t r r M 2 + C + Ku = F(t) dt dt

FATIQUE (únava materiálu)

modální superpozice, Duhamelova integrace, Newmarkova metoda

èirá empirie odpovídající metodice amerických pøedpisù ASME

OPTIMIZATION

citlivostní analýza (Taylorùv

(optimalizace geometrie

rozvoj) a metoda feasible

návrhu) THERMAL (nestacionární teplotní pole)

Pøíkazy:

r r t dT t r + Kλ T = Q Mc dt

directions pøímá integrace (implicitní metoda)


28

A_STATIC, A_STRESS, A_FREQUENCY, A_BUCKLING, A_NONLINEAR, A_FATIQUE, A_THERMAL specifikace podmínek øeení (volba metody, poèet iterací, poadovaná pøesnost, atd.) R_STATIC, R_FREQUENCY, R_BUCKLING, R_NONLINEAR, R_FATIQUE, R_THERMAL, R_DYNAMIC start analýzy.

Forma výsledkù Tiskové sestavy-tabulky posunutí a natoèení uzlù, síly a napìtí pùsobící v elementech. Tyto údaje jsou uloeny v souboru (problem.OUT), který lze vytisknout nebo pouít jako soubor vstupních dat pro vá vlastní program, který tato data vyhodnotí. Izoèáry-vrstevnice napìtí, posuvù, teplot, doby ivotnosti Grafy-èasové prùbìhy teplot, napìtí, posuvù ve vybraných uzlech Animace-pohyblivé znázornìní deformací konstrukce, tato "hraèka" slouí pouze k rychlému odhalení chyb v zadání okrajových podmínek a zatíení. Jako ilustraci zpùsobu zadávání problému napime program pro øeení pøedchozí úlohy (dvì táhla zatíená silou, viz. pøedchozí pøíklad) PT,1,0,0,0, definice vztaných bodù 1(0,0), 2(1,0), 3(0.5,0) PT,2,1,0,0, PT,3,.5,1,0, CRPLINE,1,1,3,2,2, polygon (vytvoøení køivek 1 a 2 /lines/) EGROUP,1,TRUSS2D,0,0,0,0,0,0,0, volba typu elementu (skupina 1) MPROP,1,EX,2E11,

materiálové parametry (

E,µ )

MPROP,1,NUXY,.3, RCONST,1,1,1,2,1E-4,0, prùøez prutù M_CR,1,2,1,2,1,1, generátor sítì dvouuzlových prvkù na køivkách 1 a 2 NMERGE,1,4,1,0.0001,0,1,0, slouèení uzlù se stejnými souøadnicemi DND,1,UZ,0,4,1,, nulové posuvy ve smìru z DPT,1,UX,0,2,1,UY,,, pøedepsané posuvy ve vztaných bodech 1 a 2 FPT,3,FY,-1000,3,1, zadaná síla v bodì 3 R_STATIC lineární analýza

Výsledky jsou uloeny v souboru s pøíponou OUT, uvedeme jen malý výsek: DISPLACEMENTS NODE X-DISPL. Y-DISPL. Z-DISPL. 1 0.00000E+00 0.00000E+00 0.00000E+00 2 0.00000E+00 -3.49386E-05 0.00000E+00 4 0.00000E+00 0.00000E+00 0.00000E+00

XX-ROT. 0.00000E+00 0.00000E+00 0.00000E+00

S T R E S S E V A L U A T I O N FOR S T A T I C ELEMENT FORCE STRESS NUMBER 1 -0.559017E+03 -0.559017E+07 2 -0.559017E+03 -0.559017E+07

YY-ROT. ZZ-ROT. 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00

ANALYSIS


29

Napìtí a deformace Základním výsledkem vech analýz budou tenzory napìtí (STRESS), které lze vzhledem k uèitému souøadnému systému (zpravidla lokálnímu souøadnému systému elementu) vyjádøit 9-ti hodnotami

 σ xx σ xy σ xz   t  σ = σ yx σ yy σ yz      σ zx σ zy σ zz σ xy oznaèuje napìtí ve smìru y, pùsobící na rovinu, která jest kolmou k ose x. Tenzor napìtí je symetrický ( σ xy = σ yx vyplývá z rovnice momentové rovnováhy infinitezimálního objemu) take pouze 6 jeho

kde napø.

sloek je nezávislých. V bìné praxi se tyto sloky chápou jako sloupcový vektor (i kdy to u není vektor v tenzorovém smyslu, platí pro nìj pravidla maticového a ne tenzorového poètu): T r σ = σ xx σ yy σ zz σ xy σ yz σ zx .

(

)

Nesporná výhoda tenzorového zápisu je to, e lze snadno zjistit napìtí, které pùsobí v libovolnì orientované rovinì (øezu, urèeném vektorem normály

n x ,n y ,nz ); pro sloku napìtí ve smìru i toti platí

fi =

3

∑σ

ij

j=1

nj

r t r (f = σ .n). r

Hlavní napìtí. Pro zvlátní orientaci øezu (tj. pro urèitý vektor n ) bude výsledné napìtí kolmé k rovinì øezu

r

(rovnobìné s normálou n ):

σ ni =

3

∑σ j=1

ij

t r r n j ( σn = σ .n).

r σ nazýváme hlavním napìtím a smìr n smìrem hlavního napìtí. Jak tento smìr zjistit? Kdybychom σ znali, byla by (26) vlastnì soustava tøí lineárních algebraických rovnic pro tøi neznámé souøadnice ni . Jenome tato soustava je homogenní, take nenulové øeení mùe mít pouze tehdy, pokud je její

Toto napìtí

matice singulární, tj. kdy

t t det( σ - σδ ) = 0.

A to je poadavek, který lze splnit jen pro urèité hodnoty

σ ; po rozepsání determinantu (27) se ukáe, e tato

podmínka je kubickou rovnicí, která má vdy tøi reálné koøeny a pro libovolný stav napjatosti tedy existují tøi hlavní napìtí a jim odpovídající tøi hlavní smìry (bylo by moné dokázat, e jsou vzájemnì kolmé). Zmínìná kubická rovnice pro hlavní napìtí má tvar

σ 3 - I σ σ 2 + II σ σ - III σ = 0,


30

kde I,II,III jsou tzv. invarianty tenzoru napìtí (viz. napø. esták [22])

I σ = σ xx + σ yy + σ zz = σ 1 + σ 2 + σ 3 = σ ii , II σ = σ xx σ yy + σ zz σ xx + σ yy σ zz - σ 2xy - σ 2xz - σ 2yz = σ 1 σ 2 + σ 2 σ 3 + σ 3 σ 1 =

1 2 [( σ ii ) - σ ij σ ji ], 2

t III σ = det( σ ) = σ 1 σ 2 σ 3 = ε ijk σ 1i σ 2j σ 3k . Pozn.: V pøípadì, e nìkteré hlavní napìtí je nulové (napø.

σ 3 = σ zz = 0 ), bude nulový i tøetí invariant (29) a

(28) se redukuje na kvadratickou rovnici s øeením

σ 1,2 =

σ xx + σ yy σ xx - σ yy 2 ) + σ 2xy , ± ( 2 2

viz. té Mohrovy krunice (co je vlastnì grafická konstrukce hledání hlavních smìrù napìtí pro speciální pøípad rovinné napjatosti). Hodnoty hlavních napìtí nás zajímají hlavnì proto, e jsou to velièiny nezávislé na volbì souøadného systému a mohou být tedy vodítkem pro hodnocení úrovnì napìtí. Na základì jistých hypotéz lze pomocí tìchto tøí hlavních napìtí vyjádøit jedinou velièinu, intenzitu napìtí, která pak rozhoduje napø. o tom zda pøi daném stavu napjatosti dojde k poruení materiálu nebo k jeho zplastizování. Protoe se ukazuje, e na tyto jevy nemá prakticky ádný vliv izotropní tlak (støední hodnota hlavních napìtí neboli tøetina prvního invariantu I), je výhodné izotropní tlak od tenzoru napìtí odeèíst a intenzitu napìtí vyjádøit pomocí deviátoru tenzoru napìtí

t t t 1 σ ′ = σ - Iσ δ . 3 Nejèastìji pouívané mìøítko intenzity je tzv. von Misesovo napìtí, definované prostøednictvím druhého invariantu tenzoru deviátoru napìtí (první invariant deviátoru je vdy roven nule):

σ M = 3.II σ ′ =

1 2 2 2 [( σ xx - σ yy ) + ( σ yy - σ zz ) + ( σ zz - σ xx ) + 6( σ 2xy + σ 2yz + σ 2xz )] 2 = σ 12 + σ 22 + σ 23 - σ 1 σ 2 - σ 1 σ 3 - σ 2 σ 3 .

Z poslednì uvedeného vyjádøení von Misesova napìtí je zøejmé, e pøi jednoosé napjatosti je jeho hodnota rovna jedinému nenulovému hlavnímu napìtí. Von Misesovo ekvivalentní napìtí lze interpretovat i tak, e je to napìtí pøi jednoosém tahu, kterému odpovídá stejná hustota deformaèní energie jako pøi daném stavu napjatosti (za pøedpokladu, e materiál je lineárnì elastický). Deformaci konstrukce popisujeme prostøednictvím vektoru posunutí

 u x  u     r u =  u y =  v ,      w  u z


31

který je vektorem v tenzorovém i maticovém smyslu (pøi maticových operacích ho chápeme jako sloupcový vektor). Zmìna posunutí jednotlivých bodù konstrukce podmiòuje její deformaci, která charakterizuje zmìnu vzdálenosti dvou blízkých bodù pøed a po zatíení (na deformaci nesmí mít vliv pohyb tìlesa jako celku /translace, rotace/). V pøípadì, e tyto zmìny nejsou pøíli velké, lze deformaci popsat tenzorem malých deformací, co je symetrická èást tenzoru gradientu posunutí (poznamenejme, e tyto úvahy jsou naprosto stejné jako v pøípadì tenzoru rychlosti deformace /viz. esták [22]/):

 e xx e xy  =  e yx e yy    ezx ezy

1 r t r e = ( ∆u + ∆ T u) = 2  1 ∂ u x ∂ uz  ∂ ux 1 ∂ ux ∂ u y ) + ( + ) (  2 ∂z ∂x  ∂x 2 ∂y ∂x    e xz   1 ∂u ∂u ∂ u y 1 ∂ u y ∂ uz  y x ) , + ) ( + e yz =  ( ∂x ∂y 2 ∂z ∂y    2 ∂y    ezz   ∂ uz   1 ( ∂ u x + ∂ uz ) 1 ( ∂ u y + ∂ uz )  2 ∂z ∂x ∂z  ∂y  2 ∂z

I pro tento tenzor lze nalézt tøi hlavní deformace a tøi smìry hlavních deformací zcela stejným zpùsobem jako u tenzoru napìtí. Analogicky jako intenzitu napìtí lze na základì druhého invariantu zavést i pojem efektivní deformace

ε ef =

2t t e:e = 3

2 2 [ e xx + e2yy + esub zz 2 + 2( e2xy + e2xz + e2yz )] , 3

která je potøebná pøi popisu elastoplastického chování materiálù. Definice (35) je navrena tak, aby pøi jednoosé plastické deformaci (kdy materiál povaujeme za nestlaèitelný, resp. rovna skuteèné deformaci ( e xx = ε ef

µ = 0.5 ) byla efektivní deformace (35)

, e yy = ezz = - ε ef ). 2 t r Podobnì jako byl tenzor napìtí σ nahrazen vektorem esti významných sloek σ , i tenzor deformace r t e je nahrazován estislokovým sloupcovým vektorem deformace ε


 ∂   ∂x   0    e xx     0  e yy    e    zz r r ε =   = ∆u =   ∂  2 e xy     ∂y  2 e yz     2 ezx  0     ∂   ∂z

32

0 ∂ ∂y 0 ∂ ∂x ∂ ∂z 0

0    0    ∂  ∂z   0    ∂  ∂y   ∂  ∂x 

kde poslední tøi sloky vyjadøují smykové deformace (oznaèované èasto

 u x   .  u y ,    u z

γ xy ,...). Proè byl do této "maticové"

definice zaveden (na rozdíl od definice "vektoru" napìtí) koeficient 2 u sloek smykových deformací? Vysvìtlení je vázáno na pojem Hustoty energie deformace. Objasnìme ho na jednorozmìrném pøípadì: tyè zatíená osovì silou F se prodlouí o ε. L ; práci, kterou bylo tøeba vykonat lze vyjádøit souèinem prodlouení a síly dìleno 2 (protoe pøi prodluování tyèe síla rostla lineárnì od nuly do hodnoty F). Je-li S plocha prùøezu bude energie vztaená na jednotku objemu rovna

W=

ε . L.F 1 = σ .ε [J.m-3 ] . Tato 2. L.S 2

definice hustoty deformaèní energie platí i v pøípadì vícerozmìrné napjatosti; protoe energie je skalární velièina, je ovem tøeba pouít dvojteèkový skalární souèin tenzoru napìtí a tenzoru deformace.

Pokud je vektor deformace definován dle (36) (tj. s koeficientem 2 u sloek smykových deformací) bude hustota deformaèní energie poèítaná tenzorovým zpùsobem souhlasit s výsledkem maticové operace (skalárním souèinem vektoru napìtí a vektoru deformace):

1 t t 1 W = σ :e = [ σ xx e xx + σ yy e yy + σ zz ezz + 2( σ xy e yx + σ yz ezy + σ xz ezx )] = 2 2 1 r r 1 = σ T ε = ( σ xx ε xx + σ yy ε yy + σ zz ε zz + σ xy γ xy + σ yz γ yz + σ xz γ xz ). 2 2

3.

Lineární statika Oblast lineární statiky je charakterizována lineární závislostí mezi tenzorem napìtí a deformací

(Hookùv zákon) a lineární závislostí mezi deformací (vyjádøenou tenzorem malých deformací) a posunutím; v tomto pøípadì je i napìtí lineární funkcí posunutí jednotlivých bodù konstrukce. Platí i princip superpozice, umoòující skládání posuvù a napìtí, která odpovídají samostatnì analyzovaným zatìovacím stavùm. Hookùv zákon lze pro izotropní materiál formulovat v tenzorové notaci


t σ =

33

E t Eµ r t e, ( ∆ .u)δ + 1+ µ (1+ µ )(1- 2 µ )

nebo jako vztah mezi deviátory tenzorù napìtí a deformace

t t σ ′ = 2Ge ′, G =

E . 2(1+ µ )

Pozn.: Vzpomeòte na konstitutivní rovnici newtonských nestlaèitelných kapalin (vztah mezi napìtím a rychlostí deformace je úplnì stejný jako (39), pouze místo modulu smyku G byla dynamická viskozita µ ). Tentý zákon lze formulovat i maticovì jako lineární vztah mezi vektorem napìtí a vektorem deformace

t r r r σ = D ( ε - ε 0 ),

t r D je matice elastických konstant (6 x 6 -její prvky závisí na E, µ ) a vektor ε 0 udává poèáteèní deformaci. Pro izotropní materiál a poèáteèní deformace vyvolané rozdílem teplot ( ε oxx = ε oyy = ε ozz = α∆T) kde

platí

µ µ  1- µ   µ 1- µ µ  µ 1- µ  µ   0 0  0 t E  D = (1+ µ )(1- 2 µ )   0 0  0    0 0  0 

0

0

0

0

0

0

1- 2 µ 2

0

0

1- 2 µ 2

0

0

0  0   1 0     1    0 t α∆TE  1  , D ε0 = . 1- 2 µ  0     0  0     0   1- 2 µ  2 

Metoda koneèných prvkù se opírá o princip virtuálních posunutí (øíkající, e práce vnitøních sil pøi virtuální deformaci musí být rovna práci vnìjích sil /objemových i povrchových/). Tento princip je zcela obecný (vyplývá pouze z Cauchyho rovnic rovnováhy) a pouívá se i v nelineární oblasti (napø. mimo oblast platnosti Hookova zákona). Pokud vak Hookùv zákon platí (tj. pokud je mezi vnitøními silami a deformací lineární vztah) dá se z principu virtuálních posunutí odvodit Lagrangeùv variaèní princip minima celkové potenciální energie:

Ze vech kinematicky pøípustných stavù (charakterizovaných spojitým posunutím u(x)) nastává ten, který má nejmení celkovou potenciální energii (souèet deformaèní energie a energie vnìjích sil).


34

Vycházejíce z definice hustoty energie deformace (37), Hookova zákona (40) a ze vztahu mezi deformací a posunutím (36), mùeme deformaèní energii vyjádøit integrálem, který závisí pouze na posunutích

r u vyetøované konstrukce

r E i (u) =

∫

V

1 1 rT r ε σdV = 2 2

t

r r rT D( ε -ε 0 ) ε ∫

V

dV =

1 ∫ 2

r T t r r ( ∆u ) D( ∆u - ε 0 ) dV,

V

a podobnì lze formulovat i potenciální energii vnìjích sil (jako souèin posunutí a síly):

rT r r E e (u) = - ∫ f u V

kde

r dV - ∫ pr u dS, T

Sp

r r f je objemová síla (napø. tedy gravitaèní zrychlení krát hustota pøi uvaování zatíení vlastní vahou) a p

je pøedepsané napìtí na povrchu (napø. vnìjí tlak). Z pøedchozích výrazù je patrné, e deformaèní energie je kvadratickou funkcí posunutí, zatímco potenciální energie vnìjích sil lineární funkcí posunutí a jejich souètem je parabola, viz. následující obrázek

1 r E L (u) = 2

∫

V

1 r t r rT r r T t r T r ( ∆u ) D ∆u dV - ∫ ( ε 0T D ∆u + f u) dV - ∫ p u dS. 2 V Sp Cílem øeení je nalézt takové posunutí konstrukce

r r u(x) , které by bylo

"kinematicky pøípustné" (spojité a vyhovující pøedepsaným okrajovým podmínkám) a zároveò minimalizovalo celkovou potenciální energii

r E L (u) . V

MKP se navrhne mnoina posunutí a z ní se vybere to nejlepí; pøesnìjím vyjádøením této mylenky je Ritzova metoda, spoèívající v aproximaci

r r u(x) lineární kombinací jistých základních (a r víceménì libovolnì navrených) tvarù N i (x) (bázových nebo té skuteèné funkce posunutí

interpolaèních èi tvarových funkcí), v maticové vyjádøení

r r ( ui (x) = ∑ N ij (x) q j

t r r r r u(x) = N(x)q,

N

j=1

), kde

odpovídající tøem smìrùm posunutí Sloupcový vektor

t r N(x) je matice tvarových funkcí, která má tøi øádky

r r u(x) a N sloupcù, z nich kadý odpovídá jednomu "pokusnému tvaru".

r q má N prvkù, koeficientù zmínìné lineární kombinace tvarových funkcí. Tyto koeficienty

r q mohou mít zcela reálný význam, mohou to být napøíklad pøímo posunutí nebo natoèení uzlových bodù, proto r se q nazývá vektor zobecnìných posunutí. Dosadíme-li interpolovaná posunutí (45) do (44), pøejde Lagrangeùv funkcionál na kvadratickou

funkci parametrù

1 r E L (q) = 2

r q

tr T t tr ∫ ( ∆ Nq ) D( ∆ Nq) dV -

V

t r rT t r 1 rT t D ( Nq)+ ( ∆ f Nq) dV ε 0 ∫ 2 V

∫

Sp

r pT

tr Nq dS,


35

r

a protoe vektor zobecnìných posuvù není funkcí x mùeme ho vytknout pøed integrály (na prostorové souøadnici závisí pøedevím tvarové funkce N, zadávaná zatíení a nìkdy té matice elastických konstant

r ε0) t t t t rT t 1r t 1 rT r r q [ ∫ ( ∆ N )T D( ∆ N) dV] q - [ ∫ ( ε 0T D ( ∆ N)+ f N) dV + E L (q) = 2 2 V V

nebo poèáteèní deformace

1 = 2

r qT

tT t t [ ∫ B DB V

r 1r dV] q - [ ∫ ( 2 ε

Nalézt v této fázi øeení vektor

V

T 0

t t rT t DB+ f N)dV +

∫

Sp

∫

Sp

r pT

t D

t r N dS] q =

r t pT N dS] qr =

r r 1 rT t r = q K q - F T q. 2

r q , který minimalizuje kvadratickou funkci (47) je u snadné; staèí

poloit vechny první parciální derivace funkce (47) rovny nule, co vede na soustavu lineárních rovnic

t r rr r Kq = F, kde K je matice tuhosti konstrukce (ètvercová, symetrická rozmìru N x N) a F vektor zatíení respektující teplotní zatíení (poèáteèní deformaci), a objemové i povrchové síly. Matici tuhosti i vektor zatíení lze jak patrno ze (47) získat sèítáním integrálù pøes jednotlivé elementy,

t r K, F jediného obecného prvku: Je tøeba ukázat jak v jednotlivých t typech prvkù vypadají matice diferenciálních operátorù ∆ , elastických konstant D a jak volit tvarové funkce t r N(x) . Existuje velmi mnoho zpùsobù volby tvarových funkcí a zobecnìných posunutí, take existují stovky z

take se staèí zabývat jenom výpoètem

vnìjího pohledu stejných koneèných prvkù, které se vak mohou chovat rùznì (za uèitých okolností a patologicky). Mnohost typù koneèných prvkù je zapøíèinìna i tím, e Lagrangeùv princip minima celkové potenciální energie není jediný (uvádíme ho proto, e v praxi pøevauje a program COSMOS/M ho pouívá výluènì); jeho protipólem je Castiglianùv princip minima doplòkové (komplementární) energie, kdy hledanými velièinami nejsou posuvy ale napìtí (proto se také hovoøí o deformaèní /Lagrange/ nebo o silové metodì /Castigliano/ øeení problémù lineární prunosti). Prvky, zaloené na Lagrangeovì principu, se vyznaèují tím, e náhradní koneènìprvková konstrukce je tuí ne skuteèná, zatímco u Castiglianova principu je tomu naopak; tento teoretický výsledek vak ve skuteènosti vdy neplatí, protoe matice tuhosti se nepoèítají vdy pøesnì podle pøedpisu (47) (a napø. numerická integrace prvkù matice tuhosti zpùsobuje sníení jejich tuhosti). A ji je vak pouit jakýkoliv variaèní princip, vdy musí platit, e pøi zmenování rozmìrù prvkù se øeení blíí ke skuteènosti (pøesnìji, pøedpokládanému modelu kontinua). Pozn.: Je rozumné provádìt analýzu konkrétní konstrukce vdy alespoò na dvou rùzných sítích; rozdíl výsledkù je základním vodítkem pro odhad nepøesnosti. To, které z tìchto dvou získaných øeení je lepí plyne pøímo z pouitého variaèního principu: u deformaèní metody to bude øeení s mení celkovou potenciální energií a proto je tato hodnota uvádìna ve výstupních sestavách programù.


3.1

36

Jednorozmìrné prvky (nosníky, potrubní prvky) Pod pojmem jednorozmìrný prvek máme na mysli speciální pøípad, kdy sice platí ve co bylo øeèeno o

obecné trojrozmìrné napjatosti (vèetnì vyjádøení matice tuhosti a zatíení), ale geometrie konstrukce umoòuje nahradit objemové integrály integrály køivkovými (na základì vhodných pøedpokladù o vazbì tøí sloek posunutí). Do této skupiny patøí vechny nosníkové a potrubní prvky. 3.1.1

Analytická øeení: køivé pruty (Castigliano), tah, ohyb, krut pøímých nosníkù a trubek, zvlátní pøípady: koleno, prstence. Nejdùleitìjí aplikace moných analytických øeení se týkají dostateènì jednoduchých pøíhradových

konstrukcí, výpoètu potrubních úsekù zatíených silami, momenty, dilatací a koneènì prstence (prstencové výztuhy tlakových nádob, ale té napø. i horizontální nádre s promìnným obvodovým zatíením, pokud lze ovem zanedbat vliv podélných deformací), Køupka [12]. S výjimkou pøíhradových konstrukcí lze vechny tyto pøípady povaovat za køivé pruty a jejich øeení hledat silovou metodou (viz. zmínìný Castiglianùv princip), co je v daném kontextu jednoduí ne deformaèní metoda. Napø. u potrubního úseku, který je na jednom konci zatíený známými silami je velice snadné spoèítat prùbìhy vnitøních sil a momentù podél støednice potrubí a tedy i komplentární energii vnitøních sil. Z první Castiglianovy vìty pak plyne, e posunutí pùsobitì síly lze stanovit derivací komplementární energie dle pøísluné síly (dtto pro zatíení momenty; poznamenejme, e to ve je jenom opakování základního kurzu prunosti a pevnosti [31-32]). Spoèítat prùhyby u staticky urèité prutové konstrukce je tedy jednoduché; staticky neurèitý pøípad (tøeba potrubní úsek, ukotvený na obou koncích a podrobený dilataci) se øeí uvolòováním; úloha se tak pøevede na staticky urèitou a staticky neurèité velièiny (síly/momenty v místì uvolnìní) se vypoètou z deformaèním podmínek. Pro øeení bìných úloh tedy staèí minimální teoretická výbava; staèí vìdìt jak se spoèítá komplementární energie napjatosti køivého prutu kdy je v kadém místì známa axiální síla F, posouvající síla T, ohybový a kroutící moment

M0 , M k :

1 β T 2 M 2k 1 1 F 2 M 02 ∫ + )] dl, ( ) + + ( [ E = G S J 2 E S Jp *

S [ m2 ] je plocha prùøezu, J [ m4 ] moment setrvaènosti prùøezu vzhledem k ose pùsobícího ohybového 4 momentu a J p [ m ] polární moment setrvaènosti. Souèinitel β je korekcí na skuteèný prùbìh smykového napìtí po prùøezu ( β = 1 by odpovídalo konstantnímu smykovému napìtí, co nemùe být pravdou ji jen kde

proto, e na povrchu musí být toto napìtí nulové /viz. zákon o sdruených napìtích/). Pro nejfrekventovanìjí prùøezy jsou hodnoty momentù setrvaènosti uvedeny v následující tabulce


37

J

Jp

β

1 3 bh 12

bh 2 2 (b + h ) 12

1.2

π d4 π ( D4 - d 4 ) 64

π d4 π ( D4 - d 4 ) 32

1.185

Metodiku pouití Castiglianova principu na výpoèet potrubních úsekù popisuje napø. Jao [6], zamìøíme se proto jen na nìkteré speciální pøípady, které souvisí s návrhem dùleitých potrubních systémù, kdy je tøeba uvaovat nejen vliv dilatací, osamìlých sil a pøetlaku (známé vztahy pro tangenciální a osová

σt =

membránová napìtí

pr pr , σ0 = ) nýbr i geometrických odchylek, daných buï nerovnomìrností 2s s

tlouky stìny (typické pro ohýbaná hladká kolena) nebo nekruhovitostí prùøezu (ovalita potrubí se pøi pùsobení vnitøního pøetlaku projeví vznikem ohybových napìtí). Tyto eventuality se stále jetì dají zvládnout na základì metodiky výpoètu køivých prutù, jak ukáeme na pøíkladì: Uvaujme potrubí ètvercového prùøezu zatíené vnitøním pøetlakem p, viz. následující obrázek.

Základním úkolem je stanovit pøídavnou ohybovou napjatost v kritickém místì, které je v polovinì strany; vzhledem k symetrii by staèilo uvaovat jen osminu prùøezu (a chápat stìnu jako rovný nosník), abychom se vak pøiblíili zpùsobu výpoètu tøeba eliptické trubky, nahradíme jen jednu ètvrtinu prùøezu køivým prutem, který je na jednom konci vetknutý a na druhém zatíený silou F=pH/2 a neznámým momentem

MO .

Zanedbáme-li energii odpovídající smykovým a osovým silám, zbude jen komplementární energie ohybové napjatosti vyjádøená souètem integrálù H/2 1  x2 2 E =  ( M o - p ) dx + 2EJ  ∫0 2 * i

H/2

∫ 0

H p H2 2 2  ( M o + py - ( - y ) ) dy  2 2 4 

Úhel pootoèení prùøezu v místì pùsobení hledaného ohybového momentu vyjádøíme z první Castiglianovy vìty a poloíme rovný nule (podmínka symetrie deformovaného prùøezu)

1  ∂ E *i ϕ= =  ∂ M o EJ 

H/2

∫ 0

x2 ( M o - p )dx + 2

H/2

∫ 0

( M o + py

 H p H2 2 - y ))dy  = - ( 2 2 4 

3

= Výsledkem

je

pH 1 1 ( Mo H ) ≡ 0 _ Mo = p H2 . EJ 24 24

ohybový

moment

M o , jemu odpovídá maximální pøídavné ohybové napìtí


38

pH H M o p H2 = (1+ ) . σt = 2 , resp. maximální celkové napìtí (membránové + ohybové) σ t = 4s 2s 2s Wo Nìkteré potrubní prvky by se správnì mìly poèítat jako skoøepiny (viz. následující kapitola) a ne jako køivé pruty. Nejdùleitìjím pøípadem jsou kolena (a ji hladká nebo segmentová) jejich analýzu provedl T. Kármán ji v roce 1911. Ukázal, e koleno, by s pøesnou geometrií anuloidu je podstatnì mìkèí (má mení ohybovou tuhost) ne nosník se stejným mezikruhovým prùøezem (a ohybovou tuhostí E.J; pomìr ohybové tuhosti pøímé trubky a kolena se nazývá Kármánovo èíslo k). Ani prùbìh napìtí po obvodu prùøezu neodpovídá pøedstavì ohýbaného nosníku, protoe maximální napìtí nejsou na vnìjím a vnitøním obvodu kolena a pøedevím jsou vìtí ne napìtí poèítaná z momentu odporu mezikruhového prùøezu. Podstané je vìdìt, e tyto odchylky jsou závislé na pomìru geometrického prùmìru tlouku stìny a polomìru køivosti k polomìru trubky a napø. pro Kármánovo èíslo platí

k=

kde R je polomìr køivosti kolena, r,s polomìr a síla stìny trubky. Pro Rs 10 + 12 λ 2 , = 0.35 , λ ≥ pøedstavu: koleno s polomìrem køivosti R=1 m, polomìru r=0.1 m a 1+ 12 λ 2 r2 sílou stìny s=1 cm má λ = 1 a jeho tuhost se zmení cca 1.7 krát.

3.1.2

Bernoulliho a Timoenkùv nosník, odvození matice tuhosti a vektoru zatíení V MKP o vlastnostech modelu rozhodují pøedevím matice tuhosti a proto se budeme vìnovat právì

jim. Obecný postup odvozování matice tuhosti (rovnice (42) a dalí) uplatníme na jednorozmìrný nosníkový prvek konstantního (ale jinak libovolného) prùøezu, který se bude deformovat pouze v jedné rovinì xz (neuvaujeme krut ani axiální zatíení), viz. obr. 13 V tomto pøípadì zanedbáváme napìtí esti sloek zùstávají jen dvì

σ yy = σ zz = σ zy = σ xy = 0 , take ze

σ xx , σ xz a výraz pro deformaèní energii nosníku

se zjednoduí: L

1 Ei = ∫ 2 x=0

∫ ∫ (σ

xx

ε xx + σ xz γ xz )dzdydx =

y z

1 2 ∫ ∫ ∫ (E ε 2xx + G γ xz )dzdydx = 2 ∂u ∂w 2  1  ∂u 2 ) + G( + ) dzdydx. = ∫ ∫ ∫  E( ∂z ∂x  2  ∂x =


39

kde u je posunutí ve smìru osy nosníku x, zatímco w je prùhyb (posunutí ve smìru z). Poznámka: Ve výrazu (53) je pouit Hookùv zákon ve tvaru σ xx = E ε xx , co se zdá být v rozporu s univerzálnì platným vyjádøením Hookova zákona (40,41):

E [(1- µ ) ε xx + µ( ε yy + ε zz )] . ádný rozpor zde vak není, protoe pøi jednoosé (1+ µ )(1- 2 µ ) napjatosti není ε yy + ε zz rovno nule, nýbr - 2 µ ε xx , co vyplývá z pøedpokladu o nulových napìtích ve σ xx =

smìrech kolmých na osu nosníku (ovìøte na základì (40,41)). Pro vystiení prùbìhu axiálního posunutí u(x,z) se veobecnì pøijímá pøedpoklad, e prùøez nosníku zùstává i po deformaci rovinný a pøi malých deformacích (malých úhlech natoèení) tedy platí, e

u(x,z) = zψ (x), kde ψ je úhel natoèení roviny prùøezu (vzhledem k poloze v nezatíeném stavu). Protoe ani prùhyb nosníku w(x) nezávisí na souøadnici z, lze provést naznaèenou integraci pøes celý prùøez nosníku

1 Ei = 2

 ∂ψ 2 GS ∂w )+ (ψ +  EJ( ∂x ∂x β x=0  L

∫

 )2  dx, 

kde J je moment setrvaènosti vzhledem k neutrálné ose y a S plocha prùøezu. Nyní jsou moné dva základnì odliné pøístupy; klasická Bernoulliho-Navierova teorie vychází z pøedpokladu, e rovina prùøezu se po zatíení nezkroutí a navíc zùstane kolmou k (deformované) støednici nosníku. Tento pøedpoklad znaèí, e pootoèení prùøezu ψ není nezávislé na prùhybu w , nýbr platí

∂w ψ = , ∂x

L

Ei =

∫

x=0

EJ(

∂2 w 2 ) dx. ∂ x2

Uvedená hypotéza vede jak patrno k zanedbání vlivu smykových sil, co je pøijatelné u relativnì tíhlých nosníkù. Podotknìme, e tentý pøedpoklad (zachování kolmosti prùøezu k ose) se pouívá i v teorii tzv. Kirchhoffových desek. Dalím krokem odvozování matice tuhosti je volba aproximace prùhybové èáry

w(x) ; vodítkem mùe

být analytické øeení nìkterých jednoduchých pøípadù - napø. prùhyb nosníku, zatíeného osamìlou silou je kubický polynom a ten se dá pouít i v tomto pøípadì. Obecný kubický polynom má 4 koeficienty, které musí být jednoznaènì urèeny ètyømi uzlovými parametry: u dvouuzlového nosníkového prvku hodnotami uzlových

w1 ,w2 a natoèení w1′ ,w2′ (tím bude automaticky zajitìna návaznost prùhybové èáry i na dalí nosníky, pøipojené k uzlùm 1 a 2). Prùhyb w(x) vyjádøíme jako lineární kombinaci ètyø tvarových funkcí N i (x) (kubických polynomù, splòující v uzlech 1,2 podmínky, znázornìné na následujícím obrázku; napø. N 1 (x) je kubický polynom, s jednotkovým posunutím v uzlu 1, nulovým posunutím v uzlu 2 a nulovou

prùhybù

derivací v obou uzlech):

w(x) = w1 N 1 (x)+ w1′ N 2 (x)+ w2 N 3 (x)+ w2′ N 4 (x), x 2 2x x 2 N 1 (x) = (1- ) (1+ ), N 2 (x) = (1- ) x, ... L L L


40

Po dosazení aproximace prùhybu (57) do výrazu pro deformaèní energii (56) a jeho pøepisu do maticového tvaru, zjistíme, e prvky matice tuhosti jsou vlastnì integrály ze souèinù druhých derivací tvarových funkcí

 w1   1  w1′   Ei = 2  w2     w2 ′ 

T

 ∂2 N 1   2   ∂x     2   ∂ N2 2  L  ∂x   2 ∂ N1 ∫x=0 EJ  2   ∂ x2  ∂ N 3  ∂ x2     2   ∂ N 4    ∂ x2 

2

∂ N2 ∂ x2

2

∂ N3 ∂ x2

 w1   L 2 2 2  w1′ 1 rT t r  ∂ N4 ∂ Ni ∂ N j    Kq =  = dx EJ q  K ij ∫0 ∂ x2 ∂ x2 dx . ∂ x2  2   w2    w2′

Druhý pøístup (v pøípadì nosníkù je jeho autorství pøisuzováno Timoenkovi) povauje obì velièiny

ψ ,w za nezávislé zobecnìné prùhyby. Výraz pro deformaèní energii (55) je sice ponìkud sloitìjí ne (56), ale lze pouít jednoduí aproximaèní funkce zobecnìných posuvù, protoe v integrandu se vyskytuje maximálnì první derivace (u Bernoulliho nosníku druhá derivace prùhybu). V tomto pøípadì lze pouít dokonce i lineární interpolaèní (tvarové) funkce (i kdy je jasné, e skuteèný nosník se takto deformovat nemùe - lineární závislost prùhybu w(x) odpovídá pouze nezatíenému nosníku):

w(x) = w1 N 1 (x)+ w2 N 2 (x), ψ (x) = ψ 1 N 1 (x)+ ψ 2 N 2 (x), N 1 (x) = 1kde

x x , N 2 (x) = , L L

N 1 (x), N 2 (x) jsou tvarové funkce uzlù 1, 2 a koeficienty ψ 1 ,w1 ,ψ 2 ,w2 tvoøí vektor uzlových

parametrù (poootoèení a posuv uzlù 1,2). Aproximaci zobecnìných prùhybù (59) dosadíme do vztahu pro deformaèní energii (55), který pøepíeme do pùvodní maticové notace a získáme tak tuhostní matici Timoenkova nosníku


 ψ 1   1  w1   = Ei 2 ψ  2    w2 

T

L

∫

x=0

 ∂ N1   ∂x   0     ∂ N2  ∂x   0   

N 1    ∂ N 1  EJ ∂x    0 N 2     ∂ N2  ∂x 

41

 ∂ N1 0    ∂x  GS    N1 β  

0

∂ N2 ∂x

∂ N1 ∂x

N2

 ψ 1  K 11 K 12 K 13 K 14     . . 1  w1  K 21 K 22 =   .  2 ψ   . . . .  2     . . K 44  w2   .

 ψ 1    w1 .   ψ 2    w2 

T

 ψ 1 0     w1   = dx .   ψ 2 ∂ N2    ∂x   w2 

Stejnì jako v pøedchozím pøípadì Bernoulliho nosníku jsou prvky matice tuhosti integrály souèinù tvarových funkcí, které lze vyjádøit analyticky (integrandem jsou polynomy). Napøíklad diagonální prvek

K 33

vypoèteme takto (ovìøte roznásobením matic v (60)) L

K 33 = ∫ [EJ( 0

EJ GS L dN 2 2 GS 2 )+ + , N 2 ]dx = β β 3 dx L

( N 2 (x) =

x ). L

Matice tuhosti získaná tímto zpùsobem, tj. analytickou integrací mùe za jistých okolností zpùsobit naprosto nekorektní chování prvku. Pouijeme-li ji napø. pro výpoèet prùhybu vetknutého nosníku, který je na volném konci zatíen silou F dostaneme následující vztah

    3 3 Lβ L  , w2,p_esnÇ _e_enÍ = F  Lβ + L  .  + w2 = F  GS GSL2   GS 3EJ  4EJ +   3β   Srovnání s "pøesným" øeením z bìné teorie nosníkù ukazuje, e u dlouhých nosníkù by parazitní druhý èlen ve jmenovateli (62) mohl dokonce pøeváit nad ohybovou tuhostí EJ (pøièem by výsledná tuhost nosníku podstatnì vzrostla a pøestala záviset na momentu setrvaènosti prùøezu!). Tento jev je známý pod jménem "locking" (zablokování) a souvisí s ponìkud diskutabilní volbou lineárních tvarových funkcí (u kubických tvarových funkcí, pouitých u Bernoulliho nosníku by se to stát nemohlo). Uvedený Timoenkovský prvek je moné podstatnì vylepit pouitím numerické, namísto pøesné analytické integrace; v MKP se témìø výhradnì pouívají tzv. Gaussovy integraèní formule:


1

∫

f(x)dx =

-1

kde

k

42

ζ i a gi jsou hodnoty tabelované pro rùzné stupnì formulí k (napø. Gaussovì formuli /k=1/ odpovídá g1 = 2, ζ 1 = 0 , co lze

∑ g f( ζ ) jednobodové i=1

i

i

interpretovat tak, e hodnota integrálu je aproximována souèinem délky

integraèního intervalu /2/ a hodnoty integrandu uprostøed tohoto intervalu /f(0)/). Pouijeme-li jednobodovou Gaussovu formuli, získáme integrací (61) napø. tuhost L

K 33 = ∫ [EJ( 0

K 33 ,

EJ GS L dN 2 2 GS 2 )+ , + N 2 ]dx = L dx β β 4

( N 2 (x) =

x ), L

která je mení ne tuhost (61) poèítaná analyticky (numerická integrace se projevila pouze zmenením smykové tuhosti). Co je vak dùleitìjí: jednobodová Gaussova integrace úplnì zruí parazitní smykovou tuhost (62), která byla vlastní pøíèinou zablokování prvku. Matici tuhosti Timoenkovského nosníkového prvku, získanou jednobodovou Gaussovou integrací, lze posléze zapsat jako souèet matice ohybové a matice smykové tuhosti:

1 0  t EJ  . 0 K = L . .  . .

1 1 1 1   2L 4 2L  4     . 1 1 1 - 1 0 - 2  2  2L L  L 0 0 GSL  + .  β  1 0 . 1 1  .    . 0 4 2L      . . 1  .    L2 

kde J je moment setrvaènosti a S plocha prùøezu. L je délka elementu a

β korekce na rozloení smykových

napìtí. Vektor pravé strany je tvoøen ètveøicí hodnot - zatíení dvou uzlù analyzovaného prvku. Kontrolní otázka: Co tato zatíení pøedstavují? 3.1.3

TRUSS, BEAM, PIPE, ELBOW Program COSMOS/M pouívá pro modelování jednorozmìrných prvkù typu táhlo (TRUSS2D a

TRUSS3D) prvky odvozené ji v kapitole 2, a pro modelování nosníkù (BEAM2D, BEAM3D) a pøímých trubek (PIPE) prvky, zaloené na kubických tvarových funkcích pro posunutí a kvadratických funkcích pro aproximaci prùbìhù natoèení. O tom, jakým zpùsobem je modelován prvek ELBOW (koleno) nebo TPIPE (odboèka), firma ádné blií informace neposkytuje (Kármánovské sníení tuhosti kolena vak respektováno je). Jako pøíklad zadávání parametrù uvedeme prvek ELBOW: 3 uzlové body, první dva totoné s koncovými prùøezy, tøetí uzel støed køivosti. Kadý z koncových uzlù má 6 stupòù volnosti (3 posuvy a 3


43

natoèení). RCONST: zadávané konstanty R1 - vnìjí prùmìr trubky R2 - tlouka stìny R3 - vnitøní pøetlak R4 - polomìr køivosti MPROP: materiálové parametry EX - modul prunosti ALPX - souèinitel teplotní roztanosti NUX - µ DENS- ρ . Prvek nelze pouít pro nelineární analýzu (tøeba pro uvaování plasticity), popis nosníkových prvkù viz. pøíloha.

3.2

Dvourozmìrné prvky (skoøepiny, prstencové prvky) Do této skupiny náleí pøedevím ty konstrukèní uzly, které je tøeba vyetøovat jako skoøepiny

(odboèky potrubí, plátì nádob, støení konstrukce), ale té tìlesa rotaènì symetrická (høídele, pøíruby). 3.2.1

Analytická øeení membrán a skoøepin Skoøepinové konstrukce mohou být za jistých okolností pokládány za membrány a pak není nutné

uvaovat ohybová napìtí (membrána je schopna pøenáet pouze teèné síly). Jsou to pøedevím válcové, kuelové a kulové skoøepiny s konstantní tloukou stìny vystavené pøetlaku (osamìlé síly jsou vylouèeny). Jen zcela výjimeènì lze dosáhnout membránové napjatosti v místech uloení nebo skokových zmìn geometrie skoøepiny (napø. odboèky potrubí); pozitivním pøíkladem je prùnik dvou kulových nádob, zesílený pøesnì vypoèteným ebrem (ebro, které by bylo mìkèí nebo tuí ne optimální vyvolá okamitì poruchu membránové napjatosti v okolí prùniku koulí). Podotknìme, e napø. v ohnuté trubce nebo vlnovci nemùe mebránová napjatost nastat nikdy, co vak neznamená, e by "membránový" výpoèet nebyl vhodný pro získání alespoò orientaèních hodnot napìtí (bìná taktika spoèívá v tom, e se plá nádoby nahradí ve vyetøovaném místì válcovou plochou, a ta se spoèítá jako membrána). Výpoèet mebrán a zvlátì rotaènì symetrických a symetricky zatíených mebrán je pomìrnì snadný, protoe lze pouít silovou metodu (tj. zjistit rozloení vnitøních sil jen z rovnic rovnováhy). U membrán a skoøepin se místo napìtí (Pa) hodnotí vnitøní síly velièinami, které jsou vztaené na jednotku obvodu skoøepiny: N

- normálové síly [N/m], pùsobící kolmo k mylenému øezu skoøepiny, ale v její rovinì

S

- smykové síly [N/m], pùsobící ve smìru øezu a v rovinì skoøepiny

Pouze u skoøepin pak mohou pùsobit T

- posouvající síly [N/m], pùsobící kolmo k rovinì skoøepiny

M

- ohybové a kroutící momenty [N.m/m]. U rotaènì symetrických membrán uvaujeme na ploe skoøepiny tyto dvì sloky normálových sil:

jednu ve smìru meridiánu a druhou ve smìru obvodovém (jsou to toti smìry dvou hlavních napìtí). Dvì


44

hlavní napìtí se dají vypoèíst ze dvou rovnic rovnováhy: ve smìru osy membrány (staèí uvolnit napø. vrchlík kulové skoøepiny a zjistit výslednici vnìjích sil, která naò pùsobí) a ve smìru kolmém k povrchu membrány. Tato druhá rovnice rovnováhy sil se nazývá Laplaceova a má tvar (odvození uvádí napø. Jao [6]):

Nα Nϕ + = p, Rα Rϕ v nìm symboly

Rα , Rϕ oznaèují tzv. hlavní polomìry køivosti mebrány (viz. obr.15).

Naznaèeným postupem lze najít napø. vztahy pro hlavní napìtí v kulové skoøepinì ( σ

= N / s = pr / 2s ) nebo trubce, zatíené vnitøním pøetlakem (tzv. "kotlové" vzoreèky σ t = pr / s, σ o = pr / 2s ). Znaènì sloitìjí je analýza "opravdových" skoøepin (tj. vlivu ohybových napìtí). Nejdùleitìjí jsou dva pøípady: Kruhová deska a Válcová skoøepina. Pro kruhové desky existuje øada analytických øeení (viz napø. ubrt [30]), která odpovídají rùzným pøedpokladùm o tvaru deformace prùøezu, resp. zanedbání urèitých typù vnitøních sil. Napøíklad tzv. Kirchoffova deska vychází z hypotézy o zachování kolmosti prùøezu k deformaèní ploe a zanedbává vliv membránových sil (co je ekvivalentní pøedpokladùm, které byly pøijaty pøi øeení Bernoulliho nosníku). Kirchhoffova teorie je pouitelná pro relativnì tenké desky (jinak by se uplatnily i smykové síly) a malé prùhyby (kdy se jetì neprojevují síly membránové). Umoòuje odvodit napø. následující vztahy pro prùhyb a napìtí v desce, zatíené obvodovým momentem:

w(r) =

6 MR2 (1- µ )  r 1- ( 3  R Es 

6M  )2  , σ t = σ r = 2 . s 

Pøi zatíení tlakem popisuje prùhybovou èáru Kirchhoffovy desky diferenciální rovnice 2 6(1- µ ) d 1 d dw pr, )) = (r ( E s3 dr r dr dr

kterou lze pro rùzné typy uloení okrajù desky integrovat analyticky. Napìtí, pùsobící na povrchu desky, lze pak stanovit z následujících vztahù


σt = -

45 2

Es 1 dw d w +µ 2 ) 2 ( dr 2(1- µ ) r dr

σr = -

Es d 2 w µ dw + ). ( 2 2(1- µ ) dr 2 r dr

Poznamenejme, e vechny tyto výsledky jsou pouze opakováním látky základního kurzu prunosti a pevnosti. V øadì pøípadù vak Kirchhoffova teorie nevyhoví a vede ke znaènému pøedimenzování desek; je toti omezena pøedpokladem malých prùhybù w/s < 0,3 (kdy zatíení pøenáejí pouze ohybové momenty). Souèasné uvaování ohybové i membránové tuhosti vak nutnì vede k nelineární teorii velkých posunutí (mebránová napìtí se projeví jen pøi velkých prùhybech desky), která bude nastínìna a v následující kapitole. Zde uveïme alespoò základní výsledky: Maximální prùhyb desky, zatíené konstantním tlakem je dán kubickou rovnicí její koeficienty α 1 ,α 2 jsou tabelovány pro rùzné typy uloení okrajù p R w0 w0 3 ) = α 2 ( )4 , +α1( desky. Maximálnímu prùhybu vypoètenému z rovnice (70) odpovídají E s s s tato mebránová napìtí σ m v ose desky (jsou úmìrná druhé mocninì prùhybu) a ohybová napìtí v ose nebo na okraji desky (ohybová napìtí jsou úmìrná první mocninì prùhybu):

σ m = α 3 E(

w0 2 ) , R

pøièem koeficienty

σ o = α4 E

w0 s , R2

α i jsou uvedeny v následující tabulce (ubrt [30]) α1

α2

α 3 støed

α 4 støed

α 4 okraj

pevné vetknutí

0.528

0.171

0.976

2.86

4.4

volné podepøení

0.262

0.696

0.295

1.778

0

Válcovou skoøepinu (trubku) s rotaènì symetrickým zatíením (a ji pøetlakem p nebo okrajovým momentem èi silou) je moné nahradit elasticky vázanou soustavou podélných proukù (pásových elementù, nosníèkù). Na základì této pøedstavy lze prùhybovou èáru (radiální posun stìny trubky v závislosti na souøadnici x ve smìru osy trubky) vyjádøit diferenciální rovnicí ètvrtého øádu, která popisuje prùhyb nosníku uloeného na pruném podkladì (Køupka [12]): 2 4 p* 3(1- µ ) d4 w Es3 4 , +4β w = , K= 2 , β= K dx 4 Rs 12(1- µ )

kde K je ohybová tuhost nosníèku s obdélníkovým prùøezem (jednotkové íøky a výky s) a

p* je ekvivalentní

spojité zatíení (pøi absenci teplotních napìtí a osových sil je to pøímo tlak p). Pozn.: "Pruný podklad" nosníèkù zohledòuje vazbu mezi fiktivními prouky, na nì jsme trubku "rozstøíhali".

β , jeho hodnota klesá s geometrickým prùmìrem polomìru trubky R a její tlouky s (pro nekoneènì velké R je β = 0 , a deformaci stìny trubky lze poèítat z

Míru vzájemného ovlivnìní proukù udává parametr "klasické" rovnice volného nosníku).

Protoe rovnice (72) je lineární, dá se øeit analyticky (viz. Køupka [12]) a výsledkem jsou prùhybové


46

èáry typu

w(x) = e- βx g( βx)+ w0 , kde funkce g je lineární kombinací sinù a kosinù s argumentem βx ; vystihuje prùbìh vln prùhybù, deformací a napìtí, které jsou vyvolány místní poruchou membránové napjatosti. Typickým pøíkladem je spojení dvou trubek s rùznou tloukou stìny, zatíených vnitøním tlakem: spojitost posunutí i natoèení musí zajistit ohybový moment a posouvající síla v místì styku. Poruchová napìtí se projevují jen v malé vzdálenosti od místa diskontinuity a tuto vzdálenost charakterizuje parametr 3.2.2

β [ m-1 ] ; vyplývá to pøímo z øeení (73), kde tlumení poruch popisuje èlen e- βx . Koneèné prvky pro aproximaci membrán, desek a skoøepin Uvedeným pøípadùm kruhových desek a symetricky zatíených válcových skoøepin odpovídají de facto

jednorozmìrné koneèné prvky - úseèky se dvìma uzly, které pøedstavují meridián rotaènì symetrické kuelové skoøepiny (deska i válec jsou speciální pøípady). Tyto prvky mají mnoho spoleèného s døíve analyzovanými nosníky; analogií Bernoulliho nosníku je napø. prvek pouívaný pod názvem SHELLAX. Uzlové parametry jeho dvou uzlù pøedstavuje dvojice posunutí u,w a natoèení prùøezu

ψ . Prùbìh posunutí ve smìru normály povrchu skoøepiny w(s) se modeluje kubickým polynomem (jeho 4 koeficienty jsou jednoznaènì urèeny posuvy

w1 ,w2 a natoèením prùøezù ψ 1 = - w1′ , ψ 2 = - w2′ ) a posunutí ve smìru meridiánu u(s) polynomem lineárním (2 koeficienty lineárního polynomu jsou dány posuvy

u1 ,u2 ). Matice tuhosti odpovídá matici Bernoulliho nosníku (58)

a její prvky jsou poèítány tøíbodovou Gaussovou integraèní formulí (co znamená, e pøesnì a numerická integrace zde není pouita jako nástroj umìlého "zmìkèení" prvku). Element SHELLAX je koncipován naprosto stejnì jako kvalitní nosníkový element BEAM3D, pøesto vak je jeho chování podstatnì horí. Velmi jednoduchý prvek zaloený na analogii Timoenkova nosníku (zde se pouívá název Mindlinovská skoøepina), uvaující pouze lineární prùbìhy posunutí a nezávislých natoèení prùøezu, dává ve srovnání s elementem SHELLAX a øádovì pøesnìjí výsledky. Poznámka týkající se modelování: Pøi tvorbì modelu rotaènì symetrické skoøepiny definujeme meridián jako KØIVKU, kterou pak pokrýváme elementy SHELLAX; velkou pozornost je pøitom tøeba vìnovat jemnosti dìlení a vyuívat informace, které nám poskytuje analytické øeení: Tou základní informací je rozsah poruch membránové napjatosti skrytý v koeficientu β (72). Délka elementu toti musí být v kadém pøípadì mení ne délka vlny napìtí

L perioda =

2π [m] , jinak by podstatné detaily prùbìhù napìtí mohly uniknout (je tøeba si β

uvìdomit, e tvarové funkce /polynomy/ nemohou dost dobøe aproximovat napø. sinový prùbìh vln). Pro popis rovinné deformace nebo rovinné napjatosti (membrány; napøíklad tenká stìna zatíená silami, které pùsobí v její rovinì), ale té pro modelování rotaènì symetrických tìles (høídele, tlustostìnné nádoby kdy nelze pouít nìjaký jednoduchý pøedpoklad o tvaru deformovaného prùøezu jako u skoøepin) se


47

pouívají rovinné trojúhelníkové nebo ètyøúhelníkové prvky (v systému COSMOS nazývané TRIANG a PLANE2D). Matice tuhosti se získá principiálnì stejným postupem jako u jednorozmìrných prvkù, tj. z vyjádøení deformaèní energie jednoho elementu. Napø. pro rovinnou deformaci ( ε zz = 0 ) platí vzhledem k (41) vztah (u,v oznaèují posuvy ve smìrech x,y):

E i (u,v) =

1 1 r tr ∫ ∫ ∫( ε xx σ xx + ε yy σ yy + γ xy σ xy )dxdydz = ∫ ∫ ∫ ε T Dε dxdydz, 2 2

kde

 ∂u     ∂x        1- µ  µ  ε xx    t ∂v E r    µ 1- µ  D= ε =  ε yy =  ∂y  (1+ µ )(1- 2 µ )         γ xy 0  0     ∂u + ∂v    ∂y ∂x 

  0  0 .   1- 2 µ   2 

Poznámka: Aplikaènì dùleitìjí pøípad rotaènì symetrických tìles se lií jen vektorem deformace v

u (u oznaèuje posunutí v radiálním smìru). Matice elastických rt vznikne stejnì jako matice D(3x3) (75) z obecné formulace Hookova zákona (41) vyputìním

nìm je navíc tangenciální sloka konstant

t D(4x4)

ε ϕϕ =

tìch øádkù a sloupcù, kterým odpovídají nulové deformace. Pøípad rovinné napjatosti se od varianty rovinné t deformace lií jen mírnì pozmìnìnou maticí elastických konstant D(3x3) , viz. napø. Zienkiewicz [10]. Následující postup je pro vechny uvedené modifikace vpodstatì totoný.

Prvním krokem musí být vyjádøení vektoru deformace vektorem uzlových parametrù (posunutí). Uvaujme ètyøúhelníkový prvek se ètyømi uzlovými body; kadému z nich pøísluí dva stupnì volnosti (u,v),

r q = ( u1 ,v1 ,...,u4 ,v 4 ) . Posunutí v libovolném místì uvnitø elementu (a na jeho základì deformaci) vyjádøíme prostøednictvím ètyø tvarových funkcí N i (x, y) take element bude mít celkem 8 uzlových parametrù (kadá odpovídá jednomu uzlovému bodu):


   ∂  ∂x    ε xx      = ε  0 yy        γ xy   ∂   ∂y

    0  ∂  ∂x       u ∂  .  =  0  v  ∂y       ∂ ∂   ∂x   ∂y

48

  0     N1 ∂  .  0 ∂y    ∂  ∂y 

0 N1

N 2 ... N 4 0 ...

0

 u1    v1   0  . . N 4  .  .    v 4

t r - - ∆ - - → r - - - - -N - - - - - - → t r ------------ B-------- →

Matice tuhosti je pro prvek s konstantní tloukou stìny dána pouze dvojnými (a ne trojnými) integrály

t t t tt K = s ∫ ∫ BT DBdxdy. kde B(x, y) je matice prvních derivací tvarových funkcí N i (x, y) , viz. (76). Jedinou nezodpovìzenou otázkou tak zùstává návrh vhodných tvarových funkcí

N i (x, y) a zpùsob

integrování (77). Takøka výluènì se vyuívá transformace, která pøevádí ètverec z roviny pomocných bezrozmìrných souøadnic

ξ ,η na uvaovaný ètyøúhelník (jeho strany mohou být event. zakøivené) v rovinì

x,y. Výhodou je pøedevím to, e obecná integraèní oblast v rovinì x,y je nahrazena integrací pøes ètverec v

ξ ,η . Tvarové funkce splòující poadované podmínky v uzlových bodech (tj. e funkce N i má být rovna jedné v uzlu èíslo i a ve vech ostatních nulová) se v rovinì ξ ,η sestrojí jednodue; napø. pro ètyøuzlový ètverec s délkou strany 2 jako bilineární polynomy (polynomy lineární vzhledem k ξ i η ): 1 N i ( ξ ,η ) = (1+ ξ ξ i ).(1+ η ηi ) 4 kde ξ i ,ηi jsou souøadnice i-tého vrcholu (vzhledem k volbì délky

rovinì

strany 2 jsou tyto souøadnice rovny buï 1 nebo -1).

Transformace ètverce na obecný ètyøúhelník se provádí právì prostøednictvím tvarových funkcí (78) 4

4

i=1

i=1

x = ∑ xi N i ( ξ ,η ), y = ∑ yi N i ( ξ ,η ). Z definice (79) je patrné, e vrcholy ètverce

ξ i ,ηi se správnì transformují na vrcholy ètyøúhelníku xi , yi .

Koneèné prvky, které pouívají tyté tvarové funkce pro geometrickou transformaci (79) i pro aproximaci posunutí (76) se nazývají izoparametrické.


49

Dalí postup je ji víceménì rutina; první derivace tvarových funkcí dle x,y (vystupující v matici vyjádøí pomocí derivací vzhledem k

t B ) se

ξ ,η (pouitím pravidla o derivování sloené funkce):  ∂ Ni  ∂x ∂y   ∂ N i        ∂ξ   ∂ξ ∂ξ   ∂x    =   . ,  ∂ Ni  ∂x ∂y   ∂ N i        ∂η   ∂η ∂η   ∂y  t --- J ---

t J je Jacobiho transformaèní matice. Její determinant se nazývá Jacobián, má rozmìr m2 a udává pomìr velikosti plochy elementárního kosodélníka v rovinì x,y k ploe odpovídajícího ètverce v rovinì ξ ,η . Pro transformaci plochy tedy platí dx.dy=det(J) dξ .dη (k ovìøení tohoto vztahu staèí pouèka ze støedokolské

kde

matematiky: plocha kosodélníka je rovna vektorovému souèinu jeho pøilehlých stran). Po dosazení transformaèních vztahù pro derivace tvarových funkcí a plochy do integrálù, kterými jsou vyjádøeny prvky matice tuhosti (76) shledáme, e integrandem ji nejsou polynomy, ale racionální funkce. Numerická Gaussova integrace je tedy nezbytná a pro ploné integrály se provádí takto:

t K = s

1 1

∫∫ 1 1

tt tT B ( ξ ,η )DB( ξ ,η ) det(J)dξdη = s

k

k

∑∑ i=1 j=1

tt t gi g j BT ( ζ i ,ζ j )DB( ζ i ,ζ j ) det(J( ζ i ,ζ j )),

co je pøímoèaré rozíøení Gaussovy formule (63) pro jednoduché integrály ( gi ,ζ i jsou tyté váhové koeficienty a souøadnice integraèních uzlù). k je poèet integraèních uzlù a dá se u vìtiny komerèních programù volit (program COSMOS nabízí variantu k=2 nebo k=3). Tím je moné do jisté míry ovlivnit chování prvku zhruba stejnì jako tomu bylo u prvkù nosníkových: Mení poèet integraèních uzlù sniuje tuhost prvku a napomáhá k odstranìní "shear locking" efektu (neúmìrnému zvýení smykové tuhosti, které se projeví napø. pøi modelování konzoly pásem obdélníkových elementù PLANE2D a pøi ohybovém namáhání). Pro zlepení vlastností membránových prvkù se nìkdy pouívá tzv. selektivní redukovaná integrace spoèívající v tom, e èleny matice

t K , které vyjadøují smykovou tuhost se integrují pouze jednobodovou (k=1), zatímco ostatní dvou

èi tøíbodobou formulí. Pro volbu optimálního k platí jetì jedna smìrnice: Pokud je geometrie ètyøúhelníkového elementu znaènì odchylná od obdélníkové, projeví se to výraznou závislostí Jakobiánu na souøadnicích

ξ ,η , a

pak je vhodné pouít vyí øád integrace. Vìtinou se pouívají tyto hodnoty k: - obdélníkový ètyøuzlový prvek 2 x 2 (k=2) - silnì deformovaný ètyøuzlový prvek 3 x 3 (k=3) - obdélníkový osmiuzlový prvek (uzly ve vrcholech a støedech stran) 2 x 2 a 3 x 3 - osmiuzlový prvek se zakøivenými stranami 3 x 3 a 4 x 4. Pozn.: Vyí øád Gaussovy integrace pouitý u osmiuzlových prvkù odpovídá komplikovanìjímu integrandu (tvarové funkce nejsou v tomto pøípadì bilineární, nýbr kvadratické polynomy). U mebránových prvkù existují jetì dalí zpùsoby sniování jejich tuhosti. Napø. COSMOS nabízí pod


50

tajemnou zkratkou QM6 rozíøení základního ètyøuzlového prvku PLANE2D o dalí dvì tvarové funkce (celkem tedy 6 tvarových funkcí- 4 bilineární polynomy (78) a 2 kvadratické polynomy promìnných

ξ ,η ).

Tyto dodatkové tvarové funkce jsou rovny nule ve vech ètyøech uzlech, ale u ne podél celé strany (parabolický prùbìh tvarové funkce podél jedné strany elementu toti není jednoznaènì urèen pouze dvìma uzlovými hodnotami): 2

N 5 ( ξ ,η ) = 1- ξ ,

N 6 ( ξ ,η ) = 1- η . 2

Dùsledkem je skuteènost, e posuvy u,v u nemusí být spojité na hranicích elementù a v deformované konstrukci tedy vzniknou mezi elementy "kvíry". Z èistì teoretického hlediska je to nepøípustné (a elementy tohoto typu se nazývají nekompatibilní),

ale

v

praxi

èasto

uiteèné

("rozumnost"

chování

nekompatibilních prvkù se vak musí vdy ovìøovat, zpravidla pouitím tzv. "patch" testu [1], který spoèívá v modelování pøípadù, kdy je deformace v celém tìlese konstantní; i nekompatibilní prvky musí v tomto pøípadì dávat pøesné výsledky, jinak jsou zatraceny). Shrnutí: Vechny prvky konstruované na základì Lagrangeova principu jsou tuí ne by mìly být. Zmìkèení lze dosáhnout pouitím komplikovanìjích tvarových funkcí (tj. pouitím polynomù vyího stupnì), protoe tím se rozíøí kála moných tvarù deformovaného prvku (pøesto vak bude prvek vdy o nìco tuí ne "skuteèné" tìleso). Vytvoøit modelovou konstrukci, která by byla mìkèí ne skuteèná, lze jen tak, e poruíme nìkterá pravidla hry, napø. tím, e budeme nepøesnì integrovat (redukovaná Gaussova integrace) nebo obìtujeme poadavek spojitosti posunutí (nekompatibilní prvky).

Nesymetrické zatíení rotaènì symetrických tìles Elementy SHELLAX (skoøepina) nebo PLANE2D umoòují modelovat rotaènì symetrická tìlesa nejen se symetrickým zatíením, nýbr i se zatíením, které se po obvodu mìní (rotaènì symetrická kolona ofukovaná vìtrem, vodorovná válcová nádoba zaplnìná èásteènì kapalinou, høídel na nìj pùsobí kroutící moment a hydrodynamický tlak v kluzných loiscích apod.). V tìchto pøípadech výe uvedené vztahy pro matice tuhosti neplatí, protoe je tøeba uvaovat vechny tøi sloky posunutí (v radiálním, axiálním i tangenciálním smìru) a vech est sloek deformace. To, e i tyto problémy lze chápat jako dvourozmìrné je dáno tím, e závislost vech sloek posunutí a vnìjích sil na souøadnici ϕ se nahradí sinovým nebo kosinovým Fourierovým rozvojem, napø. tlakové zatíení

1 p(r, y,ϕ ) = p0 (r, y)+ 2

∑ (p m

j=1

cj

)

(r, y) cos jϕ+ p sj (r, y) sin jϕ .

Koneènìprvkové øeení se pak rozpadne na øadu "dvourozmìrných" úloh pro jednotlivé èleny Fourierova


rozvoje j (napø. na výpoèet posunutí

pcj (r, y) ).

51

ucj (r, y),v cj (r, y),wcj (r, y) pro j-tou Fourierovu sloku zatíení

Pøi zadávání takovéto úlohy napø. v programu COSMOS je tøeba specifikovat místo jedné èíselné hodnoty tlaku p (obvykle pøíkazem PCR-Pressure on CuRve) vechny potøebné Fourierovy koeficienty rozvoje (83) (pøíkazy AXCOS nebo AXSIN pro kosinové, resp. sinové koeficienty). Známe-li pouze prùbìh

p( ϕ )

mùeme koeficienty rozvoje (83) stanovit na základì ortogonálních vlastností goniometrických funkcí; v programu COSMOS se výpoèet Fourierových koeficientù z tabelárnì zadaného prùbìhu funkce

p( ϕ ) provádí

pøíkazem AXLOAD (v posledních verzích COSMOSu byly tyto pøíkazy vyøazeny).

Skoøepinové prvky Aplikaènì jsou nesmírnì dùleité obecné skoøepinové prvky, které lze chápat jako kombinaci membránových prvkù (respektujících membránová napìtí, napø. PLANE2D) a deskových prvkù, které musí vystihnou chování elementu (relativnì tenké destièky) pøi ohybovém namáhání. Matice tuhosti skoøepinového prvku je tedy souètem matice tuhosti membránového a deskového prvku. Sestrojení kompatibilního deskového prvku, který by byl konformní s Kirchhoffovou teorií tenkých desek je nároènìjí ne pøedchozí úlohy napø. z tohoto dùvodu: Prùhybová plocha w(x,y) musí být funkce nejen spojitá, ale i hladká (musí mít spojité první derivace a to i na rozhraní elementù; jinak by toti skoková zmìna natoèení prùøezu odpovídala nekoneènì velikým vnitøním ohybovým momentùm). Sestrojit tvarové funkce, které by zajiovaly spojitost prvních derivací není jednoduché a napø. v trojúhelníkovém deskovém prvku by uvedené

poadavky

mohly

splnit

( w(x, y) = a1 + a 2 x + a 3 y + a 4 x

2

a

teprve

polynomy

pátého

stupnì

s

21

koeficienty

5

+...+ a 21 y ), které musí být jednoznaènì urèeny 21-ti uzlovými

parametry (zpravidla vechny nulté, první i druhé derivace w(x,y) ve vrcholech /3x(1+2+3)=18 parametrù/ a jetì tøi dodateèné podmínky uprostøed stran /pøedepsané normálové derivace prùhybù/). Tento prvek je nejjednoduí analogií Bernoulliho nosníkového prvku, kde pro zajitìní hladkosti staèily polynomy pouze tøetího stupnì. Aèkoliv je popsaný trojúhelníkový prvek pozoruhodnì pøesný, není v praxi pøíli pouíván, protoe druhé derivace coby uzlové parametry komplikují jeho napojení na jiné typy elementù (typické uzlové parametry jsou toti tøi posunutí a tøi natoèení). Z tohoto hlediska jsou výhodnìjí následující prvky: Mindlinovské desky Jsou analogií Timoenkova nosníku, tj. prùhyb w(x,y) a natoèení deformovaného prùøezu

ψ x (x, y),ψ y (x, y) jsou pokládány za nezávislé velièiny (normála ke støednicové ploe zùstane i po deformaci pøímková, ale ji ne kolmá k deformované ploe). Uplatòuje se tedy vliv smyku, proèe se tyto tzv. Mindlinovské desky pouívají i pro modelování tlustých skoøepin. Opìt se ovem vynoøují problémy s parazitní smykovou tuhostí (shear locking) a øeí se stejnì jako u Timoenkova nosníku, tj. selektivní Gaussovou integrací (pøi bilineární aproximaci prùhybu desky w(x,y) i natoèení

ψ x ,ψ y se pro integraci èlenù vyjadøujících

smykovou tuhost uívá jednobodová Gaussova integrace). DKT (Discrete Kirchhoff Theory) trojúhelníkové desky [1]


52

Jde o trojúhelníkový prvek v jeho uzlech se zadávají stejné parametry jako u prvku pøedchozího (prùhyb a dvì natoèení). Kirchhoffovy pøedpoklady o zachování normály k deformované ploe se uplatòují pouze ve vrcholech a ve støedech stran (odtud pojmenování diskrétní Kirchhoffova teorie), ale pøi malých rozmìrech se prvky chovají jako opravdové Kirchhoffovy desky. Prùbìhy natoèení

ψ x ,ψ y se aproximují

kvadratickými, zatímco prùhyby w(x,y) kubickými polynomy. V systému COSMOS je mono pouít tyto skoøepinové (a tedy i deskové) prvky, Seide [9]: SHELL3T (tlusté skoøepiny); trojúhelníkový prvek dle Mindlinovy teorie s lineární aproximací prùhybu i pootoèení. SHELL4T (tlusté skoøepiny); ètyøúhelníkový prvek ve tøech variantách QUAD2 - spojení dvou trojúhelníkových prvkù SHELL3T QUAD4 - spojení ètyø prvkù SHELL3T a eliminace parametrù, které odpovídají spoleènému uzlu v prùseèíku uhlopøíèek ètyøúhelníku (tzv. statická kondenzace parametrù) QUAD - bilineární aproximace w(x,y),

ψ x (x, y),ψ y (x, y)

SHELL3 (tenká skoøepina); a na malé odchylky jde o výe popsaný prvek DKT (Kirchhoff) SHELL4

(tenká skoøepina); ètyøúhelníkový prvek s variantami QUAD2,QUAD4, vzniklými skládáním

trojúhelníkových prvkù SHELL3. SHELL3L (vrstvená skoøepina, sendvièe, lamináty); trojúhelníkový prvek se stejnými stupni volnosti jako SHELL3, ale sloený a z 50-ti vrstev, které mají rozdílné materiálové (obecnì neizotropní) vlastnosti. Pro kadou vrstvu zvlá se zadávají tøi parametry: - tlouka vrstvy - èíslo skupiny materiálových parametrù (moduly prunosti ve smìrech x,y /EX,EY/, Poissonovy konstanty /NUXY,NUYZ,NUXZ/, souèinitelé teplotní roztanosti /ALPX,ALPY/ a dále té meze pevnosti v tahu a tlaku /SIGXT, SIGYT, SIGXC, SIGYC/ a meze pevnosti ve smyku, popisující soudrnost sousedních vrstev /SIGXYT,SIGXYC/1)

β , udávající smìry ortotropie v rovinì prvku (napø. EX je modul prunosti vztaený k ose, která je pootoèena o úhel β vzhledem k ose x souøadného systému elementu2).

- materiálový úhel

3.3

Tøídimenzionální prvky - TÌLESA

1

Zadávané meze pevnosti nemaji vliv na výpoèet deformací a napìtí nýbr slouí k vyhodnocení stavu napjatosti podle Tsai Wuovy

nebo Hillovy teorie, [33]. 2

Souøadný systém elementù je definován dle jistých konvencí, osa x je napø. spojnice uzlù 1 a 2.


53

Pokud nelze vechny èásti modelované konstrukce nahradit jedno nebo dvourozmìrnými prvky nezbývá ne se uchýlit k prvkùm trojrozmìrým, zpravidla ètyøstìnùm nebo estistìnùm (nazývaným krátce cihla /brick/). Po teoretické i algoritmické stránce jsou naprostou obdobou izoparametrických membránových prvkù PLANE2D: Obecnì tvarovaná "cihla" se vhodnou geometrickou transformací pøevede na jednotkovou krychli v prostoru

ξ ,η ,ζ . Tato transformace je zprostøedkována stejnými tvarovými funkcemi jako

aproximace posunutí u,v,w; pro 8-mi uzlový prvek (s pøímými hranami) jsou tvarové funkce trilineární polynomy promìnných

ξ ,η ,ζ

1 N i ( ξ ,η ,ζ ) = (1+ ξ ξ i )(1+ η ηi )(1+ ζ ζ i ), 8 zatímco pro 20-ti uzlový prvek (uzly ve vech vrcholech a støedech stran) jsou to neúplné kvadratické polynomy; napø. tvarové funkce odpovídající vrcholùm jsou dány pøedpisem

1 N i ( ξ ,η ,ζ ) = (1+ ξ ξ i )(1+ η ηi )(1+ ζ ζ i )( ξ ξ i + η ηi + ζ ζ i - 2), 8 jeho korektnost lze snadno ovìøit. Výe uvedený prvek (v systému COSMOS nazývaný SOLID /8 nebo 20 uzlù/) má stejnì jako PLANE2D urèitou nevýhodu v tom, e uzlové parametry jsou pouze posuvy a ne pootoèení- na tìleso vytvoøené prvky SOLID tedy nelze pøipojit napø. skoøepinovou desku SHELL, která by byla pevnì "pøivaøena"; spojení bude mít vdy charakter kloubu. Cihlové prvky s rotaèními stupni volnosti jsou zatím jen ve stadiu vývoje. V programu COSMOS se od verze 1.65 objevuje prvek TETRA4R - ètyøstìn s posuvnými i rotaèními stupni stupni volnosti ve ètyøech uzlových bodech, popisovaný v èlánku Pawlak 1991.

4.

Nelineární statika Linearita vztahu mezi posunutím a pùsobícími vnìjími silami mùe být naruena z tìchto dùvodù

-

dochází k velkým posunutím a natoèením prvkù konstrukce (ale deformace zùstávají malé /desky pøi vìtích prùhybech, málo klenutá víka nádrí, nìkteré pøíhradové konstrukce/)

-

projevují se velké deformace (nelze pouít tenzor malých deformací)

-

vstoupí-li do hry fyzikální-materiálová nelinearita (hyperelasticita /pry/, plasticita, teèení materiálu) tj. kdy neplatí lineární vztah mezi napìtím a deformací.

4.1

Velké deformace a posunutí


54

Jako ilustraci pøípadu, kdy se výraznì uplatòuje

vliv

velkých

posunutí

na

tuhost

konstrukce uveïme døíve ji zmínìný pøípad vzpìradla (obr.4), co je soustava dvou táhel s jediným stupnìm volnosti u (posunutí styèníku dvou táhel), kterému odpovídá jedna vnìjí síla F. Cílem je odvodit závislost této síly na posunutí u a posléze i "matici" tuhosti vzpìradla (která má v tomto pøípadì jen jediný øádek a sloupec). Protoe systém je staticky urèitý, lze vztah mezi vnìjí silou F a vnitøní silou R pùsobící v táhlech vyjádøit pøímo (viz. obr.4)

F = 2R sinϕ = 2R

H -u . L

Zamìøme se nyní na výpoèet deformace, která je dána zmìnou délky táhla z pùvodní hodnoty L na na délku l, pøi posunutí styèníku o vzdálenost u. Vztah mezi tìmito tøemi velièinami lze vyjádøit z Pythagorovy vìty (D,Hoznaèují délky odvìsen v nedeformovaném stavu)

L 2 = D2 + H 2

2

l 2 = D2 + (H - u ) , 2

L2 - l 2 = H 2 - (H - u ) = 2Hu - u2 , odkud plyne následující výraz pro relativního prodlouení (pøi smìru pùsobení síly F dle obr.4 zkrácení) prvku

ε= Pozn.:

L-l 2Hu - u2 Hu 1 u 2 L2 - l 2 _ = = - ( ). L L(L+ l) L(L+ l) L2 2 L Kdybychom poèítali deformaci táhla "klasicky" jako prùmìt posunutí u do osy táhla, dostali bychom

lineární vztah

ε=

uH u sin ϕ = 2 , který odpovídá definici tenzoru malých deformací (a je to jen první èlen L L

obecnìjího vyjádøení deformace (88)). Pøesnì toto jsme dìlali pøi odvozování matice tuhosti táhla v kapitole 2; porovnej s rovnicí (20).

Hookùv zákon (stále platí!) nyní vyjadøuje lineární vztah mezi napìtím (vnitøní silou R) a velkou deformací (88)

 Hu 1 u R = E.S.ε = E.S.  2 - ( L 2 L

 )2  , 

kde S oznaèuje prùøezovou plochu táhla. Dosazením (89) do (86) obdríme hledanou závislost vnìjí síly F na posunutí u:

 Hu 1 u F = 2.E.S.  2 - ( L 2 L

 H -u )2  ,  L

která jak patrno není lineární. Ze závislosti F(u) mùeme stanovit teènou tuhost konstrukce (pomìr malého


55

pøírùstku síly δF k odpovídajícímu posunutí δu ) derivací (90)

Kt =

dF 2ES  H u Hu 1 u 2  = ( 2 - 2 )(H - u) - [ 2 - ( ) ].  du L  L L L 2 L 

Tuto teènou tuhost lze víceménì formálnì rozloit na souèet tøí tuhostí

2ES H 2 4ES Hu 1 u 2 2R ( ) [ 2 - ( ) ]L L L L L 2 L 2R 4ES 2ES H 2 ε , Kσ = , ( ) , Ku = KL = L L L L

Kt =

K L je lineární tuhost (poèítaná "klasickým" zpùsobem na základì teorie malých deformací), K u sníení tuhosti vlivem zmìny geometrie (jak vidno je lineární funkcí deformace ε ) a K σ je zmìna tuhosti daná

kde

existencí poèáteèních napìtí. Pozn.: Pro lepí objasnìní významu K σ zvame co se odehrává v jednom malém zatìovacím kroku (na jeho poèátku existuje ji urèitý prùhyb u a v táhlech pùsobí síly R, které jsou v rovnováze s vnìjím zatíením

h /F je prùmìtem síly 2R do osy vzpìradla/). Pøi malé zmìnì posunutí δu dojde k L natoèení táhla o úhel δϕ a zmìní se tudí i hodnota prùmìtu síly R (pøi této úvaze pøedpokládáme, e se vnitøní F = 2R sinϕ = 2R

síla bìhem zatìovacího kroku nemìní)

δF = 2R sin ϕ - 2R sin( ϕ + δϕ ) = 2

h - δu 2R Rh = - 2R δu = K σ δu . L L L

Shròme: Na rozdíl od teorie malých deformací jsme pøi velkých posuvech museli uvaovat -

nelineární vztah posunutím a deformací (88)

-

rovnováhu vnitøních a vnìjích sil pro deformovanou geometrii (86). V obecném pøípadì (kdy je více ne jeden stupeò volnosti konstrukce) budou

t t t t K t , K L , K u , K σ matice

(teèná matice tuhosti, lineární matice tuhosti /tá, kterou jsme uvaovali v pøedchozí kapitole/, deformaèní matice tuhosti a koneènì matice geometrické tuhosti /té matice poèáteèních napìtí/). Východiskem pro toto obecné øeení je princip virtuálních posunutí (døíve to byl Lagrangeùv princip minima celkové energie). Co vlastnì rozumíme pod pojmem virtuální posunutí

r δq ? Jsou to nekoneènì malá posunutí bodù vyetøované

konstrukce, která vyhovují pøedepsaným okrajovým podmínkám. Napø. u výe uvedeného vzpìradla je virtuální posunutí styèníku táhel prostì diferenciál du = δu . Princip virtuálních posunutí praví, e práce, kterou konají vnitøní síly

r r r σ pøi virtuálním posunutí δq (pøi malé hypotetické deformaci) se rovná práci vnìjích sil F : r

r

r

∫ ( δε ) σdV - ( δq ) T

T

r F = 0.

V

Pozn.: Pøi zkoumání významu rovnice (93) je asi dobré se vrátit k Lagrangeovu principu a ke zpùsobu vyjádøení energie vnitøních (42) a vnìjích sil (43). Rozdíl je v tom, e vztah (93) uvauje jen malé zmìny vnitøní a vnìjí


56

energie pøi virtuálním posunutí. Moná vás napadne otázka, proè se v rovnici (93) neobjeví i variace sil Odpovìï je následující: Ano, i pøi velice malých posunutích

r δσ ?

r δq dojde k malé zmìnì rozloení vnitøních napìtí

r r r δσ , leè zmìna energie s nimi spojená δ ε T δσ je øádovì mení ne zmìna uvaovaná ve vztahu (93), tj. r r δ ε T σ . To je také dùvod proè pod pojmem virtuální posunutí rozumíme nejen libovolná myslitelná, ale

pøedevím dostateènì malá posunutí.

r r δε jsou jednoznaèným zpùsobem spojeny s virtuálními posuvy δq a tato r r závislost je lineární (i kdy vztah mezi ε a q lineární není je zachována linearita mezi nekoneènì malými Virtuální deformace

pøírùstky a v tom je vtip vech pøírùstkových metod); zapime tento vztah v maticové notaci

t r r r t r t r r r r δε = B(x,q)δq = ( B L (x)+ B N (x,q)).δq, t B je matice svazující vektor deformací a vektor posunutí (má stejný význam jako døíve, viz. napø. (47), zde t je ovem omezena jen na vztah mezi malými pøírùstky; "konstantní" èást matice B L , která nezávisí na

kde

posunutích je dokonce totoná s døíve pouívaným vyjádøením deformaèní matice, napø. (76) pro rovinné prvky PLANE2D). Dosadíme-li vztah pro virtuální deformaci do

r t r r T ( δq ) ( ∫ BT σdV - F) = 0

(93) získáme de facto

V

_

∫

V

r tT r B σdV - F = 0.

integrální rovnici rovnáhy

Tato identita musí platit i pøi libovolné variaci (virtuální zmìnì) posunutí; poznamenejme, e pro variaèní poèet platí zhruba stejná pravidla jako pro diferenciální poèet a e tøeba pro variaci souèinu platí pouèka

r δ (uv) = δuv + uδv , kterou pouijeme pøi variaci rovnice rovnováhy (95), v ní se pøi variaci posuvù δq t r mìní matice B i napìtí σ : t T r r tT r ∫ ( δB ) σdV + ∫ B δσdV = δF. V

V

t r r δσ i variaci matice B variací posunutí δq , k èemu t pouijeme Hookùv zákon, zapsaný ve stejném tvaru jako døíve (matice elastických konstant D je napø. pro

V rovnici (96) se nyní pokusíme nahradit variaci napìtí trojrozmìrnou napjatost definována vztahem (41))

tt r r t r δσ = Dδε = DBδq,

∫

V

tT r tT t t tT tT t t t t t r r r B δσdV = ∫ B DBdV.δq = ∫ ( B L + B N )D( B L + B N )dV.δq = ( K L + K u ).δq, V

t tT t t K L = ∫ B L D B L dV, V

V

t tT t t tT t t tT t t K u = ∫ ( B N D B L + B L D B N + B N D B N )dV V


kde matice tuhosti

57

t t K L je totoná s døíve analyzovanými lineárními pøípady (na rozdíl od matice K u nezávisí

r q ). t t t Variace matice B je ve skuteènosti jen variací její "nekonstantní" sloky B N a dá se nalézt taková matice K σ , na posunutích

která splòuje následující identitu

∫

V

tT r tT r t r δ B σdV = ∫ δ B N σdV = K σ δq. V

Vyjádøení (98) naznaèuje, e matice geometrické tuhosti

t r K σ je lineární funkcí napìtí σ a vùbec nezávisí na

materiálových parametrech (modulech prunosti apod.). Shrnutím vztahù (97) a (98) dostáváme avizovanou závislost mezi malou zmìnou posunutí a malou zmìnou sil (zatíení) v tomto maticovém tvaru

r t t t r ( K L + K u + K σ )δq = δF.

Poznámka pro pokroèilé: Popis struktury matice teèné tuhosti (99) je dosti veobecný; vùbec z nìho neplyne napø. návod na výpoèet matice geometrické tuhosti

t K σ . Dùvod je ten, e pouívání maticové notace

("vektorové" vyjadøování tenzorù napìtí a deformace) je zvlátì u nelineárních problémù tìkopádné, nepøehledné a pracné. V klasických uèebnicích (napø. Zienkiewicz [10]) zabírá odvozování matic tuhosti zaloené na maticové symbolice desítky stránek. Tenzorová formulace je mnohem prùhlednìjí; vychází z následujícího slokového vyjádøení tenzoru velkých deformací (Greenùv tenzor, deformace "vzpìradla" (88) je speciálním pøípadem)

1 ∂ ui ∂ u j ∂ um ∂ um + + ), ε ij = ( 2 ∂ x j ∂ xi ∂ xi ∂ x j kde

r r r ui (x) jsou sloky vektoru posunutí. Aproximujeme-li ui (x) tvarovými funkcemi N ik (x) r r ui (x) = ∑ N ik (x) q k = N ik q k , k

získáme vztah mezi tenzorem velkých deformací a vektorem zobecnìných posunutí uzlových bodù

r q

1 ∂ N mk ∂ N ml 1 ∂ N ik ∂ N jk q q . + )qk + ε ij = ( ∂ xi 2 ∂ xi ∂ x j k l 2 ∂xj Variace tenzoru deformace je vlastnì diferencování (F3) vzhledem k posunutím

r q

1 ∂ N ik ∂ N jk ∂ N mk ∂ N ml q δ qk δ ε ij = ( + ) δ qk + 2 ∂xj ∂ xi ∂ xi ∂ x j l - - - BL - - - - - BN - - - . t t t r Vztah (F4) je ji zcela konkrétní pøedpis výpoètu B L , B N (a je zøejmé, e B N je lineární funkcí posunutí q , t t r take variace δ B N je lineární funkcí variace posuvù δq ). Dosazením variace δ B N do vztahu (98)


58

získáme obecné vyjádøení prvkù matice geometrické tuhosti v závislosti na tvarových funkcích a napìtí

t σ.

  ∂ N mk ∂ N ml  t r  ∂ N mk ∂ N ml  ∫ δ BTN σdV  = ∫ δ q l σ ij dV =  ∫ σ ij dV  δ q l = K σ (kl) δ q l V  k V ∂ xi ∂ x j  V ∂ xi ∂ x j  t t Kontrolní otázka: Co jsou matice B L , B N ve výe uvedeném pøípadì vzpìradla (vechny tyto matice budou mít samozøejmì jen jeden prvek)? Rovnice (99) je vlastnì ji urèitý návod na praktické øeení problémù velkých deformací; místo nekoneènì malých variací se uvaují koneèné pøírùstky (buï posunutí nebo zatìovacích sil), opakovanì se poèítá teèná matice tuhosti a øeí pøísluná soustava rovnic. Zpùsobù øeení je vícero, blíe se o nich zmíníme a pozdìji. Jen malá poznámka: Matice geometrické tuhosti se v programu COSMOS mùe pouít i pøi lineární analýze, která se skládá ze dvou fází: v první se øeí soustava rovnic s "klasickou" lineární maticí tuhosti, na základì tohoto øeení se vypoète matice geometrické tuhosti a celé øeení se opakuje (s maticí tuhosti

t t K L + K σ ). Tento postup se provádí v pøípadì, e je nastaven parametr elementu (pøepínaè) In_plane_stiffnes.

V systému COSMOS jsou matice geometrické tuhosti a velkých deformací pøipraveny pro témìø vechny prvky; výjimkou jsou PIPE,TPIPE,ELBOW (potrubní systémy), TRIANG, TETRA (trojúhelníkové membránové prvky a ètyøstìny).

K dispozici je tedy m.j. i prvek táhlo TRUSS2D v nelineárním provedení,

take lze ukázat postup nelineární analýzy vzpìradla (ji jsme ho poèítali dvakrát, ale vdy na základì teorie malých deformací). Definice geometrie, materiálových parametrù, okrajových podmínek, volba typù elementù i generování sítì je naprosto stejné jako v programu uvedeném na str.25. Lií se jen v tìchto ohledech: - je tøeba definovat pouitou metodu øeení - je nutné specifikovat "èasový" prùbìh zatìování nebo posunutí. Uvozovky byly pouity proto, e se nejedná o nestacionární dìj v nìm by se uplatòovaly setrvaèné síly, ale èas hraje roli jen jakéhosi poèitadla zatìovacích krokù. Text vzorového programu je následující:

PT,1,0,0,0, Souøadnice 3 vztaných bodù PT,2,1,0,0, PT,3,.5,.01,0, Souøadnice styèníku (0.5, 0.01) m CRPLINE,1,1,3,2,2, Vytvoøení dvou køivek (táhla) EGROUP,1,TRUSS2D,0,0,0,0,0,1,0, Typ elementu, poslední èíslice 1 avizuje velké deformace MPROP,1,EX,2E11, Modul prunosti RCONST,1,1,1,1,1E-4, M_CR,1,2,1,2,1,1, Generátor sítì (pokrytí køivek 1 a 2) vdy jediným elementem NMERGE,1,4,1,0.0001,0,1,0, DND,1,UZ,0,4,1,, DPT,1,UX,0,2,1,UY,, TIMES,0,100,1, Poèet èasových krokù 100, velikost kroku 1 CURDEF,TIME,1,1,0,0,100,1, Definice funkce è.1, lineární závislost od 0 do 1 ACTSET,TC,1, Aktivujeme funkci è.1 FPT,3,FY,-65,3,1, Zatíení ve styèníku silou Fy=-65 x aktivovaná_funkce(t) NL_RESGRAPH,1,2,2, Ukládání hodnot posuvù UY v uzlu èíslo 2 pro vykreslení grafu UY(t) R_NONLINEAR Sputìní výpoètu pro implicitní metodu a pøesnost øeení

Výsledky jsou uvedeny na následujícím obrázku (prùbìhy posunutí uzlu èíslo 2 /styèníku/ v závislosti na "èase" - zatíení, a lineární prùbìh zatìovací síly):


59

Zobrazené grafy byly získány následujícími pøíkazy: ACTXYPLOT,0,1,TIME,DISP,2,2; XYPLOT,1; závislost w(t) -posunutí styèníku CURPLOT,TIME,1 zatìovací funkce F(t) Teoretické øeení (90) pøedpovídá, e kritická síla F pøi ní by mìlo dojít k pøeskoku ("propadnutí" vzpìradla) je 64.47 [N] a odpovídající hodnota posunutí w=0.004226 [m]; z výsledkù programu COSMOS je patrné, e jetì pøed dosaením tohoto kritického bodu se objevily potíe (divergence iterací v 94-tém zatìovacím kroku, tj. pøi zatíení F=61.1 [N]) a výpoèet byl pøedèasnì pøeruen. Uvedeným postupem tudí nelze sledovat co se s konstrukcí dìje po pøekroèení kritického zatíení (postbuckling analysis); pro studium tìchto dìjù je moné zvolit místo metody zatìovacích krokù metodu pøírùstkù posuvù (jako "èasovou" funkci bychom v tomto pøípadì definovali posunutí styèníku a výsledkem by byl odpovídající prùbìh zatíení).

4.2

Elasticko-plastické chování materiálu Pøekroèí-li v nìkterém místì konstrukce efektivní napìtí mez kluzu

σ K , pøestane platit Hookùv zákon

a po odtíení se konstrukce nevrací do pùvodního stavu. Toto elasticko-plastické chování je typické pro tlakové nádoby (k místní plastizaci dochází v místech poruch membránové napjatosti) a nìkdy je ho tøeba uvaovat i u potrubních systémù (obvykle v souvislosti s vlivem teplotních dilatací). Nejjednoduí analytické metody jsou zaloeny na pøedstavì ideálnì plastického materiálu, který je charakterizován tím, e napìtí zùstává po pøekroèení meze plasticity konstantní a nedochází ke zpevòování materiálu. Dle tohoto modelu by se napøíklad trubka zatíená pøetlakem, kterému odpovídá obvodové napìtí rovné mezi kluzu, okamitì "roztekla" (bylo by dozaeno mezního stavu vzhledem k nepøípustnì velkým deformacím). Pøedpoklad konstantnosti napìtí podstatnì zjednoduuje výpoèet mezního zatíení konstrukce; známým pøíkladem je stanovení mezního ohybového momentu, který zpùsobí úplné zplastizování prùøezu nosníku (vytvoøení plastického kloubu). Napø. pro nosník obdélníkového prùøezu (íøka b výka h) platí

Mp=

σ K bh h bh2 , . =σK 2 2 4

co je 1.5 krát více ne ohybový moment, pøi nìm je pouze v krajních vláknech dosaeno meze kluzu. Poznamenejme, e pomìr "plastického" a "elastického" momentu závisí na geometrii prùøezu nosníku a napø. u I-profilu je pouze cca 1.15. Pøi výpoètu pøíhradových konstrukcí a potrubních sítí je nejvìtím problémem zjistit v kterém místì se plastické klouby vytvoøí; analýzy tohoto druhu jsou pomìrnì sloité a opírají se o tzv. simplexovou optimalizaèní metodu (ejnoha 1990). V pøípadì obecné elasto-plastické napjatosti existují vpodstatì dva postupy øeení: -

Deformaèní teorie plasticity, zaloená na pøedpokladu, e existuje jednoznaèný vztah mezi aktuálním

napìtím

σ a deformací ε (tj. e nezáleí na zpùsobu zatìování - cestì, kterou se soustava do daného stavu

dostala). S touto teorií (Hencky, Nádai) jste byli seznámeni v základním kurzu prunosti a pevnosti (Anonym [32]), kde byla pouita napø. pro výpoèet elasto-plastického stavu silnostìnné trubky. Pro pøípad jednoosé napjatosti mùeme tuto hypotézu vyjádøit vztahem

ε = εe+ ε p = (

1 1 + )σ , E Hs


kde

60

E je Youngùv modul prunosti a H s seèný modul plasticity (který ovem není materiálová konstanta, ale

závisí na okamité hodnotì plastické deformace). Poznamenejme, e vzhledem ke svým omezením se deformaèní teorie plasticity v MKP prakticky nepouívá i kdy její jednoduchost láká (Schneider 1990). -

Pøírùstková teorie plasticity se snaí sledovat celý prùbìh zatìování a formuluje vztahy pouze mezi

malými pøírùstky napìtí a deformace. Napø. pro jednoosou napjatost vztahem

dε = d ε e + d ε p = ( kde

1 1 + )dσ , E Ht

H t je teèný modul plasticity /podíl pøírùstku napìtí k

pøírùstku plastické deformace/, který ji vìtinou lze povaovat za materiálovou konstantu (alespoò pro bilineární materiály, jejich pracovní diagram odpovídá obrázku 20). V pøípadì trojrozmìrné napjatosti je tøeba vechny výe uvedené koncepty zobecnit a pøedevím stanovit kritérium, které by na základì okamité hodnoty tenzoru napìtí a okamitém stavu materiálu (napø. míøe jeho zpevnìní) rozhodlo, zda bylo èi nebylo dosaeno meze plasticity. Protoe toto kritérium nesmí záviset na tom jak byl zvolen souøadný systém sloek napìtí, lze pouít buï tøi hlavní napìtí nebo tøi invarianty tenzoru napìtí . Formálním zápisem takové podmínky je rovnice, popisující v prostoru hlavních napìtí plochu, nazývanou plocha plasticity:

F( I σ ,II σ ,III σ ) = σ K ( κ ), kde

σ K oznaèuje okamitou mez kluzu materiálu, která závisí na pøedchozí deformaci ( κ je skalární parametr

zpevnìní materiálu). Pro rùzné materiály se osvìdèují rùzné plochy plasticity, tj. funkce F. Pro kovy (s výjimkou litiny) se pouívá pøedevím Von Misesova válcová plocha plasticity definovaná toliko prostøednictvím druhého invariantu deviátoru napìtí (deviátoru proto, e vliv izotropního tlaku na plastizaci tvárných materiálù je zanedbatelný)

F( II σ ′ ) = 3 II σ ′ = =

3 3 t t σ ′:σ ′ = σ i ′j σ j ′i = 2 2

1 2 2 2 [( σ 1 - σ 2 ) + ( σ 1 - σ 3 ) + ( σ 2 - σ 3 ) ] = 2 =σK ,

kde funkce plasticity F není nic jiného ne døíve ji definovaná intenzita napìtí (32).


61

Pozn.: Dalí pomìrnì èasto uívané podmínky plasticity jsou - Trescovo kritérium (podobné Misesovu kriteriu, pouze místo válce je estiboký hranol tomuto válci vepsaný) - Druckerovo Pragerovo kriterium a Stassiho kriterium (plocha plasticity je kuel event. rotaèní paraboloid v prostoru napìtí); tato kriteria jsou vhodná pøedevím pro stavaøe (beton, ale té napø. litina). Výe uvedené modelové pøedstavy pøedpokládají, e pokud pøi plastizaci materiálu dochází ke zpevnìní (zvìtení meze kluzu), je toto zpevnìní izotropní (to znamená, e se napø. zvìtí polomìr Misesovy válcové plochy plasticity, ale její osa se nezmìní a stále prochází poèátkem souøadného systému

σ 1 ,σ 2 ,σ 3 ).

Nìkteré materiály (i oceli) vak vykazují spíe posunutí ne "nafouknutí" plochy plasticity, pøièem smìr posunutí se obvykle odvozuje od smìru výslednice tøí hlavních napìtí. Tento jev se nazývá Bauschingerùv efekt a pokud jde pouze o posunutí plochy plasticity hovoøíme o kinematickém zpevnìní materiálu. Pro pøiblíení: U materiálu tyèe namáhané tahem se zvýí mez kluzu v tahu (na hodnotu danou nejvìtím pøedchozím tahovým napìtím) a o stejnou hodnotu se sníí mez kluzu v tlaku (pøi izotropním zpevnìní by se i tato mez zvýila). Zabývejme se nyní návrhem obecného konstitutivního modelu elastoplasticity, který by vyjadøoval souvislost mezi malým pøírùstkem tenzoru napìtí

t r r dσ = Dep dε ,

r r dσ a tenzoru deformace dε . Ukáeme, e tento vztah lze

vyjádøit formálnì stejnì jako Hookùv zákon (40-41), tj. kde

t t Dep je matice (6 x 6) analogická matici elastických konstant D (41), jenome

závisí nejen na modulech prunosti, ale i na okamité hodnotì tenzoru napìtí. Vztah (105) je klíè k dalímu øeení MKP, protoe na jeho základì lze sestavit teènou matici tuhosti

t Kt =

∫

V

tT t t B Dep BdV,

a pak postupovat úplnì stejnì jako v pøedchozím pøípadì velkých deformací (tj. pøírùstkovou metodou).

Elastoplastický model (105) vychází vdy z urèitých hypotéz; ta základní øíká, e tenzor pøírùstku plastické deformace má stejné hlavní smìry jako tenzor celkových napìtí. Situace v urèitém ji

r

zplastizovaném místì konstrukce x je na zaèátku zatìovacího kroku plnì charakterizována tenzorem napìtí, jeho tøi vzájemnì kolmé hlavní smìry mùeme povaovat za souøadné osy prostoru napìtí (Haighùv Westergaardùv prostor). Aktuální hodnoty hlavních napìtí na ploe plasticity

σ 1 ,σ 2 ,σ 3 jsou slokami vektoru, který se nachází

F( σ 1 ,σ 2 ,σ 3 ) = σ K . Pøírùstek zatíení vyvolá zmìnu tenzorù elastické i plastické

deformace, a protoe pøedpokládáme shodnost hlavních smìrù napìtí a plastické deformace mùeme tenzor pøírùstku plastické deformace plnì charakterizovat jistým vektorem Haighova Westergaardova prostoru. Poznámka: Nezarazilo vás tvrzení, e symetrický tenzor urèený 6-ti nezávislými èísly lze nahradit vektorem, tj. pouze 3-mi hodnotami? Vysvìtlení: dalí 3 informace jsou ukryty ve smìrech souøadných os HW prostoru (který je sdílen jak tenzorem napìtí tak tenzorem plastické deformace; je-li znám tenzor celkových napìtí jsou tím dle hypotézy è. 1 automaticky dány tøi hlavní smìry tenzoru pøírùstku plastických deformací). A nyní pøijde druhá základní hypotéza urèující smìr vektoru pøírùstku plastické deformace jako normály k ploe plasticity (vektor pøírùstku plastické deformace má dle této hypotézy smìr gradientu plochy plasticity v HW prostoru, co jsou dalí dvì podmínky pro tenzor pøírùstku plastické deformace):

∂F r r dλ , d ε p = ∆Fdλ = ∂σ


62

kde dλ udává zatím neznámou absolutní hodnotu pøírùstku plastické deformace. Pozn.: Není to toté co pøedchozí hypotéza; uvìdomme si, e i kdy by tenzory pøírùstku plastické deformace a celkových napìtí mìly

t

t

stejné hlavní smìry, neznamená to, e by si byly podobné ( d ε p necσ ). Rovnice (107) je vyjádøením asociativního "zákona", pøièem slovo asociativní vyjadøuje skuteènost, e rovnice (107) spojuje dvì principiálnì nezávislé vìci: Podmínku urèující kdy dojde k plastizaci (funkce plasticity F) a vlastní konstitutivní rovnici. Volba kritéria plasticity F tedy dle hypotézy (107) plnì urèuje chování materiálu v plastickém stavu. Velikost pøírùstku deformace (skalární hodnoty dλ ) vyplyne z podmínky, e i na konci uvaovaného zatìovacího kroku se výsledné napìtí nachází na ploe plasticity, co lze vyjádøit poadavkem aby úplný diferenciál (103) byl nulový

dσK ∂F T r dκ = H t dλ . ( r ) dσ = dκ ∂σ Komentáø k rovnici (108): uvìdomte si, e maticový souèin øádkového vektoru (transpozice gradientu F) se sloupcovým vektorem (pøírùstku napìtí) je skalár (matice 1 x 1). Dále: Kdyby materiál nevykazoval zpevnìní,

H t (o ní pozdìji dokáeme, e má význam teèného modulu plasticity) rovna nule (a vektor pøírùstku napìtí by leel na nemìnící se ploe plasticity F = σ K ). byla by novì zavedená konstanta

Nyní vyjádøíme tenzor pøírùstku celkové deformace jako souèet elastické èásti (kde platí Hookùv zákon) a plastické sloky (pouijeme asociativní zákon):

t -1 r ∂F r r r dε = d ε e + d ε p = D dσ + r dλ . ∂σ Z rovnic (108) a (109) mùeme vypoèíst hodnotu dλ tak, e pùvodnì vektorovou rovnici (109) pøevedeme na skalární vynásobením øádkovým vektorem

t T ( ∆F ) D

∂F T t r ∂F T r ∂F T t ∂F ∂F T t ∂F ( r ) Ddε = ( r ) dσ + ( r ) D r dλ = H t dλ + ( r ) D r dλ ∂σ ∂σ ∂σ ∂σ ∂σ ∂σ ∂F T t r ( r ) Ddε ∂σ dλ = ∂F T t ∂F Ht + ( r ) D r ∂σ ∂σ a dosazením do rovnice (109) získáme finální vztah, který platí pro libovolnou funkci plasticity F:

t ∂F ∂F T t   r( r ) D  r D  t t r r t r t ∂F ∂ ∂σ σ  dε = Dep dε . dσ = Ddε - D r dλ =  D t ∂F T ∂F ∂σ  H t + ( r ) D r   ∂σ ∂σ

Zbývá jetì objasnit dùvod, proè byla konstanta

H t v rovnici (108) oznaèena za teèný modul plasticity


63

(a povaována za známý materiálový parametr). Není to opìt nic jiného ne dalí hypotéza (Hillova teorie), která souvisí s interpretací zpevnìní materiálu κ ; dle této hypotézy je pøírùstek zpevnìní pøímo roven hustotì dissipované energie pøi zmìnì plastické deformace

r r r ∂F dκ = σ T d ε p = dλ σ T r = dλF = dλ σ K . ∂σ Posloupnost identit (112) vypadá jednodue, ale není samozøejmá; napø. to, e skalární souèin vektoru napìtí a gradientu funkce F v HW prostoru je pøímo roven funkci plasticity (a tedy mezi kluzu) plyne z Eulerova teorému, který øíká, e pro vechny homogenní a lineární funkce F(x) platí

x

dF =F. dx

Tuté práci dκ vykonanou pøi plastické deformaci lze vyjádøit i souèinem skalární intenzity napìtí a efektivní plastické deformace (viz. diskuse týkající se významu pojmu intenzity napìtí na str. 26)

dκ = σ M d ε p = σ K d ε p Porovnáním (113) a (112) vyplyne, e

d ε p = dλ a z definice (108) vysvitne i význam velièiny H t

Ht =

d σ K dκ d σ K d σ K . = = dεp dλ dκ dλ

Konstitutivní rovnice (111) je zatím stále jetì obecná; konkrétní podoby nabude volbou funkce plasticity F. Pouitím Misesovy podmínky ( F = σ M ) v rovnici (111) bychom získali následující vyjádøení matice elastoplasticity (Prandtlùv Reussùv model)

t t Dep = D -

9 G2 r rT 2 σ ′σ ′ , (H + 3G) σ M

r σ ′ je deviátor tenzoru napìtí, σ M intenzita napìtí, G modul prunosti ve smyku a H modul plasticity3 dσ dσ (smìrnice pracovního diagramu materiálu v plastické oblasti ), který souvisí s teèným modulem ( ) dle dεp dε kde

vztahu plynoucího z obr. 20

dε = (

1 1 1 )dσ = dσ , + H E Ht

Ht =

H . H 1E

Pozn.: Úplné odvození Prandtlova Reussova modelu je pomìrnì zdlouhavé a nebudeme ho uvádìt. Nicménì jeho nejdùleitìjí èást, vyjádøení gradientu funkce plasticity (a tedy smìru plastické deformace) je pomìrnì uiteèné cvièení operací s tenzory. Ukate, e výsledek této operace (derivováním (104)) je

3

r ∂F 3 σ ′ r= , ∂σ 2 σ M

U uhlíkových ocelí je modul plasticity H cca desetkrát mení ne modul prunosti E, zatímco u ocelí austenitických a dvacetkrát

(nerezové materiály vykazují mení zpevnìní).


64

co znaèí, e pøírùstek tenzoru plastické deformace je prostì skalární násobek tenzoru deviátoru napìtí! V systému COSMOS se pouívá výe uvedený postup a von Misesovo nebo Druckerovo Pragerovo kriterium plasticity s výjimkou skoøepin (SHELL3T a SHELL4T), kde se dosaení plastického stavu hodnotí jen na základì skoøepinových sil a momentù (Iljuinovo kritérium). Volitelné varianty se oznaèují kryptogramy: VMI-von Mises Isotropic hardening, VMK-von Mises Kinematic hardening, DP-Drucker Prager, a zadávají se jako parametry elementù typu TRUSS2D, TRUSS3D, BEAM2D, BEAM3D, PLANE2D, SOLID, SHELL3T, SHELL4T pøíkazem EGROUP. Pro vlastní prùbìh øeení je tøeba specifikovat stejnì jako u metody velkých deformací zatìovací køivku a numerickou metodu øeení (poèty iterací v kadém zatìovacím kroku). Postup zadání vysvìtlíme na pøíkladu elastoplastického zatìování nosníku (vetknutý nosník, zatíený na volném konci postupnì se zvyující silou). Pro øeení této jednoduché úlohy lze volit rùzné typy elementù; napø. nosníkový prvek BEAM2D, skoøepinový prvek SHELL4T nebo membránový prvek PLANE2D (ten pouijeme). Jako model plasticity zvolíme VMI (von Mises izotropní zpevnìní) a bilineární prùbìh pracovního diagramu materiálu (elastický modul prunosti kluzu

E = 2,1.1011 [Pa] , plastický modul H = 1010 [Pa] a mez

σ K = 200[MPa] ). Prùbìh zatíení (síly F) zvolíme tak aby modeloval dìj, nazývaný pøizpùsobení

konstrukce (Jao [6]). V èem spoèívá? Pøedstavme si periodicky tlakovanou nádobu (se stejnými prùbìhy tlaku v kadém cyklu). Pokud by v ádném místì nádoby nebylo dosaeno meze plasticity byl by prùbìh napìtí v kadém zatìovacím cyklu stejný. V opaèném pøípadì to neplatí: Dojde-li bìhem prvního cyklu k plastické deformaci, projeví se to zpevnìním materiálu (izotropním nebo kinematickým zvýením jeho meze kluzu), ale souèasnì i vznikem zbytkových deformací a napìtí v odtlakované nádobì. Nový zatìovací cyklus pak ji vlastnì pùsobí na zmìnìnou nádobu a výsledkem tého zatíení je jiná distribuce napìtí. Jestlie v tìchto následujících cyklem nedojde k pøekroèení meze kluzu znamená to, e se konstrukce pøizpùsobila a nehrozí nebezpeèí kumulace plastických pokození. Souèinitel pøizpùsobení (shakedown factor) udává pomìr zatíení na mezi elasticity pøizpùsobené a pùvodní konstrukce. Poznamenejme, e u ideálních elastickoplastických materiálù bez zpevnìní mùe být souèinitel pøizpùsobení maximálnì 2 (co odpovídá pøípadu, kdy zbytková napìtí jsou stejnì velká jako napìtí, pøi nich bylo dosaeno meze kluzu, ale mají opaèný smìr). Základní metody výpoètu souèinitele pøizpùsobení vycházejí z Melanova nebo Koiterova teorému (viz. Jao [6]). Text programu, v nìm je zaifrován prùbìh zatíení zobrazený na obr. 21 je následující

C* C* COSMOS/M Geostar V1.65L C* CRPCORD,1,0,0,0,1,0,0,1,0.005,0,0,0.005,0,0,0,0, definice 4 køivek, stran obdélníka (bokorys nosníku) EGROUP,1,PLANE2D,0,1,0,0,1,0,0, volba typu elementu, 1-oznaèuje VMI (von Mises Isotropic) MPROP,1,EX,2.1E11, Youngùv modul prunosti E MPROP,1,GXY,79E9, modul prunosti ve smyku G MPROP,1,ETAN,1E10, modul plasticity H (ve v [Pa]) MPROP,1,SIGYLD,200E6, mez kluzu RCONST,1,1,1,2,.01,0, tlouka prvku (íøka nosníèku) SF4CR,1,1,2,3,4,0, plocha vymezená 4-mi køivkami M_SF,1,1,1,4,5,6,.5,1, generátor sítì 5 x 6 elementù DCR,4,ALL,0,4,1, vetknutí vech uzlù na køivce 4 TIMES,0,8,1, 8 zatìovacích krokù (pøírùstkù zatíení) CURDEF,TIME,1,1,0,0,1,.9,2,1.0,3,.9,4,0.,5,0,6,.9,7,1.0,0 zatìovací funkce (funkce èíslo 1) PRINT_NDSET,1,37,37, pro zkrácení výstupní sestavy pøedepíeme tisk posuvù pouze v uzlu 37 NL_DEFPLOT,1,2,3,4,5,6,7,8, v krocích 1,2,.. se budou ukládat informace o rozloení napìtí ACTSET,TC,1, aktivace funkce èíslo 1 FND,37,FY,-13,37,1, zatíení v uzlu è.37 silou Fy=-13 x aktivovaná_funkce_è.1 NL_SOLN,1,0, volba metody NR-Newton Raphson


A_NONLINEAR,S,1,1,50,0.001,0,N,0, R_NONLINEAR

65

parametry metody øeení : 50-iterací sputìní výpoètu

Výsledky uloené do souboru obsahují pro kadý zatìovací krok tyto údaje: - posunutí a natoèení konstrukce ve vybraném uzlu 37 (pùsobitì síly); implicitnì by to bylo ve vech 42 uzlech - informace o vech 30-ti elementech; u kadého prvku se uvádí hodnoty vech sloek napìtí ve vech integraèních bodech (pouívá se Gaussova integrace 2 x 2, take pro kadý element se uvádí 4 tenzory napìtí, 4 intenzity napìtí a zpráva o tom, zda je daný bod elementu v plastickém nebo elastickém stavu). Ukázka èásti výstupní sestavy (údaje o elementu èíslo 30, který je v nejvíce exponovaném místì vetknutí nosníku): NONLINEAR STATIC ANALYSIS Time step increment . . . . . . . . . . . . . = 1.0000 Total number of solution steps . . . . . . . . = 8 Total number of equations . . . . . . . . . . = 70 Solution Technique used: Regular Newton-Raphson Stiffness matrix reformation interval = 1 Number of steps between equilibrium iterations = 1 Maximum number of equilibrium iterations allowed . . . = 50 Displacement tolerance for equilibrium iterations . . = 0.10000E-02 Displacements for step 1 (Time = 0.1000E+01) Total number of equilibrium iterations performed = 2 Node_______X1__________X2__________X3__________R1__________R2__________R3 37 0.61878E-03 -.16303E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 Element Stress von Mises Number_IPT_State____Stress-XX___Stress-YY___Stress-XY___Stress-ZZ_____Stress___________Principal_Stresses 30 1 Plastic 2.1987E+08 4.6764E+07 -5.2144E+04 0.0000E+00 2.0061E+08 2.1987E+08 4.6764E+07 0.0000E+00 2 Elastic 1.8555E+08 4.7219E+07 -5.2023E+04 0.0000E+00 1.6702E+08 1.8555E+08 4.7219E+07 0.0000E+00 3 Plastic 2.0564E+08 8.7691E+06 -4.8039E+04 0.0000E+00 2.0140E+08 2.0564E+08 8.7691E+06 0.0000E+00 4 Elastic 1.8563E+08 4.1039E+06 -5.1957E+04 0.0000E+00 1.8361E+08 1.8563E+08 4.1039E+06 0.0000E+00 Displacements for step 6 (Time = 0.6000E+01) Total number of equilibrium iterations performed = 1 Node_______X1__________X2__________X3__________R1__________R2__________R3 37 0.64745E-03 -.17321E+00 0.00000E+00 0.00000E+00 0.00000E+00 0.00000E+00 Element Stress von Mises Number_IPT_State____Stress-XX___Stress-YY___Stress-XY___Stress-ZZ_____Stress___________Principal Stresses 30 1 Elastic 1.9306E+08 4.8566E+07 3.6695E+05 0.0000E+00 1.7394E+08 1.9307E+08 4.8565E+07 0.0000E+00 2 Elastic 1.8937E+08 4.8465E+07 3.8177E+05 0.0000E+00 1.7039E+08 1.8938E+08 4.8464E+07 0.0000E+00 3 Elastic 1.9865E+08 1.6454E+07 -1.9741E+05 0.0000E+00 1.9095E+08 1.9865E+08 1.6453E+07 0.0000E+00 4 Elastic 1.9632E+08 8.1838E+06 -1.8021E+05 0.0000E+00 1.9236E+08 1.9632E+08 8.1836E+06 0.0000E+00

Pro kvalitativní posouzení je vhodnìjí obrázek, na nìm jsou zobrazeny izoèáry sloky napìtí

σ xx a to pro tøi

význaèné zatìovací stavy: maximum zatíení v prvním cyklu, odlehèený stav (rozloení zbytkových napìtí) a opakované zatíení ve druhém cyklu.


66

Z výsledkù je zøejmé, e k jistému pøizpùsobení nosníku dolo: Pøi elastickém zatìování nosníku bez zbytkových napìtí by meze kluzu bylo dosaeno silou

Fy=σK

B H2 = 8.333[N] , která je mení ne síla 6L

F y = 11.3[N] ; èást elementu è.30 blií povrchu nosníku se do

pøedepsaná ji v prvním zatìovacím kroku

plastického stavu opravdu dostala (viz. protokol výsledkù, krok 1). Zatíení naprosto stejnou silou v následujícím zatìovacím cyklu (krok 6) ji vyvolá pouze elastické deformace. To je dobrá zpráva; ta patná souvisí s tím, e z uvedených výsledkù nelze dostateènì pøesnì vypoèítat skuteènou hodnotu souèinitele pøizpùsobení, protoe chybí informace o míøe zpevnìní. Jediné co víme je to, e pøizpùsobená konstrukce se jetì pøi zatíení silou F 6 = 11.3 [N] chová elasticky, zatímco v následujícím kroku ( F 7 = 13 [N] ) se dostává opìt do plastického stavu. Souèinitel pøizpùsobení se tedy bude nacházet nìkde v mezích

kp=

F7 F6 = 1.56 . = 1.36 a 8.333 8.333

Závìreèná poznámka: Smysl hledání souèinitele pøizpùsobení tkví v tom, e umoòuje stanovit takové poèáteèní zatíení konstrukce (napø. zkuební tlak), které vede k "nejlepímu" pøizpùsobení. Vpodstatì jde o zjitìní mezního stavu únosnosti (tj. zatíení, pøi nìm se zplastizuje celý prùøez stìny).

4.3

Teèení materiálu - creep Zvlátì pøi vyích teplotách ( 370

o

C u uhlíkových a 430 o C u vysoce legovaných ocelí) dochází k

èasové zmìnì deformace i pøi konstantním zatíení; to znamená, e materiál se chová jako vazká kapalina. Pøesnìji øeèeno viskoelastická, protoe výsledná deformace je souètem elastické deformace (vypoètené pro dané napìtí z Hookova zákona) a èasovì promìnné sloky (v oblasti tzv. sekundárního creepu je tato sloka deformace pøímo úmìrná èasu a analogie s kapalinami je dokonalá). Jakkoliv jsou jevy elastoplasticity a elastovazkosti principiálnì odliné, modelují se èasto stejným zpùsobem. Vzpomeòme na základní konstitutivní rovnici newtonských nestlaèitelných kapalin:

kde

t t σ ′ = 2 µd,

t t σ ′ [Pa] oznaèuje deviátor napìtí, d [1/s] rychlost deformace a µ [Pa.s] je dynamická viskozita (v


67

pøípadì kapalin nenewtonských není konstantní, ale závisí na invariantech rychlosti deformace). Pro modelování creepu se pouívá právì tato rovnice v pøírùstkovém tvaru

dt r r σ ′, dεc= 2µ co je ovem prakticky toté co Prandtlùv Reussùv model plasticity (117), který také tvrdí, e pøírùstek nevratné deformace je úmìrný deviátoru tenzoru napìtí. Pøedstavitelem tìchto modelù creepu je napø. Nortonùv model mocninové závislosti mezi rychlostí deformace a napìtím (je úplnì stejný jako mocninový model kapalin se dvìma parametry: koeficientem konzistence a indexem toku). Program COSMOS pro analýzu creepu pouívá rovnì mocninový model4, ale v ponìkud jiném významu (uvauje toti explicitnì i závislost konzistence materiálu na èase, tzv. èasové zpevnìní):

r dεc r = c0 c2 σ cM1-1 t c2-1 σ ′. dt V modelu (120), který má tøi parametry ( c0 ,c1 ,c2 ) oznaèuje

σ M intenzitu napìtí.

Vzhledem k podobnosti modelù creepù a plasticity, je postup výpoètu MKP témìø stejný, pouze místo zatìovacích krokù jsou kroky èasové5. Demonstrujme to na pøíkladì zatìování tlustostìnného válce vnitøním pøetlakem (pouijeme dvourozmìrné izoparametrické elementy PLANE2D s 8-mi uzlovými body, rotaènì symetrickou variantu matice tuhosti a mocninový model creepu). Geometrie a èasový prùbìh zatìovacího tlaku jsou uvedeny na obr. 22. Text programu: SF4CORD,1,.16,0,0,.25,0,0,0.25,1,0,.16,1,0,ukázka toho, e povrch i vztané body lze vytvoøit jedním pøíkazem EGROUP,1,PLANE2D,0,0,1,0,0,0,1, volba elementu 1-axisymetrický, 1-creep MPROP,1,CREEPC,6.4E-18,4.4,1, koef. mocnin. modelu c0,c1,c2 (c2=1, tedy bez èasového zpevnìní) MPROP,1,EX,20E6, Modul prunosti E MPROP,1,NUXY,.3, Poissonova konstanta M_SF,1,1,1,8,3,1,1,1, sí 3 elementù s 8-mi uzly DND,1,UZ,0,18,1,UY,RX,RY,RZ,, pøedepsané posuvy a natoèení CURDEF,TIME,1,1,0,365,2.2,365,2.6,-365,5.6,-365, zatìovací funkce ACTSET,TC,1, aktivace zatìovací funkce TIMES,0,5,.1, t0,tmax,èasový krok A_NONLINEAR,S,1,200,200,0.001,0,N,0, statika, 200 iterací PCR,3,1,3,1,1, pøetlak na vnitøním povrchu (dán prùbìhem pøedepsané funkce) NL_SOLN,0,0, modifikovaná Newtonova Raphsonova metoda øeení R_NONLINEAR sputìní výpoètu

4

ε c = c0 e 5

Je jetì jedna alternativa, exponenciální model teèení, který lze v jednorozmìrném pøípadì popsat rovnicí

c1σ

σ

c4

[1- ec2( c3 ) t ]+ c5 ec6 σt .

V programu COSMOS vak u ádného typu prvku nelze uvaovat creep a plasticitu souèasnì.


68

Na obrázku jsou znázornìny základní výsledky; èasové prùbìhy radiálního posunutí uzlù 8 a 11 na vnitøním a vnìjím povrchu válce (to je standardní operace) a analogické prùbìhy intenzit napìtí (graf tohoto typu ji není podporován programem COSMOS a soubor dat intenzit napìtí bylo nutné vytvoøit pomocí vlastního programu na zpracování výsledkového souboru). Z výsledkù øeení je patrné, e teèení oceli se projevilo blahodárnì, protoe "odbouralo" vysoké pièky intenzity napìtí na vnitøním povrchu válce a pøispìlo k celkovému vyrovnání napìtí po prùøezu. Poznámka: Stejnì pøíznivì se teèení materiálu projevuje u potrubních systémù relaxací dilataèních napìtí.

69


Algoritmizace pøírùstkové metody V tomto odstavci se krátce zmíníme o numerických metodách elastoplastických výpoètù, které jsou ostatnì podobné metodám analýzy velkých deformací i teèení materiálù. Témìø vechny prakticky pouívané postupy jsou variantami Newtonovy metody (metody teèen); princip øeení je patrný z obr.23 na nìm je znázornìn schematický prùbìh závislosti zatíení F na posunutí q. Kadý zatìovací krok (pøírùstek zatíení ∆F ) se skládá z tìchto fází 1 1-2

r r q 1 , F 1 (existuje rovnováha) t tT t t -výpoèet matice teèné tuhosti K t1 = ∫ B Dep1 BdV a -výchozí stav

øeení soustavy rovnic pro pøírùstek posuvù pøi pøedepsaném zatíení 2-3

t r r K t1 ∆ q 1 = ∆ F 1 .

-pro takto vypoètenou deformaci se opraví napìtí tak, aby leela na ploe plasticity; tím ovem bude poruena rovnováha vnitøních sil odpovídající

3-4

-nerovnováné síly pøírùstek

r ∆ F2 .

t r r ∆ F 2 pøedstavují zatíení pro dalí iteraci K t2 ∆ q 2 = ∆ F 2 . Pokud je vypoètený

r ∆ q 2 pøíli velký, vracíme se k pøedchozímu kroku 2-3.

Výe popsaný postup je klasická Newtonova Raphsonova metoda øeení soustavy nelineárních rovnic; v programu COSMOS oznaèovaná zkratkou NR. Implicitnì je COSMOSem pouívána Modifikovaná Newtonova Raphsonova metoda (oznaèení MNR), která se lií tím, e teèná matice tuhosti se v kadém pøírùstkovém kroku poèítá jenom jednou (ve fázi 12) a v následujících iteracích (fáze 3-4) se nemìní. Úspora operací je vykoupena sníením rychlosti konvergence (napø. pøi øeení elastoplastického nosníku se MNR neosvìdèila). Tøetí COSMOSem nabízenou moností je metoda BGFS (Broyden, Goldfarb, Fletcher, Shanno) patøící do skupiny tzv. quazinewtonských metod, které vycházejí z aproximace matice teèné tuhosti, kterou poté v následujících iteracích zpøesòují. Metoda je povaována za velmi efektivní a jsou jí vybaveny jen ty nejlepí programy MKP (napø. NASTRAN).

5.

Vyhodnocení napjatosti Hlavním problémem analytických postupù dimenzování je volba zatíení navrené konstrukce.

Dùleitou roli zde hrají pøedpisy, které obvykle stanoví hlavní projektové zatíení, napø. výpoètový pøetlak pøi kontrole tlakové nádoby (ASME 1989). Specifikace dalích vlivù (dodateèné reakce, zatíení vìtrem, snìhem a seismicitou, dilatace, seismicita) je povaována za vìc projektanta (výrobce). Výe popsanými postupy MKP lze pro víceménì libovolnou konstrukci a libovolné zatíení stanovit v kterémkoliv místì konstrukce tenzor napìtí (tedy i intenzitu napìtí), tenzor deformace a posunutí. Otázka zní co si s tìmito èísly poèít? Staèí vybrat z tohoto seznamu prostì tu nejvìtí intenzitu napìtí a porovnat ji napø. s mezí


70

pevnosti pouitého materiálu? Ano, v øadì pøípadù to staèit mùe, nicménì u exponovaných konstrukcí, jako jsou napø. tlakové nádoby je nutné uvaovat i "charakter" napìtí (viz. 5.1), maximální pøípustné deformace (záleí na kontextu, obvykle 1%), rychlost teèení (napø. dle ASME maximálnì 0.01% za 1000 hodin) nebo "únavu materiálu" (rezervu ivotnosti, kapitola 5.2).

5.1

Kategorizace napìtí dle ASME Codu

Pøíkladem postupu vyhodnocování výsledkù MKP jsou Americké pøedpisy pro navrhování tlakových nádob ASME-Boiler and Pressure Vessel Code. Asi nejdùleitìjím prvkem tohoto pøedpisu je kategorizace napìtí na tøi základní typy, pro nì jsou stanoveny rùzné koeficienty bezpeènosti vzhledem k normativnímu dovolenému napìtí

σ D ≤ max(

σK σP , ) . Smyslem kategorizace je rozumný pøedpoklad, e projektant 1.5 2.4

bude schopen provést analýzu napìtí jen v elastickém oboru (tedy bez uvaování vlivu elastoplasticity).

1-primární napìtí - jde o "staticky urèitá" napìtí, která dostaneme pøímo z podmínek rovnováhy. Základní pøíklady primárních napìtí jsou:

σ m (kotlové vzoreèky); primární napìtí jsou nejnebezpeènìjí, protoe pøi pøekroèení meze kluzu by okamitì dolo k neomezené deformaci ( σ m ≤ σ D ) - membránová napìtí

- prostá ohybová napìtí; protoe víme, e pøekroèení meze kluzu napø. v krajních vláknech nosníku jetì nevyèerpá jeho únosnost, zvyuje se dovolené napìtí 1.5 krát ( σ b ≤ 1.5 σ D , tato hodnota odpovídá zvìtení ohybového momentu pøi plném zplastizování obdélníkového prùøezu nosníku). 2-sekundární napìtí - jde o "staticky neurèitá" zpravidla ohybová napìtí, vzniklá poruchou membránové napjatosti (pøechody válec-kuel, zmìna tlouky stìny,...) a dále teplotní napìtí vyvolaná zmìnou teplot podél plátì nádoby. Protoe v místech, kde tato napìtí pùsobí se pøedpokládá monost pøizpùsobení konstrukce (tzn. vznik výhodných zbytkových napìtí po prvních zatìovacích cyklech) zvyuje se dovolené napìtí na trojnásobek. 3-pièková napìtí - úzce lokální napìtí s dosahem kratím ne je tlouka stìny nádoby (proto jsou nezachytitelná výe uvedenými skoøepinovými prvky). Jde pøedevím o napìtí vrubová nebo napìtí vyvolaná zmìnou teplot po tlouce stìny. Pøi statickém dimenzování se vùbec neuvaují (s výjimkou posouzení køehkého lomu), ovlivòují pouze ivotnost konstrukce. Pozn.: Skoøepinové prvky (SHELLAX, SHELL3, ...) urèují vdy pouze souèet primárních a sekundárních napìtí, celková napìtí (tj. vèetnì pièkových) je teoreticky moné stanovit jen pomocí elementù PLANE2D (u 2D konstrukcí) nebo SOLID (u 3D). Urèitý problém spoèívá v tom jak tato souètová napìtí rozloit do tøí výe uvedených kategorií; u tlakových nádob se zpravidla pouívá následující postup: primární napìtí se poèítají jako prùmìrná napìtí po prùøezu, sekundární linearizací vypoèteného profilu napìtí a pièková tvoøí zbytek, viz. obr.24.

71


5.2

ivotnost konstrukce dle ASME Codu Únava konstrukce je podmínìna opakovanou plastickou deformací materiálu, k ní dochází i tehdy,

jsou-li vechna vypoètená nominální napìtí pod mezí kluzu. Jak to, e i v tomto pøípadì mohou vzniknout plastické deformace? Souvisí to s existencí mikrotrhlin, které jsou silnými koncentrátory napìtí by s malým dosahem. Tyto jevy lze vyetøovat v zásadì dvojím zpùsobem, analýzou mechanizmu íøení trhliny nebo experimentálnì na základì únavových zkouek. Jsou-li napø. defektoskopickými metodami lokalizovány konkrétní trhliny, dá se pomocí lineární lomové mechaniky (LLM) alespoò pøiblinì odhadnout, pøi jakém nominálním napìtí dojde k nestabilnímu íøení trhliny - lomu. Pøedstavme si, e v jinak homogenním materiálu se náhle objeví kruhová trhlinka, její rovina je kolmá na smìr pùsobícího nominálního tahového napìtí

σ , obr.25. Pøi vzniku trhlinky bylo tøeba investovat jistou

energii na pøekonání soudrnosti materiálu a celkem logicky lze pøedpokládat, e tato energie je úmìrná ploe trhlinky ( G s je konstanta úmìrnosti)

E s = G s π r 2 [J]. Souèasnì se vznikem trhliny dochází ke zmìnì elastické deformaèní energie tìlesa, protoe "ploky" trhliny se pùsobením napìtí polomìru trhlinky r, napìtí

σ rozevøou o jistou vzdálenost u, která je pøímo úmìrná

σ a nepøímo úmìrná modulu prunosti E: u = fσ

r . E

Vztah (122) (v nìm f je jen jakási konstanta úmìrnosti) vyplývá z rozmìrových úvah; zvate, e posunutí u nemùe být napø. kvadratickou funkcí polomìru vady, protoe by pak mìlo rozmìr energie odpovídající rozevøení trhliny je úmìrná souèinu síly ( π r

E i = fπ r 3

2

[ m2 ] . Zmìna deformaèní

σ ) a posunutí u:

σ2 . E

Jestlie je energie potøebná na poruení soudrnosti materiálu ( E s ) mení nebo rovna energii deformace ( E i ) dochází k rozbìhu trhliny

E i = fπ r 3 σ 2 / E ≥ E s = G s π r 2

σ

fr = K I ≥

1/ 2 G s E = K c , kde K c [Pa.m ]

je lomová houevnatost, mìøitelná materiálová charakteristika. Porovnáním vypoèteného souèinitele intenzity napìtí KI s namìøenými hodnotami lomové houevnatosti se tudí zjistí bezpeènost konstrukce s ohledem na rozbìh detekované vady. Uvedená analýza je pomìrnì hrubá a nevystihuje skuteènost, e napìtí v okolí trhliny zdaleka není konstantní; ve vrcholu (koøeni) trhliny dochází k výrazné koncentraci napìtí, které je pøímo úmìrné souèiniteli

72


intenzity napìtí KI (nezamìòovat se souèinitelem koncetrace napìtí!) a je nepøímo úmìrné druhé odmocninì vzdálenosti od koøene trhliny. Význam souèinitele KI je tedy té v tom, e umoòuje vyjádøit prùbìh napìtí v blízkosti koøene explicitními vztahy, viz. napø. Irwin 61. Pozn.: Kromì souèinitele intenzity napìtí KI, který odpovídá normálovému zatíení trhliny, se definují souèinitele intenzity napìtí i pro smyková zatíení rovnobìná (KII) nebo kolmá (KIII) s rovinou trhliny. V programu COSMOS je pro výpoèet souèinitele intenzity napìtí KI pouíván pøíkaz CRACK, èíslo_trhliny, uzel_1, uzel_2, uzel_3 kterým se definuje trhlina buï na rozhraní osmiuzlových elementù PLANE2D (mezi uzly 1 a 2, uzel_3 se nezadává), eventuálnì mezi dvacetiuzlovými elementy SOLID u trojrozmìrných modelù. Vypoètená hodnota K se vztahuje na pomìry panující ve pièce trhliny (uzel_1). Vrame se nyní zpìt k metodám "vìtìní" ivotnosti z výsledkù únavových zkouek. V pøedpisem ASME Code jsou uvedeny tzv. konstrukèní návrhové køivky ivotnosti urèující závislost amplitudy støídavého napìtí

σ A na poètu cyklù N kdy dochází k únavovému lomu (pro rùzné materiály a rùzné teploty jsou

pøirozenì jiné návrhové køivky, pøíklad je na obr. 26); informace v nich ukryté se vyuívají následujícím zpùsobem: 1.

Specifikuje se historie zatìování, pøièem jsou uvaovány pouze zmìny zatíení (bez ohledu na dobu

jejich trvání). Definují se základní události (EVENT: tlaková zkouka,...) a seznam zatìovacích stavù i=1,2,...,L (LOADING: pracovní tlak,...) s poètem opakování 2.

ni .

Pro kadý zatìovací stav se provede kompletní analýza MKP. Výsledkem jsou rozloení celkových

napìtí. 3.

Protoe se návrhová køivka týká zmìn a ne absolutních hodnot napìtí, proetøí se vechny moné

kombinace zatìovacích stavù I,J (poèet takových dvojic je L(L-1)/2). V místech konstrukce posuzovaných na ivotnost se vypoèítá rozdíl tenzorù napìtí testované dvojice

r r r σ IJ = σ I - σ J a urèí se hlavní napìtí

r σ IJ1 ≤ σ IJ2 ≤ σ IJ3 . Tenzor σ IJ je rozkmit napìtí (dvojnásobek amplitudy) take s návrhovou køivkou ( σ A - N ) se srovnává intenzita amplitudy napìtí definovaná vztahem ( σ IJ3 - σ IJ1 ) . Výsledkem je zjitìní, 2 e konstrukce mùe pøeít N IJ opakování zatìovacích stavù I a J. Podíl skuteèného poètu cyklù kombinace zatíení I,J k meznímu poètu cyklù je faktor dílèího pokození D IJ = n IJ / N IJ . 4.

Celkové pokození konstrukce je dle Minerova pravidla dáno souètem dílèích pokození. Poèítá se

takto: Smìrem od nejvyích amplitud napìtí se postupnì probírá seznam vech L(L-1)/2 kombinací zatíení a napø. u kombinace I,J se k celkovému pokození D pøiète dílèí pokození zakalkulovaných cyklù

D IJ = min( n I ,n J ) / N IJ . Poèet

n IJ = min( n I ,n J ) se potom odeète od zbývajích poètù cyklù zatíení I a J, která se

vyskytují v jetì nevyhodnocených kombinacích (a zatíení s mením poètem cyklù ze seznamu vypadne). Výsledkem souètu dílèích pokození je podíl vyèerpání celkové ivotnosti konstrukce

D= ∑

n IJ ( ≤ 1) . N IJ

Poznámka: Minerovo pravidlo tvrdí vpodstatì to, e nezáleí na èasové posloupnosti zatìovacích stavù. Je to hypotéza krajnì nedùvìryhodná (stejnì jako zdùvodnìní, e "pøesné poøadí zatìovacích stavù tak jako tak

73


obvykle neznáme"), nebo zcela ignoruje pøizpùsobování konstrukce v prùbìhu zatìování (a ji vlivem plastizace nebo teèení materiálu). Jako ilustraci uvedeme pøíklad únavového výpoètu rotaènì symetrické trysky, zatíené mìnícím se vnitøním pøetlakem, obr.26. Budeme pøedpokládat, e v zatìovací historii se vyskytují tøi zatìovací stavy odpovídající tlakùm

p1 = 0, p2 = 10000, p3 = 20000 [psi] (anglosaské jednotky, napø. psi, je výhodné

pouívat proto, e návrhové køivky v programu COSMOS jsou pøevzaty z pøedpisù ASME; lze ovem

σ A - N ). Pokud není pøekroèena mez plasticity staèí metodou koneèných prvkù analyzovat jen jediný zatìovací stav (napø. pro p2 = 10000 psi ), nebo napìtí pøi

nadefinovat i své vlastní køivky formou tabulky

jiném zatíení se získají pøenásobením pomìrem tlakù. Zatìovací historie trysky je tvoøena opakováním dvou událostí: tlaková zkouka (kombinace zatíení 1 a 3) s poètem opakování 2000 a nábìh na provozní tlak (kombinace 1 a 2) s poètem cyklù 10000. Pøi øeení programem COSMOS pouijeme pro statickou analýzu trysky prstencové prvky PLANE2D. Sí prvkù vytvoøíme tak, e nejprve popíeme obrys meridiánu trysky (CONTOUR) a jím vymezenou oblast (REGION) pokryjeme ètyøúhelníkovými prvky pomocí automatického generátoru (MA_RG). CRPCORD,1,3,0,0,8,0,0,8,2,0,6,2,0,5,3,0,5,7,0,4,8,0,4,10,0,3,10,0,3,0,& 0, Definice obrysových køivek èíslo 1 a 9 (souèasnì se zadávají souøadnice vztaných bodù) EGROUP,1,PLANE2D,0,1,1,0,0,0,0, Budou pouity elementy PLANE2D: varianta QM6, osová symetrie MPROP,1,EX,3E7, Modul prunosti E [psi] MPROP,1,NUXY,.3, Poissonova konstanta CTNU,1,9,1,4,2,1,3,1,4,1,5,3,6,1,7,3,8,2,9,8, Definice obrysu (vyjmenování køivek, z nich se skládá a poètù elementù na køivách) RG,1,1,1, Definice oblasti vymezené jediným obrysem èíslo 1 MA_RG,1,1,1,5, Generování trojúhelníkových prvkù (5 "vyhlazovacích" iterací) MARGCH,1,1,1,Q,4,1,0.4, Konverze trojúhelníkù na ètyøuzlové (Q) ètyøúhelníky NCOMPRESS,1,50,1, Pøeèíslování uzlù ECOMPRESS,1,50,1, Pøeèíslování elementù DCR,1,UY,0,1,1,, Nulové posuvy v axiálním smìrù na køivce èíslo 1 PCR,9,10000,9,1,10000, Zatíení tlakem 10000 psi podél køivky 9 R_STATIC Statický výpoèet FT_EVENT,1,10000, FT_EVENT,2,2000, FT_LOAD,1,1,0, FT_LOAD,2,1,1,1, FT_LOAD,3,2,1,2, FT_LOAD,4,2,0, FT_LOC,1,17,1,1,1, ACTSET,LOC,1, A_FATIGUE,0,1,2, R_FATIQUE

Událost èíslo 1 (normální najetí), poèet cyklù 10000 Událost èíslo 2 (tlaková zkouka), poèet cyklù 2000 Zatìovací stav 1, náleí události 1, nulový tlak Zatìovací stav 2, rovnì událost 1, napìtí z výpoètu, násobení mìøítkem 1 Zatìovací stav 3, patøí události 2, napìtí z výpoètu, násobené mìøítkem 2 Zatìovací stav 4, událost 2, nulové zatíení Vyhodnocení ivostnosti jen v uzlu 17 (hrana trysky) Pouití standardní návrhové køivky (napìtí-poèet cyklù) dle obr. 26 Výpoèet na únavu

F A T I G U E C A L C U L A T I O N R E S U L T S:

CYCLES ALTER. PARTIAL LOADING (EVN) LOADING (EVN) USED/ALLOWED ------------- ------------- ------------------- -----1 1 3 2 2000 7413. 42034. .26979 1 1 2 1 8000 0.7820E+05 21017. .10231 Cumulative Fatigue Usage Factor = 0.372098

STRESS

FACTOR


74

Z uvedených výsledkù (týkají se pouze jediného nejnamáhanìjího uzlu na hranì trysky) je patrné, e bìhem pøedpokládané zatìovací historie dojde k vyèerpání ivotnosti trysky na 37%.

75


6.

Stabilita Stabilitu konstrukce je moné hodnotit dle prùbìhu celkové potenciální energie E; obecnì lze øíci, e

konstrukce je stabilní tehdy, kdy její celková potenciální energie je minimální (kdy má absolutní minimum, tj. kdy prostì neexistuje nìjaký jiný deformaèní tvar, kterému by odpovídala mení hodnota E). Naopak, kdy se konstrukce nachází ve stavu, který je sice rovnováný (odpovídá mu tedy lokální minimum energie E), ale vedle toho existuje jetì jiný rovnováný stav s mení hodnotou E (globální minimum), staèí malý popud k tomu aby konstrukce "pøeskoèila" do onoho energeticky "chudího" stavu.

6.1

Analytická øeení stability nosníkù a válcù Pro základní konstrukèní prvky jako je pøímý prut, válcová nebo kulová skoøepina existují analytická

øeení, která poskytují údaje o kritickém zatíení (tlakem, silou nebo momentem) pøi nìm dochází ke ztrátì stability. Jako pøíklad uveïme tlaèenou válcovou skoøepinu (trubku namáhanou osovou silou), která se mùe chovat rùzným zpùsobem v závislosti na její tíhlosti (pomìru L/R): Velmi tíhlé válce se z hlediska stability podobají tyèím a pro kritickou sílu platí známý Eulerùv vztah, odvozený ji v úvodní kapitole:

F krit =

π3 3 2 Es R L

(

L R ), > 2.57 R s

kde s je tlouka stìny. "Støednì tíhlé" válce se u chovají jako typické skoøepiny a dochází pouze k místní ztrátì stability ("prolomení" stìny). Kritická síla v tomto pøípadì nezávisí na délce a dokonce ani na polomìru válce, nýbr pouze na tlouce stìny s

F krit = 3.8E s2

(1.22

R L s ). < < 2.57 s R R

Teprve u "tlustých/krátkých" válcù se délka válce znovu uplatní, co je dùsledek spolupùsobení koncù válce:

F krit = 5.65E

R s3 . L2

Tyèe zatíené tlakem se podobají tlaèeným válcùm a i u nich nelze kritickou sílu vyjádøit jediným univerzálním vztahem (platným pro vechny tíhlosti); historie výzkumu problému stability tyèe je zajímavá, Gere [36]: První øeení

F krit = (

π 2 ) EJ, L

publikoval matematik Leonhard Euler ji v roce 1744 vpodstatì jako høíèku, její platnost ovìøil Considére a po 140 letech (1889). Teprve tehdy se zjistilo, e kritické síly poèítané dle Eulerovy rovnice souhlasí s

76


experimentem jen pro velmi tíhlé tyèe (cca L/R>100), jinak je skuteèné kritické zatíení mení. Vysvìtlení, zaloené na pøedstavì zplastizování tyèe, podal Engesser (nìmecký stavitel), který navrhl pouít v Eulerovì vztahu místo modulu prunosti E modul plasticity H. Prakticky souèasnì (píe se rok 1895) navrhuje polák Jasinski teorii redukovaného modulu, která vychází z pøedstavy, e pøi vyboèení tyèe bude "vnìjí" odlehèená strana v elastické oblasti, zatímco "vnitøní" vyboèením pøitíená strana v oblasti plastické; tomu odpovídá Eulerùv vztah s modulem prunosti, který je kombinací modulu elasticity a plasticity (odtud název teorie "redukovaného modulu"). Problémem vzpìru s uvaování plasticity a velkých deformací se nezávisle na svých pøedchùdcích zabýval i Theodore von Kármán (1908), leè pøesto a do roku 1946 nebylo zøejmé, která z teorií je správná. Teprve tehdy toti napsal americký profesor aeronautiky Shanley krátký jednostránkový èlánek s názvem "The column paradox" v nìm rozpory teorie plastického a redukovaného modulu objasòuje. Poznámka: Tzv. Tetmayerùv vztah, oblíbený v zemích koruny èeské, ádnou teorií není - je to prostì jen empirický "fit" namìøených dat.

6.2

Vyetøování stability MKP Výpoèet kritických zatíení MKP vychází z rovnic pro velké deformace konstrukce, tj. rovnic

pøírùstkového typu

r t t t r ( K L + K u + K σ )∆q = ∆F. - - - Kt - - r

Mez stability odpovídá stavu, kdy nekoneènì malá zmìna vnìjích sil ( ∆F ) zpùsobí koneènì velká posunutí

r r q . To nastává tehdy, je-li matice teèné tuhosti K t singulární, protoe potom existuje nenulové øeení ∆q i

pro nulový vektor pravé strany. Toto konstatování je souèasnì smìrnicí pro volbu korektní metody vyetøování stability konstrukce: staèí pouít výpoètový modul velkých deformací, zvolit monotónnì rostoucí zatíení a sledovat deformaèní odezvu; dosaení meze stability se projeví náhlým rùstem posunutí a zpravidla i numerickými potíemi s konvergencí Newtonovy metody. Uvedený pøístup je výpoètovì nároèný (jako vechny nelineární metody MKP), proèe se takøka výluènì dává pøednost metodì, která vychází z rovnice (129) v ní je vynechána matice tuhosti teèné tuhosti

t K u . Matice

t r K t pak pøestane záviset na posunutích a je pouze lineární funkcí napìtí σ . Oznaèíme-li

r t σ , K σ referenèní napìtí a referenèní matici geometrické tuhosti (odpovídající víceménì libovolnì t zvolenému referenènímu zatíení), bude matice geometrické tuhosti pøi λ -násobku zatíení λ -násobkem K σ .

symboly

Na mezi stability tedy musí dle lineární teorie stability platit

t t r ( K L + λ k K σ ) q k = 0.

Rovnice (130) je formulací vlastní úlohy, jejím øeením je obecnì N (=poèet stupòù volnosti) rùzných hodnot

r λ k (vlastních èísel) a stejný poèet odpovídajících vektorù zobecnìných posunutí q k (vlastních vektorù). Z tìchto N-øeení nás pøirozenì zajímá jen jediné, to s nejmením vlastním èíslem λ 1 násobkù zatíení

(protoe udává nejmení kritické zatíení).


77

Pro øeení vlastní úlohy byla vyvinuta øada metod, z nich uvedeme jen jedinou (Seide [9]): Inverzní iterace je metoda, umoòující stanovení právì jen jediného vlastního èísla (a vlastního vektoru). Spoèívá v opakovaném vyèíslení následujícího souèinu matic

t -1 t r(m-1) r q (m) q , 1 = λ 1 K L Kσ 1

t r r q 1T(m) K L q (m) 1 λ 1 = - r T(m) t r(m) , m = 1,2,... q 1 Kσ q 1

r q 1 (obvykle se volí jednotkový vektor). Z rovnice r(m) r r (131) je zøejmé, e pokud iteraèní proces konverguje (tj. kdy q 1 = q 1 ) je q 1 opravdu øeením vlastní úlohy pøièem výsledek nezáleí na poèáteèní aproximaci vektoru

(130). Co ale na první pohled patrné není: vlastní vektor vypoètený inverzní iterací odpovídá nejmenímu vlastnímu èíslu (kritickému násobku zatíení), co právì pøi øeení stabilitní úlohy potøebujeme. Pøíklad: Rotaènì symetrická nádoba, její meridián je znázornìn na obrázku 27 je zatíena referenèním vnìjím pøetlakem p=1 [Pa]. Geometrický model (za povimnutí snad stojí pouití pøíkazu CRSPOLY, který umoòuje velmi jednodue kreslit èáry sloené z pøímých úsekù a krunic), volba typu prvkù (SHELLAX) a materiálových vlastností (z databáze), okrajových podmínek (nulové posuvy ve smìru z /obvod/ ve vech uzlech a nulové posunutí ve smìru osy nádoby v jediném uzlu 1) jsou popsány v následujícím programu: CRSPOLY,1,0.3,0,0,L,0.3,0.01,0,A,0.29,0.02,0,A,0.28,0.03,0,L,0.27,0.19& ,0,A,0.26,0.2,0,A,0.19,0.19,0,A,0.13,0.18,0,A,0.12,0.19,0,A,0.11,0.2,0& ,A,0.1,0.19, 0,L,0.1,0.19,0, úseèky L, kruhové oblouky A (s hladkými pøechody!) EGROUP,1,SHELLAX,0,0,0,0,0,0,0, elementy SHELLAX RCONST,1,1,1,1,.001, tlouka stìny 1 mm PICK_MAT,1,PC_STEEL,SI, pouitý materiál (uhlíková ocel); jednotky SI M_CR,1,3,1,2,2,1, generátor sítì na køivkách 1 a 3 (2 elementy/køivku) M_CR,4,4,1,2,6,1, dtto pro køivku 4 (6 stejnì dlouhých dvouuzlových elementù) M_CR,5,5,1,2,2,1, dtto pro køivku 5 (2 elementy) M_CR,6,7,1,2,4,1, dtto pro køivky 6 a 7 (po 4 elementech) M_CR,8,11,1,2,2,1, dtto pro køivky 8 a 11 (po dvou elementech) NMERGE,1,41,1,0.0001,0,1,0, slouèení uzlù na rozhraní køivek NCOMPRESS,1,41,1, pøeèíslování uzlù DND,1,UY,0,1,1,, nulový posun ve smìru osy symetrie uzlu 1 DCR,1,UZ,0,11,1,RY,RX,, nulové obvodové posuvy a pøísluná natoèení PCR,1,-1,10,1,-1, tlak p=-1 (vnìjí pøetlak, viz. obr.27) A_BUCKLING,1,I,16,0,0,0,1e-05,0,1e-06,0, metoda øeení (inverzní iterace), 16 iterací R_BUCKLING analýza (výpoèet vlastního èísla a vlastního tvaru) Ve výsledcích je dùleité vlastnì jen jediné èíslo, první vlastní hodnota: B U C K L I N G E I G E N V A L U E (S) by INVERSE POWER METHOD EIGENVALUE EIGENVALUE NUMBER 1 -0.3382092E+06


78

Z výsledkù je zøejmé, e ke ztrátì stability dojde pøi tlaku 338 kPa, ale ne vnìjím nýbr vnitøním kdy se "vyboulí" dno nádoby (by vyztuené lemem). První vlastní tvar deformovaného meridiánu je uveden rovnì na obrázku 27.


7.

79

Dynamika Dynamiku konstrukcí je zpravidla tøeba vyetøovat v pøípadech, kdy se uplatní:

- vliv vìtru na nádoby a potrubní sítì (periodické odtrhávání vírù) - seismicita (zemìtøesení, vliv dopravy, technologie) - nedokonale vyváené pohony. Cílem je zjistit, zda posuzovaná konstrukce splòuje poadavky, kladené napø. na seismickou odolnost rùznými pøedpisy, napø. ANSI/ASME (USA), RTM 108.020.37-81 MEM (SSSR). Pøi provìøování se hodnotí nejenom napjatost, vyvolaná napø. maximálním výpoètovým zemìtøesením (na rozdíl od statické kontroly se uvaují ponìkud vìtí hodnoty souèinitelù bezpeènosti), ale i pracovní zpùsobilost konstrukce (tím se rozumí napø. kontrola maximálních výchylek kmitajícího potrubí, která nesmí zpùsobit kolize s okolím, pøeruení kabeláe, naruení funkènosti armatur /ventilù/, hermetiènosti pøírubových spojù ap.). Druhým úkolem je optimalizace konstrukce, tj. nalezení rozumného kompromisu mezi poadavky statického a dynamického výpoètu: dilataèní napìtí v potrubní trase lze napø. úèinnì sniovat zmenováním tuhosti (aplikací kompenzaèních smyèek, odstranìním pevného ukotvení), jenome to má za následek i sníení vlastních frekvencí trasy a zvýení dynamických napìtí.

7.1

Vlastní tvary a vlastní frekvence kmitání konstrukcí Popis kmitání je moné zaloit na d'Alembertovì principu (který praví, e dynamickou úlohu lze chápat

jako statickou pokud jsou respektovány setrvaèné síly) a vyjít tudí ze základních rovnic lineární statiky (47,48):

rT t tr r t r t tt Kq = F, K = ∫ BT DBdV, F = ∫ f NdV, V

V

r f [N.m-3 ] je vektor objemových sil. Uvaujeme-li pouze zrychlující síly platí r r r t d2u d2 q f = -ρ 2 = -ρN(x, y,z) 2 , dt dt

kde

kde

t r N je døíve ji zavedená matice tvarových funkcí a q(t) vektor zobecnìných posunutí uzlových bodù

modelované konstrukce. Dosazením (133) do "statické" rovnice (132) získáme soustavu obyèejných diferenciálních rovnic druhého øádu pro neznámé èasové prùbìhy uzlových posunutí:

r tT t t r tT t d2 q ∫ B DBdV q = - V∫ N ρNdV dt 2 V t t --- K ----- M ---t t kde K je matice tuhosti a M matice hmotnosti (závisí pouze na tvarových funkcích a hustotì materiálu).


80

Øeením této soustavy rovnic (soustava je lineární, homogenní a s konstantními koeficienty) jsou harmonické funkce

r r q(t) = q k eiω t , k

jejich dosazením do soustavy diferenciálních rovnic (134) dostaneme vlastní problém, který

je navlas stejný jako pøi vyetøování lineární stability:

t t r (K - ω 2k M) q k = 0 .

Ilustrativní pøíklad: Tyèový prvek, uvaujeme pouze posunutí ve smìru osy táhla x. Dvìma uzlovým bodùm odpovídají v tomto pøípadì pouze dvì uzlová posunutí

q T = ( u1 ,u2 ) . Pouijeme-li lineární tvarové funkce N(x)

x L  u1 t r u(x) = ( N 1 (x), N 2 (x)).  = N.q  u 2 t t dN 1 d N 2 ), D = E, B=( , dx dx N 1 (x) = 1-

x , L

N 2 (x) =

získáme tyto matice tuhosti a matice hmotnosti prvku: L

t tT t t K = ∫ B DBdV = S ∫ V

0

 dN 1     dx  ES  1 - 1  dN 1 dN 2     .(E). dx , dx  dx = L   1 1  dN 2     dx 

L  N 1 t tT t ρSL  2 1 M = ∫ N ρNdV = S ∫   .( ρ ).( N 1 , N 2 )dx =  . 6  1 2 V 0  N 2

Poznámka: Vimnìte si, e seètením vech prvkù matice M obdríme celkovou hmotnost tyèe

ρSL . Matice

hmotnosti vypoètená dle pøedpisu (138) odpovídá spojitému rozloení hmoty podél prvku (konzistentní matice hmotnosti). Kdybychom hmotnost táhla soustøedili do jeho dvou uzlových bodù dostali bychom soustøedìnou matici hmotnosti, která je diagonální (její diagonální prvky jsou souèty vech prvkù konzistentní matice hmotnosti v pøísluném øádku):

t ρSL  1 0  . Ms= 2  0 1 Smysl zavedení nekonzistentní matice hmot spoèívá v tom, e ponìkud zjednoduuje a zrychluje výpoèet vlastních èísel; u niích vlastních frekvencí není vliv nekonzistence na chybu výsledku pøíli výrazný.


81

Pouijeme-li pro øeení vlastních frekvencí (tj. pro øeení vlastní úlohy (136)) metodu inverzní iterace,

ω 1 (a odpovídající vektor amplitud posunutí uzlových bodù r q 1 ). Je zajímavé, e inverzní iterace není nic jiného ne Stodolova metoda výpoètu frekvence vlastních kmitù,

obdríme jen jedinou (nejmení) vlastní frekvenci kterou pouívá i ÈSN 690010.

Potøebujeme-li stanovit i nìkterou vyí vlastní frekvenci (víme, e ná diskretizovaný model má N rùzných vlastních frekvencí) lze pouít inverzní iteraci s frekvenèním posunutím

t t r [K - ( ω 2 + ω 2s )M]q = 0 t t -1 t r(m) r q (m+1) = ω i2 (K - ω 2s M ) M q i , i

jejím øeením je vlastní vektor a vlastní frekvence

ωs :

ω i , která je nejblíe zvolené frekvenci ω s .

Poznámka: Iteraèní výpoèet spoèívá v opakovaném vyèíslování pravé strany rovnice (140), kde se vak kromì hledaného vlastního vektoru vyskytuje i neznámá vlastní frekvence. I kdy se to moná bude zdát

ω i dosadit cokoliv, protoe cílem iterací je zjistit jen "smìr" a ne velikost r vlastního vektoru q i (obvykle se stejnì normalizuje na jednotkovou délku). Vlastní frekvenci, která odpovídá podivným, lze za tuto hodnotu

nalezenému tvaru kmitù, zjistíme a výpoètem Rayleighova kvocientu (pøesvìète se vynásobením rovnice (140) transponovaným vektorem

rT t r q q i K i ω i2 = - ω 2s + r T t r . qi M qi

r q i zleva a uvìdomte si, e qT.K.q i qT.M.q jsou skalární velièiny)

V pøípadì, e bychom chtìli získat nìkolik vlastních frekvencí najednou, mùeme pouít nìkterou z metod, nabízených napø. programem COSMOS (inverzní iterace je implicitní): - Iterace podprostoru; zobecnìní metody inverzní iterace na skupinu vlastních vektorù - Lanczosova metoda; jedna z nejefektivnìjích iteraèních metod, hledající vlastní frekvence jako koøeny charakteristického polynomu - Jacobiho metoda; klasický postup iteraèního hledání vech vlastních vektorù najednou (tato metoda je vhodná pro úlohy s relativnì malým poètem stupòù volnosti). Pøíklad: Výpoèet vlastních tvarù kmitù a vlastních frekvencí skøínì pøevodovky. Výpoèty tohoto typu se uplatní napø. v oblasti bezdemontání diagnostiky: Opotøebení ozubení pøevodovky nebo loisek je moné monitorovat piezoelektrickými akcelerometry tj. z namìøené hladiny vibrací (a úrovnì jednotlivých harmonických spektra). Aby vak bylo moné správnì navrhnout umístìní vibraèního detektoru a frekvence, které bude pøipojený analyzátor vyhodnocovat, musíme znát vlastní frekvence samotného zaøízení (zde napø. skøínì pøevodovky). Geometrii symetrické skøínì pøevodovky silnì zjednoduíme (obr. 28) a omezíme se jen na model jedné poloviny (to by bylo naprosto korektní pøi statickém výpoètu, kdy bychom poèítali napìtí vyvolané silami v ozubených kolech, ale pøi dynamickém výpoètu to znamená, e výsledkem budou toliko symetrické tvary


82

kmitù a odpovídající vlastní frekvence). Pouijeme skoøepinové prvky SHELL4 (Kirchhoffovské desky) ve variantì QUAD4 (kombinace 4 trojúhelníkù DKT), a pro výpoèet tøí nejmeních vlastních frekvencí metodu iterace podprostoru: CRPCORD,1,0,0,0,0.39,0,0,0.39,0.24,0,0,0.24,0,0,0,0, definice obdélníka (øezu skøínì) v rovinì z=0 (ose symetrie skøínì) CRFILLET,5,1,2,.05,1,0, zaoblení rohù s polomìrem R=0.05 m CRFILLET,6,2,3,.05,1,0, CRFILLET,7,3,4,.05,1,0, CRFILLET,8,4,1,.05,1,0, CRPCOORD,9,0.05,0.05,0.35,0.34,0.05,0.35,0.34,0.19,0.35,0.05,0.19,0.35& ,0.05,0.05,0.35, definice obdélníka v rovinì z=0.35 m (boèní strana skøínì) SFEXTR,1,8,1,Z,.3, vytvoøení povrchu taením obrysu skøínì o 0.3 m ve smìru z CRJOIN,29,10,20,19,25,1,1, prostorová køivka è.29, která je teènou køivek 10 a 20 spojuje body 19 a 25 CRJOIN,30,11,18,19,24,1,1, tímto zpùsobem se definují hladké pøechodové køivky ve vech "rozích" skøínì CRJOIN,31,21,12,26,20,1,1, CRJOIN,32,23,11,27,20,1,1, CRJOIN,33,24,9,28,17,1,1, CRJOIN,34,12,15,17,22,1,1, CRJOIN,35,10,14,18,21,1,1, CRJOIN,36,17,9,23,18,1,1, SF4CR,9,9,10,11,12,0, definice povrchu è. 9 ze ètyø køivek: 9,10,11,12 (rovinný bok skøínì) SF4CR,10,19,31,11,29,0, podobnì dalí rovinné povrchy SF4CR,11,22,33,12,32,0, SF4CR,12,9,34,13,35,0, SF4CR,13,10,36,16,30,0, SF3CR,14,27,32,31,0,zaoblené rohy definuji jako plochy vymezené tøemi køivkami SF3CR,15,33,28,34,0, SF3CR,16,36,35,25,0, SF3CR,17,29,30,26,0, EGROUP,1,SHELL4,1,0,0,0,0,0,0,0, volba typu elementu SHELL4 RCONST,1,1,1,6,.005,0,0,0,0,0, tlouka stìny skøínì s=0.005 m M_SF,1,17,1,4,3,3,1,1, parametrické generování sítì (kadý povrch 3 x 3 elementy) NMERGE,1,612,1,0.0001,0,1,0, slouèení uzlù na rozhraní povrchù NCOMPRESS,1,605,1, pøeèíslování uzlových bodù PICK_MAT,1,PC_STEEL,SI, volba materiálu: PC_STEEL DND,44,ALL,0,44,1, okrajové podmínky: pøiroubování v uzlech 44 DND,41,ALL,0,41,1, a 41 (nulové posuvy i natoèení) DCR,1,UZ,0,8,1,RX,RY,, podmínky symetrie na køivkách 1 a 8 (tvoøí øez skøínì osou symetrie) A_FREQUENCY,3,S,16,0,0,0,0,1e-005,0,1e-006,0,0,0, volba metody (S-Subspace Iteration), 3 vlastní hodnoty R_FREQUENCY sputìní analýzy

Výsledky jsou uvedeny ve výstupním souboru s pøíponou OUT; vidíte, e obsahují údaje i o hmotnosti skøínì (18.44 kg) nebo momentech setrvaènosti, základem jsou vak tøi poadované vlastní frekvence NUMBER OF EQUATIONS . . . . . . . . . . . . . .(NEQ) = 876 MASS .184438E+02 |VOLUME .235442E-02 | IX .385042E+00 |IY .563600E+00 |IZ .513221E+00 | FREQUENCY FREQUENCY FREQUENCY PERIOD NUMBER (RAD/SEC) (CYCLES/SEC) (SECONDS) 1 .1059183E+04 .1685743E+03 .5932104E-02 2 .1194810E+04 .1901599E+03 .5258732E-02 3 .1456128E+04 .2317500E+03 .4314995E-02

Geometrie skøínì a tvary prvních tøí vlastních kmitù jsou uvedeny na obrázku 28.


7.2

83

Vynucené kmitání a seismicita

Pøípady vynuceného kmitání s buzením vnìjí silou (tlakem,...) nebo pohybem urèitého bodu (podpìry) lze popsat soustavou obyèejných diferenciálních rovnic, která je celkem pøímoèarým zobecnìním (134)

t d 2 qr t dqr tr r + Kq = F(t), M 2 + C dt dt t kde C je novì zavedená matice tlumení, která se obvykle (napø. v COSMOSu) aproximuje lineární kombinací

matice tuhosti a matice hmot

t t t -1 C = α M M + α K K, a α M [ s ], α K [s] jsou Rayleighovy koeficienty tlumení (jejich hodnoty se pro

konkrétní konstrukci zjiují experimentálnì, protoe to nejsou výhradnì materiálové charakteristiky i kdy vyjadøují vliv vnitøního tlumení pouitého materiálu; podstatné je to, e experimenty potvrzují pøedpoklad frekvenèní nezávislostí tlumení). Pro zpøehlednìní výkladu ji v dalím textu t tlumící èlen

C uvaovat nebudeme (i kdy Rayleighovo tlumení (143) nepùsobí ádné zásadní komplikace).

Harmonická budící síla (HARMONIC_ANALYSIS) Nejjednoduí a souèasnì nejbìnìjí je pøípad harmonické budící síly Pøíklad: Válec (komín, kolona, potrubí), pøíènì obtékaný proudem vzduchu, je rozkmitáván silou, která pùsobí kolmo na smìr proudìní. Pøíèinou tohoto jevu je periodické odtrhávání vírù v úplavu (Kármánova vírová stezka); bylo ovìøeno, e frekvence odtrhávání vírù je pøímo úmìrná rychlosti proudìní U ∞ a nepøímo úmìrná prùmìru obtékaného válce D. Budící sílu, pùsobící na èelní plochu S=D.H, lze aproximovat harmonickou funkcí (Køupka [34])

F(t) = kde

1 U∞ , ρ c D S U 2∞ sin ωt, ω = 2 10πD

c D ≈ 0.3, U ∞ ∈(8 - 14 m / s) a ρ je hustota vzduchu. Uvaujeme-li stacionární kmitání (tj. stav, kdy ji odeznìly vechny pøechodové dìje), budou odezvy


84

konstrukce (posunutí, deformace i napìtí) rovnì harmonické funkce (a se stejnou frekvencí jako budící síly)

r r F(t) = F ω eiωt ,

r r q(t) = q ω eiωt .

Dosazením (145) do soustavy diferenciálních rovnic (142) a krácení harmonickou funkcí získáme soustavu algebraických rovnic pro neznámý vektor amplitud zobecnìných posunutí

t t t r r (- ω 2 M + iωC + K) q ω = F ω , t r r Kω q ω = F ω , t r kde K ω je matice dynamické tuhosti. Ze soustavy (146) je patrné, e pøi nenulovém tlumení je q ω vektor

komplexních èísel, take odezva konstrukce se vzhledem k budícím silám fázovì posune (komplexní amplitudu lze chápat jako fázový posuv). Typická sekvence pøíkazù v programu COSMOS pro harmonickou analýzu je tato (COSMOS ovem nepouívá výe uvedenou metodu analýzy, nýbr postup zaloený na rozvoji øeení do tzv. modálních tvarù, /viz. následující odstavec/): PTYPE,5 NFREQ FCURVE FF R_DYNAM

(volba typu analýzy, harm. analýza) (modální analýza, poèet vlastních frekvencí) (závislost amplitudy budící síly na frekvenci) (zatíení ve specifikovaných uzlech)

Obecné prùbìhy budících sil (MODAL_ANALYSIS + TIME_HISTORY_ANALYSIS) Pro naprosto libovolnou pravou stranu (zatíení) lze soustavu diferenciálních rovnic (142) øeit pøímou integrací. Pøímá integrace znamená nahrazení derivací diferencemi a výpoèet odezev po jednotlivých èasových krocích; pøíkladem jsou napø. formule Eulerovy nebo Runge Kutta, které se vak pro øeení rovnic s druhými derivacemi pøíli nehodí a nahrazují se mnohem efektivnìjí Newmarkovou nebo Wilsonovou metodou, Bathe [7]. Velkou nevýhodou pøímé integrace rovnic zapsaných ve tvaru (142) je to, e v kadé z N-rovnic se vyskytuje vech N prvních i druhých derivací, take se v kadém integraèním kroku musí øeit velké soustavy algebraických rovnic. Obvykle (v programu COSMOS vdy) jest proto hledán zpùsob jak soustavu (142) transformovat na ekvivalentní soustavu vzájemnì nezávislých diferenciálních rovnic (kdy v kadé rovnici je toliko jediná hledaná funkce). Tuto transformaci umoòují vlastní vektory, které jsou øeením pøedchozího problému volného netlumeného kmitání (136). Zmínìné vlastní vektory (tvary vlastních kmitù

r q k ) mají tu zajímavou vlastnost, e jsou K-

ortogonální a M-ortonormální, co znamená, e platí

r tr q iT K q j = ω i2 δ ij ,

r t r q iT M q j = δ ij .

Proè to tak je: Napime (136) pro k=i a k=j a násobme tyto rovnice transponovanými vektory j-tého a i-tého vlastního kmitu zleva (výsledkem jsou skaláry!):


85

r tr r t q Tj K q i - ω i2 q Tj M qi = 0 r tr r t r q iT K q j - ω 2j q iT M q j = 0. Protoe matice tuhosti

t t rT t r rT t K i matice hmot M jsou symetrické, platí napø. q i K q j = q j K q i a vzájemným

odeètením rovnic (148) získáme vztah

rT t r (- ω i2 + ω 2j ) q i M q j = 0

rT t r ω i ≠ ω j musí být skalární hodnota q i M q j rovna nule a pouze pro r i=j mùe být nenulová. Nu a protoe na "velikosti" vektorù q i nezáleí lze je normalizovat tak aby r t r q iT M q i = 1 . t K-ortogonalitu vlatních kmitù lze dokázat analogicky (eliminací M z rovnic (148)). Seøadíme-li nyní vech N-vlastních vektorù vedle sebe, získáme ètvercovou matici Q(NxN) t r r r Q = ( q 1 ,q 2 ,..., q N ), která má vzhledem k (147) tyto vlastnosti tT t t t 2 tT t t t Q KQ = Ω , Q MQ = I, t t2 2 2 kde Ω je diagonální matice se ètverci vlastních frekvencí na diagonále ( Ωii = ω i ) a I je matice jednotková. který øíká, e pro rùzné vlastní frekvence

V praxi nikdy nepoèítáme vech N, ale pouze M<
Ωb maticí diagonální (a

z (151) plyne, e napø. Rayleighovo tlumení je speciálním pøípadem proporcionálního útlumu). Diagonální prvky matice

Ωb mají rozmìr frekvence a nazýváme je frekvencemi útlumu ω bi i-tého vlastního tvaru kmitání.

Pomìr frekvence útlumu a pøísluné vlastní frekvence je tzv. koeficient pomìrného útlumu i-tého vlastního

ζ i = ω bi / ω i (vechny koeficienty útlumu lze vyjádøit pomocí dvou Rayleighových koeficientù a vlastních frekvencí vztahem ζ m = ω bm / ω m = ( α M / ω m + α K ω m ) / 2 , který plyne pøímo z výe uvedených tvaru

definic). Matice vlastních vektorù

tr r q(t) = Qp(t).

t Q je základem transformace zobecnìných uzlových posunutí

Vztah (152) znaèí, e vektor uzlových posunutí

r q(t) je nahrazen lineární kombinací

r q i (pi(t) jsou koeficienty lineární kombinace); to je výhodné m.j. proto, e "skuteèné" kmity r r q(t) lze obvykle docela dobøe nahradit váeným souètem jen nìkolika prvních vlastních tvarù kmitù q i (tvary,

vlastních vektorù

odpovídající vyím frekvencím se zpravidla rychle tlumí). Místo N-diferenciálních rovnic pro funkce

q1 (t),..., q N (t) pak staèí øeit jen nìkolik málo rovnic pro modální funkce p1 (t),..., p M (t), M»N . Tyto rovnice získáme násobením pùvodní soustavy (142) (s vynechaným tlumením) transponovanou maticí vlastních


86

kmitù, substitucí (152) a vyuitím ortogonality matic tuhosti a hmot

r tT t t d 2 p tT t t r tT r Q MQ 2 + Q KQp = Q F, dt 2 d pm + ω 2m pm = f m (t), m = 1,2,... dt 2

Soustava rovnic (153) je ji separovaná, tj. tvoøená nezávislými rovnicemi,

které

lze

øeit

oddìlenì.

Tato

fáze

øeení

se

v

COSMOSu

skrývá

v

oddíle

NORMAL_MODE_ANALYSIS. Navazující øeení diferenciálních rovnic a zpìtná rekonstrukce èasových prùbìhù posunutí uzlù z èasových prùbìhù modálních posuvù se nazývá TIME_HISTORY_ANALYSIS. Protoe jsou rovnice (153) lineární a mají konstantní koeficienty, dá se najít jejich analytické øeení ve tvaru Duhamelova integrálu:

pm (t) =

1 ωm

t

∫ f ( τ )sinω (t -τ )dτ . m

m

0

Poznámka: (154) je nejjednoduí moný tvar odpovídající nulovým poèáteèním podmínkám. I pøi uvaování tlumení by modální funkce byly vyjádøeny Duhamelovými integrály typu (154), ale se sníenou frekvencí

ωm

(tlumení se projevuje sníením vlastních frekvencí, Bittnar [1]). V této fázi by ji mìlo být zøejmé jak øeit konkrétní pøípady, kdy je napø. zadán èasový prùbìh síly F(t) v nìkterém uzlu: 1)

Provede se analýza volných kmitù, jejím výsledkem je matice vlastních tvarù Q.

2)

Prùbìh "uzlových" síl F(t) se rozloí na jednotlivé "modální" síly

3)

Numerickou integrací Duhamelových integrálù (154) nebo pøímou integrací diferenciálních rovnic

(153) se vypoètou èasové prùbìhy "modálních" posuvù

r tT r f(t) = Q F(t) .

r p(t) (COSMOS pouívá pøímou integraci rovnic (153)

Newmarkovou nebo Wilsonovou metodou, které jsou zhruba rovnocenné; detailní popis tìchto algoritmù i jejich porovnání uvádí napø. Okrouhlík [24]). Èasový krok numerické integrace samozøejmì ovlivòuje pøesnost øeení; zpravidla by nemìl být delí ne desetina periody vlastních kmitù s nejvyí uvaovanou vlastní frekvencí. 4)

"Skuteèná" posunutí uzlových bodù

vlastních tvarù

r q(t) se rekonstruují násobením modálních posuvù maticí

tr r q(t) = Qp(t) . Z tìchto posunutí je pak moné vypoèítat deformace a napìtí stejným zpùsobem

jako døíve. Typická sekvence pøíkazù programu COSMOS: PTYPE,2 (volba typu analýzy, TIME_HISTORY) NFREQ (modální analýza, poèet vlastních frekvencí) TCURVE (èasový prùbìh budících sil...) FF (aplikovaných v uzlech ..) TSTEP (èasový krok, metoda integrace,...) R_DYNAM


87

Kinematické buzení (BASE_MOTION) V pøedchozím odstavci jsme analyzovali dynamiku systému s pevnými podporami, na který pùsobily èasovì promìnné síly nebo tlaky. Pøípad, kdy je buzení systému vyvolané pohybem základù, se øeí podobnì, ale je pøece jen o nìco sloitìjí. Omezíme se na pøípad, kdy vechny body ukotvení mají stejné èasové prùbìhy zrychlení, dané skalární funkcí

q&&b (t) (UNIFORM_BASE_MOTION); celková posunutí libovolného bodu

konstrukce mùeme vyjádøit jako souèet relativního posunutí qr(t) a posunutí qb(t), které odpovídá pohybu celé konstrukce jako tuhého tìlesa

r r r q(t) = q r (t) + I b qb (t),

kde vektor Ib urèuje, kterých stupòù volnosti se pohyb konstrukce jako celku týká (je-li napø. konstrukce buzena jen ve smìru x, budou prvky Ib rovny jedné u vech uzlových posunutí ve smìru x, jinak nuly). Protoe samotné posunutí qb(t) nevyvolá ani tøecí síly, ani zmìnu tvaru konstrukce, ale jen síly setrvaèné, získáme dosazením do obecné rovnice dynamické rovnováhy uzlù rovnici

r r t d 2 qr t d qr tr t r M + C + K q r = - M I b q&&b (t) 2 dt dt

Násobíme-li rovnici (156) zleva transponovanou maticí vlastních tvarù Q a nahradíme-li vektor uzlových

r

tr

relativních posunutí qr vektorem relativních modálních posunutí ( q r = Q pr ), získáme rovnici

r r r t d 2 pr tT t r t d pr t2 r && p q q&&b (t). I Γ + = M (t) = + 2 Q I Ω Ω b b r b dt dt 2

Vzhledem k tomu, e matice Ω jsou diagonální, pøedstavuje pøedchozí maticový zápis soustavu separovaných obyèejných diferenciálních rovnic pro relativní modální posuvy

d prm d 2 prm + 2 ω bm + ω 2 prm = - Γ m q&&b (t), 2 dt dt Symbolem

m = 1,..., M.

r Γ m je zde oznaèen tzv. participaèní faktor m-tého vlastního tvaru (vektor Γ urèuje jakým

zpùsobem se rozloí síly, vyvolané pohybem základù, na jednotlivé vlastní tvary). Povimnìte si, e participaèní faktor má rozmìr hmotnosti a pøedstavuje vlastnì podíl hmoty systému, která se uplatní pøi m-tém vlastním tvaru kmitání; Eukleidovská norma participaèního vektoru by správnì mìla být rovna hmotnosti celého systému a pøípadná odchylka je projevem zanedbání vyích tvarù kmitù. Øeení rovnic (158) je shodné s postupem, popsaným v pøedchozím odstavci (kmitání, vyvolané vnìjími silami). Nedeterministické budící síly - Seismicita (RANDOM_VIBRATION_ANALYSIS) Výe uvedené postupy pøedpokládaly, e je pøesnì znám èasový prùbìh budících sil nebo pohybu


88

podpor. Zjednoduené zpùsoby analýzy chování konstrukcí pøi chvìní základù vak nejsou zaloeny na konkrétních záznamech seismografù (èasových prùbìzích zrychlení) nýbr na skuteènosti, e budící síly mají náhodný charakter a k jejich popisu (a souèasnì k popisu odezvy konstrukce) jsou spíe ne èasové prùbìhy vhodnìjí frekvenèní (spektrální) charakteristiky. Základem je vyjádøení budících sil Fourierovými integrály

1 F(t) = 2π kde

∞

∫ F ( ω ) e dω , -iωt

ω

-∞

F ω ( ω ) je Fourierova transformace èasového prùbìhu budících sil, pro ni platí (Press [35]) ∞

Fω (ω ) =

∫ F(t)e

iωt

dt.

-∞

r F(t) je nahrazen souètem nekoneènì mnoha harmonických prùbìhù F ω e-iωt , kde r 2 * F ω je amplituda odpovídající frekvenci ω . Ètverec této amplitudy S F ( ω ) = F ω F ω =|F ω | je výkonová hustota (PSD Power Spectral Density) budicí síly pøi frekvenci ω a pro zcela náhodný signál je konstantní Interpretace: Èasový prùbìh

(co je pøedpoklad, který se osvìdèuje napø. právì u seismického buzení). Poznámka: Termín výkonová hustota je pøejatý z elektrotechniky, kde se výkonová hustota urèuje jako elektrický výkon signálu, který proel pásmovým filtrem (a obsahuje pouze frekvence v intervalu

ω ,ω + ∆ω ); pøesnìji, výkon signálu dìlený

délkou intervalu ∆ω . Poznámka: Velmi èasto (napø. v manuálech COSMOSu) se setkáváme s odlinou definicí výkonové spektrální hustoty signálu SF( ω ). Zní takto: Výkonová spektrální hustota èasového prùbìhu F(t) je rovna Fourierovì transformaci autokorelaèní funkce RF(t) signálu ∞

R F (t) =

∫ F( τ + t)F( τ )dτ ,

∞

S F (ω ) =

-∞

∫ R (t)e dt. F

iωt

-∞

I kdy tato definice PSD vyhlíí úplnì jinak ne pøedchozí, vyjadøuje toté díky identitì, která praví: Fourierova transformace autokorelaèní funkce nìjakého signálu F(t) je rovna absolutní hodnotì souèinu Fourierových transformací této (reálné) funkce F(t) ∞ ∞ ∞  iωt   iω (t+τ ) -iωτ     dt = (t) dt = τ τ τ τ τ τ F(t + )F( )d dt = F(t + )e F( )e d R e e F ∫ ∫ ∫ ∫ ∫   -∞ -∞  -∞ -∞  -∞ ∞

∞

iωt

∞  ∞  iωτ =  ∫ F( τ )e dτ   ∫ F( τ )e-iωτ dτ  = F ω ( ω ) F ω (-ω ) = F ω ( ω ) F *ω ( ω )  -∞   -∞  Hvìzdièka jako horní index Fourierovy transformace F oznaèuje komplexnì sdruenou funkci a protoe souèin komplexního èísla s tím samým, ale komplexnì sdrueným èíslem, je roven druhé mocninì absolutního hodnoty, je zøejmé, e obì definice výkonové spektrální hustoty jsou opravdu rovnocenné (poslední rovnost v pøedchozím vztahu je správná pouze pro reálné funkce F(t), co ale v daném kontextu platí).


89

Základem dalího postupu je snaha o nalezení vztahu mezi známou výkonovou spektrální hustotou budící síly a spektrální hustotou odezvy (posunutí konstrukce). Odezvu konstrukce na harmonickou budící sílu

r r r -iωt -iωt ji urèit umíme; vektor zobecnìných posunutí bude ve tvaru q(t) = q ω ( ω )e , a jeho F ω ( ω )e r dosazením do základní rovnice vynuceného kmitání získáme Fourierovu transformaci odezvy (posunutí q ω ) t t t r r (- ω 2 M - iωC + K) q ω ( ω )e-iωt = F ω ( ω )e-iωt t t t t r t r q ω ( ω ) = - K ω-1 F ω ( ω ), kde K ω = ω 2 M + iωC - K. Poznámka: Pøedchozí úvaha moná není pøíli srozumitelná. Alternativní odvození, vycházející pøímo z definice transformace (160), je toto: Fourierovu transformaci uzlového posunutí získáme aplikací Fourierovy

&& + Cq& + Kq = transformace na jednotlivé èleny výchozích diferenciálních rovnic (142) ( Mq

F ) a vyuitím

pouèky o transformaci derivací: ∞

∫

-∞

∞ +∞ dq iωt iωt - iω ∫ q eiωt dt = - iω qω ( ω ) e dt = [q e ] -∞ dt -∞

∞

d 2 q iωt ∫-∞ dt 2 e dt = - ω 2 qω ( ω ). Výsledkem bude opìt vyjádøení transformace odezvy vztahem (163). Ze známé Fourierovy transformace odezvy lze nyní urèit její výkonovou spektrální hustotu (pøesnìji matici výkonových spektrálních hustot

t Sq )

t t * -1 r* r T t -1 T t * -1 t t -1 T r* r T S q ( ω ) = q ω q ω = ( K ω ) F ω F ω ( K ω ) = ( K ω ) S F ( ω )( K ω ) . Rovnice (165) pøedstavuje poadovaný vztah mezi zadávanou spektrální výkonovou hustotou buzení a odezvy. Jakým zpùsobem se ale dostat zpátky od výkonové spektrální hustoty posunutí ke skuteèným (støedním nebo maximálním) posunutím uzlových bodù a tedy i k napìtím? Zamysleme se nad významem diagonálních prvkù matice výkonových spektrálních hustot; Sqii je dle definice PSD rovno Fourierovì transformaci autokorelaèní funkce posunutí qi(t) a její èasový prùbìh lze získat zpìtnou Fourierovou transformací (159)

1 (t) = ∫ q ( τ ) q ( τ +t)dτ = ∫ S ( ω ) e dω 2π ∞

Rqii

∞

i

qii

i

-∞

-iωt

-∞

∞

2 Rqii (0)= ∫ qi ( τ -∞

1 )dτ = ∫ S ( ω )dω . 2π ∞

qii

-∞


90

Rqii(0), støední hodnota kvadrátu amplitudy odezvy posunutí qi (Mean Square Response), je tudí mìøítkem výsledných posuvù a základním výsledkem analýzy (podobnì se poèítají rychlosti i zrychlení). Pøedloený postup má tu nevýhodu, e vyaduje vyèíslování velkého poètu integrálù (166) (ètvercová matice S má rozmìr N = poèet stupòù volnosti celé soustavy). Proto se v analogii s metodami pøímé integrace hledá monost, jak vyuít výsledky modální analýzy ke zmenení rozmìru matice výkonových spektrálních hustot pouze na hodnotu M (M = poèet uvaovaných vlastních tvarù kmitù). Východiskem je zavedení substituce (152), tj. q=Q.p a separace soustavy diferenciálních rovnic (152) vyuitím M-ortogonality matice vlastních tvarù Q

r r r t d2 p tT r t dp t2 r I 2 + 2 Ωb + Ω p = f = Q F(t) dt dt

N dpm d 2 pm 2 p f Qmi F i (t) . + 2 + = (t) = ω ω bm m m ∑ m dt dt 2 i=1

Stejným zpùsobem jako døíve (Fourierovou transformací rovnic (167)) obdríme Fourierovy transformace modálních odezev pm

pωm ( ω ) =

1 f ( ω ) = H m ( ω ) f ωm ( ω ), - ω - 2i ω bm ω + ω m2 ωm 2

kde Hm je pøenosová funkce (pøesnìji diagonální prvek matice pøenosù). Výkonová spektrální hustota (PSD) modální odezvy je souèin Fourierovy transformace modální odezvy (168) s odezvou komplexnì sdruenou, co v maticové notaci zapíeme takto (výsledná ètvercová matice PSD má dimenzi toliko M)

r r* t t t t* t* S p ( ω ) = H( ω )f( ω ) f ( ω ) H ( ω ) = H( ω ) S f ( ω ) H ( ω ) . Výkonovou spektrální hustotu (PSD) modálního buzení f vyjáøíme z PSD uzlových sil F

t ∞ tT r t rT R f (t) = ∫ Q F( τ ) F ( τ +t)Qdτ -∞

t tT t R F (t)Q

=Q

t tT t t S f ( ω ) = Q S F ( ω )Q . Analogicky lze vyjádøit PSD uzlových posunutí z ji vypoètených PSD modálních posuvù (169):

tt tT t S q ( ω ) = Q S p ( ω )Q .

Støední hodnoty ètvercù amplitud posunutí (MSR, hodnoty autokorelaèních funkcí pro t=0) stanovíme numerickou integrací PSD modálních posuvù, viz.(166) a podobnostní transformací

t t tT t Rq (t = 0) = Q R p (t = 0)Q .


91

V programu COSMOS se pro integraci PSD pouívá Gaussova dvou nebo tøíbodová integrace na subintervalech, které pokrývají uvaovaný rozsah frekvencí; uivatel programu zadává poèet subintervalù mezi sousedními vlastními frekvencemi a míru jejich zhuování v okolí vlastních frekvencí (èím je mení tlumení systému

ζ,

tím detailnìji by mìly být vystieny prùbìhy PSD v tìsné blízkosti vlastních frekvencí). COSMOS nabízí i úspornou variantu, která neuvauje vzájemné spolupùsobení modù (mimodiagonální prvky matice modálních PSD) a navíc pøedpokládá, e v okolí vlastních frekvencí je prùbìh PSD konstantní.

Typická sekvence pøíkazù: PTYPE,4 (typ analýzy, RANDOM_VIBRATION) NFREQ (poèet vlastních frekvencí) MDAMP (souèinitel tlumení jednotlivých modù kmitání) CURTYP,1 (typ funkce, frekvenèní závislost) FCURVE (zadání PSD) BSTYPE,BSDIR (typ kinematického buzení /posun èi zrychlení základù atd/) FRANGE (rozsah budicích frekvencí) R_DYNAM Analýza na základì spektra odezvy (SPECTRA_RESPONSES) Tento postup (jeden z nejèastìjích, pouívaný ve vìtinì pøípadù seismických analýz) je opìt zaloen na modální analýze a jeho princip je prostý: Výsledkem modální analýzy je vlastnì nahrazení pùvodní sloitì provázané mechanické soustavy, souborem nìkolika systémù s jedním stupnìm volnosti, elementárních oscilátorù (kombinací hmota+pruina+tlumiè). Pro daný èasový prùbìh zatíení (vibrací) je moné pomìrnì snadno stanovit tzv. spektrum odezvy elementárního oscilátoru (spektrum odezvy je závislost maximální odezvy /výchylky, rychlosti nebo zrychlení/ versus vlastní frekvence soustavy bod+pruina+tlumiè, pøièem parametrem tìchto køivek je souèinitel tlumení). Prùbìhy spektra odezvy oscilátoru pak nahrazují informace o konkrétním èasovém prùbìhu budicích sil (spektra jsou publikována v èasopisech nebo je prodávají specializované firmy) a umoòují odhadnout maximální odezvy konkrétní vyetøované konstrukce napø. jako souèet spekter odezev elementárních oscilátorù pro hodnoty frekvencí, které byly zjitìny modální analýzou systému. U tohoto postupu (stejnì jako pøi analýze zaloené na výkonových spektrálních hustotách) tedy odpadá nutnost èasové integrace diferenciálních rovnic, ale na druhé stranì zase není jednoznaèné, jakým zpùsobem nebo s jakou vahou by se mìla spektra odezev oscilátorù, odpovídající vlastním frekvencím soustavy, sèítat program COSMOS nabízí øadu metod: SRSS, CQC, NRL o nich se zmíníme pozdìji. Praktický postup výpoètu je následující: modální analýzou se vypoètou vlastní frekvence konstrukce a matice vlastních tvarù Q (kadý sloupec této matice odpovídá jednomu vlastnímu tvaru). Rozmìr této matice je N (poèet øádkù = poèet stupòù volnosti celé konstrukce) krát M (poèet vlastních èísel, obvykle podstatnì mení

ω m a jí korespondující modální tlumení ζ m odeèteme z grafu spektra odezvy oscilátoru hodnotu Sd( ω m ,ζ m ) a násobením participaèním faktorem (viz. odstavec, popisující kinematické buzení) získáme maximální modální posunutí odpovídající frekvenci ω m ne N). Pro kadou vlastní frekvenci


92

pm,max = Γ m S d ( ω m ,ζ m ) Posunutí v jednotlivých uzlových bodech, odpovídající m-té vlastní frekvenci, je dáno m-tým sloupcem matice Q (m-tým vlastním tvarem kmitù) krát modální posunutí, tedy

qi,m,max = Qi,m pm,max ,

i = 1,..., N m = 1,..., M.

Výsledné posunutí (nazývá se S.M.R., Structure Maximum Response) se zjiuje superpozicí odezev jednotlivých modù a to buï -jako souèet absolutních hodnot M

qi,max= ∑| qi,m,max |, m=1

-metodou SRSS (Square Root Sum of Squares)

qi,max=

M

∑q

2

m=1

i,m,max

-metodou CQC (Complete Quadratic Combination)

qi,max=

M

M

∑∑ qi,m,max qi,n,max r mn ( m=1 n=1

ωm ωn

)

kde rmn je vzájemný korelaèní koeficient (závisí na pomìru frekvencí a tlumení) -nebo metodou NRL (Naval Research Laboratory), která je vlastnì kombinací metody SRSS a metody superpozice absolutních hodnot spekter odezev

qi,max=| qi,M,max |+

M

∑q m=1

2

i,m,max

-q

2

i,M,max

,

kde qi,M,max oznaèuje spektrum "nejvìtí" odezvy v i-tém uzlovém bodì (nejvìtí mezi vemi modálními spektry odezev). Typická sekvence pøíkazù COSMOS: PTYPE,1 (typ analýzy, RESP_SPECTRA) NFREQ (modální analýza, poèet vlastních frekvencí) RPRINT,1 (co se má tisknout /posuvy, rychlosti,zrychlení) CURTYP,1,1 (typ køivky; frekvenèní závislost, kinematické buzení) FCURVE (bodovì zadané spektrum odezvy) BSTYPE, BSDIR (typ kinematického buzení /napø. zrychlení/ a jeho smìr) MCOMB (volba metody skládání odezev SRSS,CQC,NRL)


RESPR FRANGE R_DYNAM

93

(tisk spektra odezvy) (jednotky, v nich se vyjadøuje frekvence)

Pøíklad: Uvaujme potrubní úsek ve tvaru písmene S, který je na obou koncích pevnì vetknutý. Pokusíme se naznaèit postup analýzy chování tohoto potrubního úseku v pøípadì zemìtøesení. Tuto situaci bychom mohli modelovat tím, e bychom pøedepsali èasový prùbìh zrychlení (nebo posunutí èi rychlosti) v místech vetknutí a pouili pøímou integraci pro výpoèet odezvy (posunutí a napìtí v libovolné èísti potrubí). Vhodnìjí (a jednoduí) je vak pøedpokládat, e otøesy pùdy mají zcela náhodný charakter a v prùbìhu zrychlení jsou zastoupeny vechny frekvence. Místo závislosti zrychlení na èase proto zadáme závislost støední hodnoty kvadrátu zrychlení na frekvenci, tj. prùbìh spektrální výkonové hustoty PSD (pozn.: V programu COSMOS se pod pojmem PSD rozumí ji odmocnina této hodnoty, vpodstatì tedy støední hodnota zrychlení). Øeení ovem musí zaèít modální analýzou potrubního úseku, tj. výpoètem alespoò nìkolika základních vlastních frekvencí a tvarù kmitù (R_FREQUENCY). Výsledkem dalího øeení jsou prùbìhy PSD (posunutí, rychlostí, zrychlení, napìtí jako funkce frekvence) v libovolnì zvoleném místì potrubí. Postup øeení je uveden v následujícím programu; upozornìme na nìkteré delikátní body: 1. V pøíkazu PD_ATYPE je tøeba specifikovat rozsah analyzovaných frekvencí (zde 1 a 100); vodítkem mohou být výsledky modální analýzy, tj. hodnoty vlastních frekvencí potrubní sítì. 2. Pøíkazem PD_MDAMP se zadává hodnota vnitøního tlumení (jako pomìr skuteèného ke kritickému tlumení); je to vìcí empirie a obvykle se propoèítává nìkolik hodnot tlumení (zde byla pouita hodnota 0.1 stejná pro vechny uvaované vlastní kmity) 3. Pøíkazem PD_CURDEF je definována P.S.D. zrychlení chvìjících se základù; uvaujeme konstantní hodnotu a=1

m.s-2 .

CRPCOORD,1,0,0,0,10,0,0,20,0,0,20,10,0,30,10,0,30,10,0, Geometrie potrubního úseku (4-køivky) EGROUP,1,PIPE,0,0,0,0,0,0,0, Pouité elementy: PIPE RCONST,1,1,1,3,.2,.005,0, Prùmìr potrubí (0.2 m) a tlouka stìny (5 mm); nulový pøetlak PICK_MAT,1,PC_STEEL,SI, Materiál: uhlíková ocel M_CR,1,4,1,2,2,1, Generátor sítì : na kadém úseku 2 elementy NMERGE,1,12,1,0.0001,0,1,0, Slouèení uzlù na rozhraní køivek NCOMPRESS,1,12,1, Pøeèíslování uzlù DND,1,ALL,0,1,1, Vetknutí v konvových uzlech 1 a 9 DND,9,ALL,0,9,1 A_FREQ,10,S,30,0,0,0,0,1e-05,0,1e-06,0,0,0,0, Poadavek na výpoèet prvních 10-ti vlastních frekvencí (iterace podprostoru) R_FREQ Výpoèet vlastních frekvencí a vlastních tvarù PD_ATYPE,4,4,0,1,100,0,0,11,2,3,1e+10, Výpoèet výkonové spektrální hustoty pøi seismickém buzení v rozsahu frekvencí 1 a 100 PD_PRINT,1,1,1,0,0, Poadavek na tisk posuvù, rychlostí i zrychlení PD_MDAMP,1,1,4,.1, Vnitøní tlumení (souèinitel tlumení=0.1) PD_CURTYPE,1,1,1, Funkce è. 1 bude typu P.S.D. (spektrální hustota budicích kmitù: zrychelní - frekvence) PD_CURDEF,1,1,0,1,1500,1, Definice funkce è.1 pouze dvìma body: Konstantní P.S.D.=1 v intervalu frekvencí 0 a 1500 (rad/s) PD_BASE,1,0,0,1,0, Buzení nikoliv silou, nýbr pohybem ukotvení (vibrace ve smìru z) R_DYNAMIC Provedení výpoètu PD_RXYSET,1,5,3, Poadavek na vykreslení P.S.D. v uzlu è.5 (zrychlení ve smìru z) PD_PREPARE,1 Pøíprava souborù pro vykreslení grafu INITXYPLOT ACTXYPLOT,0,1,FREQUENCY,ACCELERATION,5,3,12,5N3 Definice prùbìhù v grafu XYPLOT Vlastní vykreslení prùbìhu spektrální hustoty odezvy v uzlu è. 5


WCREATE,2; ACTNUM,ND,1; DPLOT; EPLOT;

94

Ukázka toho jak vytvoøit 2 okénka Aktivace èíslování uzlù Oznaèení míst ukotvení Vykreslení elementù

Výsledný prùbìh PSD zrychlení v uzlu 5 (ohyb potrubní smyèky) je uveden na obrázku 29 Z prùbìhu PSD je dobøe patrné jak se rovnomìrnì rozloená energie budicích kmitù "pøelila" do oblasti malých frekvencí (cca 5 rad / s).

8.

Teplotní pole Teplotní pole v analyzovaných aparátech (tuhých tìlesech) popisuje Fourierova Kirchhoffova rovnice

vedení tepla:

ρ cp

∂T ∂ ∂T ∂ ∂T ∂ ∂T = (λx )+ ( λ y )+ ( λ z )+ Q, ∂t ∂x ∂x ∂y ∂y ∂z ∂z

s pøíslunými okrajovými podmínkami na povrchu tìlesa (pøedepsané teploty, hustoty tepelného toku q nebo souèinitel pøestupu tepla

α ).

Metoda koneèných prvkù je zpravidla zaloena na hledání funkce (tøeba prùbìhu posunutí), která minimalizuje jistý funkcionál (tøeba celkovou potenciální energii). Tento princip lze pouít i pøi výpoètu teplotního pole T(t,x,y,z), protoe lze ukázat, e funkce, která vyhovuje pøedepsaným okrajovým podmínkám prvního druhu (tj. pøedepsaným teplotám na pøísluné èásti hranice

E T (T) = ∫ V

S T ) a která minimalizuje funkcionál

t ∂T 1 T2 T ) ( ( ∆T Λ∆T - QT + ρ c p T )dV + ∫ qTdS + ∫ α ( - TT ∞ )dS ∂t 2 2 S Sα q

 λ x 0 0    Λ = 0 λ y 0    0 0 λ z


95

splòuje Fourierovu Kirchhoffovu diferenciální rovnici (179) vèetnì okrajových podmínek s pøedepsanou

S q (q oznaèuje sloku vektoru hustoty tepelného toku kolmou k povrchu) i na èásti hranice S α kde je pøedepsán souèinitel pøestupu tepla α a teplota vnìjího prostøedí T ∞ .

hustotou tepelného toku na hranici

Chováte-li k výe uvedenému tvrzení nedùvìru a nebojíte-li se napø. Gaussovy vìty, mùete ekvivalenci F.K. rovnice a variaèního principu (180) ovìøit následujícími úvahami: 1.

Pøedpokládejme, e T(t,x,y,z) je funkce, která je øeením naeho variaèního problému, tj. nabývá

S T a souèasnì minimalizuje funkcionál E T (T) . Zvolíme si variaci (poruchu) tohoto teplotního pole jako souèin malé skalární hodnoty ε a víceménì libovolné funkce W(x,y,z): δT = εW . Protoe "variované" teplotní pole T + δT musí splòovat okrajové podmínky prvního druhu na hranici S T je tøeba aby funkce W(x,y,z) byla na této èásti hranice nulová. Pro zjednoduení budeme uvaovat konstantní tepelnou vodivost materiálu λ x = λ y = λ z = λ a vyjádøíme hodnotu funkcionálu pøedepsaných hodnot (teplot) na èásti hranice

variovaného teplotního pole

∂T 1 T T T E T (T + εW) = ∫ { λ [ ( ∆T ) ∆T + ε ( ∆T ) W + ε 2 ( ∆T ) ∆W] - Q(T + εW)+ ρ c p (T + εW) }dV + ∂t 2 V + ∫ q(T + εW)dS + ∫ Sq

2.

Sα

1 2 ε2 2 α ( T + εTW + W - T T ∞ - εW T ∞ )dS. T 2

E T jen jako funkci parametru ε a urèeme její derivaci

Chápejme nyní

∂ E T (T + εW) ∂T T T = ∫ [ λ ( ∆T ) ∆W + 2λε ( ∆T ) ∆W - QW + ρ c p W ]dV + ∂ε ∂t V + ∫ qWdS + ∫ α (TW + εW - T ∞ W)dS Sq

Sα

Pokud T opravdu minimalizuje funkcionál

∫

T [ λ ( ∆T ) ∆W - QW + ρ c p W

V

E T musí být tato derivace rovna nule pro nulovou hodnotu ε , tedy

∂T ]dV + ∫ qWdS + ∫ α (T - T ∞ )WdS = 0. ∂t S Sα q

Objemový integrál v souèinu gradientù teploty a funkce W nyní upravíme pouitím Gaussovy vìty (která dí, e objemový integrál divergence vektoru se rovná povrchovému integrálu normálové sloky tohoto vektoru):

∫

V

r r r λ∆ .(W∆T)dV = ∫ Wn. λ∆TdS = ∫ Wn. λ∆TdS + ∫ Wn. λ∆TdS. S

Sq

Sα

Souèasnì vak platí identita (o divergenci souèinu, která se poèítá podobnì jako derivace souèinu funkcí)

∫

V

T λ∆ .(W ∆T) dV = ∫ λ [( ∆W ) ∆T +W ∆ 2 T] dV, V


96

take dosazením (184,185) do rovnice (183) získáme integrální rovnici

- ∫ W( λ ∆ 2 T + Q - ρ c p V

∂T r r )dV + ∫ W(q + λn. ∆T)dS + ∫ W[ α (T - T ∞ )+ λn. ∆T]dS = 0, ∂t S Sα q

která musí platit pro libovolnou funkci W. To ovem znamená, e musí být splnìna Fourierova Kirchhoffova rovnice (179) (ta toti anuluje první, objemový integrál), dále rovnice

r q + λn. ∆T = 0

(která zajistí splnìní okrajové podmínky druhého druhu na hranici

r α (T - T ∞ )+ λn. ∆T = 0.

S q ) a koneènì na hranici S α rovnice

(co je okrajová podmínka tøetího druhu). Ekvivalence variaèního principu (180) s F.K. rovnicí (179) a pøíslunými okrajovými podmínkami je tím dokázána. Poznámka: Z uvedeného odvození je zøejmé, e variace se nevztahují na parciální derivaci teploty dle èasu, take èasovou derivaci teploty je tøeba urèovat jinak ne z variaèního principu. Dále: pokud je T(t,x,y,z) jen aproximací teplotního pole, splòuje pøesnì pouze okrajové podmínky prvního druhu (pøedepsané teploty na hranici), zatímco okrajového podmínky druhého i tøetího druhu stejnì jako vlastní diferenciální rovnice platí jen pøiblinì, co vyplývá z (188). Proto se øíká, e okrajové podmínky prvního druhu jsou silné a okrajové podmínky, v nich figurují tepelné toky, slabé. Koneènìprvková aproximace teplotního pole se volí jako lineární kombinace tvarových funkcí

r

N

r T(t,x, y,z) = ∑ T i (t) N i (x,y,z)= N T T i=1

kde symbolem

r r T jsme oznaèili vektor uzlových teplot (je funkcí èasu) a N je vektor korespondujících

tvarových funkcí (závisí toliko na prostorových souøadnicích x,y,z). Dosazením aproximace (189) do funkcionálu (180) a jeho minimalizací vzhledem k vektoru uzlových teplot (minimalizace se netýká nestacionárního èlenu, viz. pøedchozí poznámka) získáme soustavu obyèejných diferenciálních rovnic pro uzlové teploty (detaily snad není tøeba uvádìt, protoe postup je naprosto stejný jako u Lagrangeova variaèního principu):

r t dT t r r + K λ T = Q, Mc dt

r rT r Tt r r rT r r r t t r M c = ∫ ρ c p N N dV K λ = ∫ ( ∆N ) Λ∆NdV - ∫ αN N dS Q = ∫ QNdV - ∫ qNdS + ∫ α T ∞ NdS V

V

Sα

V

Sq

Sα

Øeení této soustavy obyèejných diferenciálních rovnic je moné získat vpodstatì stejným postupem jako pøi


97

analýze vynuceného kmitání, tj. modální analýzou (pøevedením na soustavu nezávislých diferenciálních rovnic). Øeením tìchto "separovaných" rovnic jsou Duhamelovy integrály nebo lze pouít pøímou integraci. Protoe v tomto pøípadì jde o soustavu diferenciálních rovnic pouze prvního øádu jsou zcela na místì metody Runge Kutta nebo metody Eulerovy. Program COSMOS pouívá právì tuto poslední (a nejjednoduí) monost - explicitní Eulerovu integraci, která má vak tu nevýhodu, e èasový integraèní krok nelze volit libovolnì a musí být dostateènì malý, nemá-li øeení divergovat. COSMOS nenabízí ádné prostøedky pro zjitìní maximální tolerovatelné délky èasového kroku a proto je tøeba se spokojit alespoò s orientaèním odhadem

∆t <

∆ x2 , λ 2 ρ cp

kde ∆x je charakteristický rozmìr pouitých elementù (stabilitní omezení (191) platí pro jednorozmìrné øeení rovnice vedení tepla metodou sítí /viz. pøedmìt základy numerické matematiky/). Pøíklad: Elektrické vyhøívání dna nádoby. Uvaujme extrémnì zjednoduené schema: dno pánve je nerezový

sd = 2 mm , λ = 15W.m-1 .K -1 ) spojený øadou paralelních eber se základovou deskou, její teplota je známá (napø. T 0 = 0 K). Vzdálenost eber je H=0.1 m a jejich výka B=0.05 m. Z jedné strany je dno plech (tlouky

ohøíváno topnými tìlesy (v kadé sekci je mezi ebry jedno tìleso, které do dna pánve vyzaøuje teplo s konstantní hustotou tepelného toku

q = 5000 W.m-2 ) a na druhé stranì dna je kapalina se støední teplotou

T ∞ = 60 , pøièem souèinitel pøestupu tepla stìna-kapalina je α = 1000 W.m-2 .K -1 . Celý dìj budeme povaovat za kvazistacionární a pokusíme se zjistit jaká èást výkonu topného tìlesa se pøenáí do kapaliny a hlavnì jaký bude prùbìh teploty podél dna (co je dùleitá informace pøedevím z hlediska nebezpeèí napékání vsádky na topných plochách pánve). Na základì vypoèteného teplotního pole lze potom stanovit teplotní napìtí ve dnì i v ebrech a i jeho superpozici s napìtím, které odpovídá hydrostatickému tlaku v nádobì. Vzhledem k symetrii problému staèí uvaovat jen jednu polovinu jedné topné sekce dle obrázku 30. Nejprve ukáeme, e variaèní princip (180) lze vyuít k získání pøibliného øeení i bez programù pro MKP. Zanedbáme zmìny teploty napøíè stìnou a budeme uvaovat lineární prùbìh teploty podél ebra i podél dna. Oznaèíme-li T1 teplotu v místì styku ebra a dna, bude teplotní pole popisováno vztahy

T d (x) = T 1 + cx, T _ (y) = T 1 (1-

y ). B


98

V aproximaci (192) jsou pouze dva neurèené parametry: T1, c, jejich "optimální" hodnoty stanovíme minimalizací funkcionálu (180), který se v naem pøípadì redukuje na souèet tìchto integrálù

sd E T (c,T 1 ) = 2

H/2

∫ 0

B

s_ T 12 λ c dx + ∫ λ 2 dy + q 4 0 B 2

H/2

∫

H/2

( T 1 + cx)dx + α

0

∫ 0

1 2 [ ( T 1 + cx ) - ( T 1 + cx)T ∞ ]dx. 2

Stanovením minima (193) (øeením soustavy dvou lineárních rovnic) získáme tyto výrazy pro parametry teplotních profilù

α sd λ + ) λ s_ T 1 48 H 2 . T1 = 2 , c= α sd λ α sd α λ s _ 2 4B ( + ) H ( λ 2 + )( +α )48 H 2 16 H 12 BH ( α T ∞ - q)(

Vyèíslením tìchto výrazù pro nae zadání získáme teplotu dna v místì ebra T 1=63,62 K a nejvyí teplotu uprostøed mezi ebry T=65,67 K (uvádìné hodnoty jsou vztaeny k referenèní nulové teplotì paty ebra). Pro sestavení koneènìprvkového modelu pouijeme elementy PLANE2D; viz text programu: PT,1,0,0,0,Definice vztaných bodù ve vrcholech obrysu PT,2,.001,0,0, PT,3,.05,0,0, PT,4,.05,.002,0, PT,5,.001,.002,0, PT,6,.001,.05,0, PT,7,0,.05,0, PT,8,0,.002,0, PTTOL,.000001, Zmenení limitní vzdálenosti bodù (body s mení odlehlostí se "slijí") SF4PT,1,1,2,5,8,0, Vytvoøení 3 ploch (obdélníkù): rozhraní ebro - dno SF4PT,2,2,3,4,5,0, dno SF4PT,3,8,5,6,7,0, ebro EGROUP,1,PLANE2D,0,1,0,0,0,0,0, volba typu elementù RCONST,1,1,1,1,1, "tlouka" elementù (uvaujeme jednotkovou délku topné sekce) MPROP,1,KX,15, tepelná vodivost (nerez) MPROP,1,DENS,7700, hustota MPROP,1,C,450, mìrná tepelná kapacita MPROP,1,EX,2E11, modul prunosti MPROP,1,NUXY,.3, Poissonova konstanta MPROP,1,ALPX,1.3E-5, souèinitel teplotní roztanosti M_SF,1,3,1,4,3,3,1,1, generování sítì 3 x 3 ètyøuzlových elementù na kadé ploe NMERGE,1,108,1,0.0001,0,1,0, slouèení uzlù na rozhraní ploch NCOMPRESS,1,108,1, pøeèíslování uzlù CECR,1,1000,60,5,4,0, souèinitel pøestupu tepla (1000) a teplota mladiny (60) na køivce 1 a 5 NTCR,8,0,8,1, nulové teploty na køivce 8 (základna ebra) HXCR,6,5000,6,1, hustota tepelného toku 5000 na køivce 6 DCR,8,ALL,0,8,1, "vetknutí" ebra na základnu DCR,9,UX,0,9,1,, nulové posuvy ve smìru x podél køivek 1,7,9 DCR,1,UX,0,1,1,, DCR,7,UX,0,7,1,, PCR,5,1E4,5,1,10000, Dno pánve zatíeno hydrostatickým tlakem p=0.1 bar A_THERMAL,S,0.001,1,1,20, Poadavky na teplotní analýzu (statika) R_THERMAL Výpoèet teplotního pole A_STATIC,T,0,0,1e-06,1e+10,0,0,0, Poadavky na statický výpoèet (T-uvaování teploty) R_STATIC

Statický výpoèet

Výsledky výpoètu teplotního pole i statiky jsou uvedeny na následujícím obrázku 31


99

Z hodnot zobrazených izoterem je zøejmá pozoruhodnì dobrá shoda se jednoduchouèkým analytickým øeením (194). V okénkách èíslo 2 a 3 jsou zobrazeny prùbìhy intenzit napìtí pro zatíení hydrostatickým tlakem (maximální napìtí cca 5,4 MPa) a napìtí teplotní; ta jsou pomìrnì velmi malá (cca 1 MPa), ale jen díky tomu, e byla odstranìna okrajová podmínka v ose topné sekce, která zabraòovala volné dilataci dna ve smìrech x i y (co zhruba odpovídá realitì). Pro ilustraci: Kdyby se dilataèní napìtí uplatnit mohla (co lze modelovat zablokováním posuvù ve smìru x na køivce 7), nalezli bychom ve výsledcích i napìtí dosahující 200 MPa (porovnejte s hrubým odhadem

σ = Eα∆T _170MPa ).

Poznamenejme, e maximální napìtí zpùsobené hydrostatickým tlakem lze spoèítat analyticky (øeením

p H2 H 2 Mo ) ) pøípadu oboustranì vetknutého nosníku se spojitým obtíením σ max = = = p( 2 sd W o 24W o σ max = 6,25MPa , co je cca o 15% vyí hodnota ne výsledek MKP (obr. 31). Vzhledem k tomu, e prùøez dna byl nahrazen pouze 3 x 3 elementy PLANE2D, které nejsou urèeny k modelování deformací nosníkù, zdá se být dosaená pøesnost vypoètených napìtí pøijatelná.

Co øíci k tomuto pøíkladu závìrem: asi tolik, e vechny základní informace o rozloení teplot i o velikosti napìtí bylo moné urèit i bez metody koneèných prvkù, pøedevím díky relativnì jednoduché geometrii a protoe lo o stacionární pøípad bez komplikací zpùsobených sáláním èi promìnností termofyzikálních parametrù. A právì komplikacím tohoto druhu bude vìnován následující pøíklad (fázové zmìny, Stefanùv problém): Chladnutí rotaènì symetrického hliníkového odlitku, jeho geometrie je znázornìna na obrázku 32 (osy symetrie y i x). Pøedpokládáme, e na poèátku (v èase nula) je forma zaplnìna tekutým hliníkem o teplotì

T = 700 o C a v prùbìhu chlazení je udrována konstantní teplota povrchu T 0 = 500 o C (realistiètìjí


100

okrajovou podmínku tøetího druhu bychom mohli modelovat stejnì jako v pøedchozím pøíkladì). Fázová zmìna, k ní pøi chladnutí bude docházet, se dá ve Fourierovì Kirchhoffovì rovnici (179) postihnout uvaováním teplotní závislost mìrného tepla, asi tak jak je uvedeno na obr. 32. Mìrné teplo je zde pro obì fáze (kapalný i tuhý hliník) povaováno za konstantní

-1

c p = 1000J. kg .K -1 a jenom v malém okolí

∆T teploty fázového pøechodu T f se objeví "puls", jeho plocha je rovna skupenskému teplu tuhnutí -1

l = 401kJ. kg . Potom toti bude splnìna podmínka, e skupenské teplo odpovídá zmìnì mìrné entalpie pøi o fázovém pøechodu (uvìdomte si, e platí h = ∫ c p TdT ). Zvolíme-li tedy íøku pulsu ∆T = 5 C , vyjde jeho l -1 = 80200 J. kg .K -1 . výka c p = ∆T Vdy kdy modelujeme nestacionární pøechodový dìj je tøeba nejprve urèit dobu jeho trvání. V pøípadì

dx f [m / s] , kterou se bude pohybovat rozhraní tuhé a kapalné dt fáze. Uvaujeme-li toti jednorozmìrný pøípad (tuhnutí kapaliny o teplotì T f , která zaplòuje poloprostor x>0) tuhnutí odlitku lze napø. odhadnout rychlost

a pøijmeme-li pøedpoklad, e teplotní profil v tuhé fázi je lineární, platí (dokate):

q=

d xf T f - To cp λ = ρ [l + ( T f - T o )] , 2 dt xf

x f (t) =

2( T f - T o )λt . cp ρ [l + ( T f - T o )] 2

Poznamenejme, e ke stejnému výsledku lze dojít mnohem komplikovanìjím zpùsobem, pøibliným øeením Stefanova problému pouitím erfc-funkcí pro popis teplotního profilu, viz. Heømanský [38]. Pro ná konkrétní pøípad ( ρ = 2380kg

/ m3 , λ = 103W.m-1 .K -1 ) bude èasová závislost odlehlosti rozhraní fází od povrchu rovna x f (t) ≈ 0.0054 t . Vzhledem k charakteristických rozmìrùm odlitku (vnìjí polomìr

0.2 m a výka 0.16 m) mùe být pøimìný èasový krok pro monitorování pohybu rozhraní 10 s a celkový èas 100 s (pro tento èas je odhadnutá poloha rozhraní ve vzdálenosti cca 5,5 cm od povrchu. Text programu, pouívajícího opìt elementy PLANE2D (ale ve variantì rotaènì symetrických prstencù), je následující: CURDEF,TEMP,1,1,500,1000,654.5,1000,655.5,80200,659.5,80200,660.5,1000& ,700,1000, teplotní závislost mìrného tepla (funkce èíslo 1) ACTSET,TP,1, aktivace této funkce MPROP,1,C,1, mìrné teplo = 1 krát aktivovaná teplotní závislost ACTSET,TP,0, aktivace funkce èíslo 0 (tj. ádné funkce) MPROP,1,DENS,2380, hustota MPROP,1,KX,103, tepelná vodivost CRSPOLY,1,0,0,0,L,0.2,0,0,L,0.2,0.08,0,A,0.18,0.1,0,L,0.16,0.1,0,A,0.1& 1,0.08,0,A,0,0.04,0,L,0,0,0, obrys meridiánu CT,1,0,.03,1,1, volba obrysu a støední velikosti elementù 0.03 RG,1,1,1, definice oblasti vymezené obrysem èíslo 1 MA_RG,1,1,1,0, generování trojúhelníkù v oblasti 1 MARGCH,1,1,1,Q,4,1,0.4, zmìna trojúhelníkù na ètyøúhelníky EGROUP,1,PLANE2D,0,1,1,0,0,0,0, volba typu elementù (1-rotaènì symetrické prstence PLANE2D)


101

NTCR,2,500,6,1, teploty na køivkách 2 a 6 AUTOSTEP,600,10, automatická volba èasového kroku, ukládání výsledkù vdy po 10-ti sekundách TIMES,0,100,1, poèáteèní a koneèný èas, základní velikost integraèního kroku (1 sekunda) TUNIF,700, poèáteèní teploty stejné ve vech uzlech (700 stupòù) NL_SOLN,1,0, metoda øeení (MNR)

A_THERMAL,T,0.001,1,1,20, R_THERMAL

poadavek na provedení nestacionární analýzy, max. 20 iterací v kadém èasovém kroku sputìní výpoètu

Pozornost si zasluhuje pøedevím pøíkaz AUTOSTEP,Tref,dt_print, který byl navren speciálnì pro pøípady vyetøování fázových zmìn. Tref je referenèní teplota slouící k výpoètu optimálního èasového kroku; tato vìta samozøejmì mnoho informací nesdìluje, ale v manuálech [33] nic bliího uvedeno není. Výsledky výpoètu jsou ve formì izoterm uvedeny pro tøi hodnoty èasu na obrázku 32. Obrázek 32 byl vytvoøen pøíkazy WCREATE,4, vytvoøení 4 oken ACTTEMP,1,TEMP teploty v èasovém kroku èíslo 1 (10 sekund) SETPLOT,3,650,570,0 3-rozliovací schopnost, 750-770 - interval teplot TEMPPLOT,1; vykreslení izoterm èarami (1) aby mohly být vytitìny na èernobílé tiskárnì

COSMOS v procesním inženýrství a strojnictví

Recommend Documents