Obsah

Nelineární statika

Nelineární statika

Schéma výpočtu

Postup výpočtu

1. Příprava výpočtu

Řešení této úlohy musí předcházet řešení lineární elastostatické úlohy, přičemž platí následující podmínky:

Použít se mohou jen isoparametrické prvky.
Nelze využít obecnou plochu symetrie nebo periodicitu.
Materiálové vlastnosti se zadají podle kapitoly Materiálové vlastnosti a přílohy Nelineární materiál v Referenční příručce.
Zatěžovací stavy se zadají podle odstavce Posloupnost zatěžovacích stavů.

Program: RMD2, RPD2, SRH2, FEFS (2D úloha) nebo RMD3, RPD3, SRH3, FEFS (3D úloha)
Vstupy: name.i1, name.i2, name.i3, name.i4
Protokol: name.o1, name.o2, name.o3, name.o4
Výstupy: binární soubory
Detaily: Lineární elastostatika, Organizace výpočtu, Referenční příručka: Vstupy, Nelineární materiál

2. Zpracování historie zatížení a výpočtového modelu

Vstupní data se zapíší do textového souboru name.iP. Při definici historie zatížení se postupuje podle odstavce Posloupnost zatěžovacích stavů. Pokud již výpočet proběhl (tj. nelineární úloha byla úspěšně vyřešena), lze připojit další zatížení s klíčem $\mathtt{KREST}=2$ (viz odstavec Posloupnost zatěžovacích stavů).

Program: HPP2 (2D úloha) nebo HPP3 (3D úloha)
Vstupy: name.iP, binární soubory z předchozího výpočtu
Protokol: name.oP
Výstupy: binární soubory
Detaily: Organizace výpočtu, Posloupnost zatěžovacích stavů, Specifikace výpočtového modelu, Referenční příručka: Vstupy

3. Řešení soustavy rovnic

Vstupní data se zapíší do textového souboru name.iL. Doporučeným algoritmem řešení je metoda BFGS ($\mathtt{KMET}=2$). Pokud je řešení předčasně ukončeno v důsledku překročení maximálního počtu iterací $\mathtt{NITER},$ je nutné program znovu spustit. Poslední iterace přerušeného řešení je automaticky použita jako výchozí aproximace pro nový běh programu. Před opětovným spuštěním programu je možné (ale nikoliv nutné) změnit vstupní data v souboru name.iL.

Program: HPLS (2D i 3D úloha)
Vstupy: name.iL, binární soubory z předchozího výpočtu
Protokol: name.oL
Výstupy: binární soubory (řešení je v souboru name.PLS)
Detaily: Organizace výpočtu, Řídící parametry řešiče, Referenční příručka: Vstupy

4. Výpočet deformací a napětí

Vstupní data se zapíší do textového souboru name.i5, přičemž klíč typu úlohy $\mathtt{KPROB}=2.$ Vždy se musí zadat $\mathtt{ILC}$ s ohledem na skutečný počet vyřešených stavů. V souborech name.STR a name.STB je efektivní plastická deformace označena jako SCALAR1 a efektivní creepová deformace jako SCALAR2.

Program: STR2 (2D úloha) nebo STR3 (3D úloha)
Vstupy: name.i5, binární soubory z předchozího výpočtu
Protokol: name.o5
Výstupy: name.STR (volitelně), name.STB (volitelně)
Detaily: Organizace výpočtu, Referenční příručka: Vstupy

Posloupnost zatěžovacích stavů

Odezva nelineárních materiálů závisí nejen na velikosti a směru působících sil, ale také na pořadí, v jakém jsou zatěžovací účinky přikládány. Proto je třeba zadat časový sled zatížení. Vycházíme z obvyklého pojmu zatěžovací stav, kterým máme na mysli teplotu a veškeré síly působící na těleso v daném časovém okamžiku. Definice všech zatěžovacích stavů, jež budou využity v nelineárním výpočtu, se provádí běžným způsobem programem RPD2/RPD3. Z takto předem připravených zátěží se potom pomocí programu HPP2/HPP3 vybere a sestaví posloupnost potřebná pro modelování historie zatížení, jak je schematicky znázorněno na následujícím obrázku.

Automaticky se přitom předpokládá, že přechod z $L_i$ do $L_{i+1}$ probíhá rovnoměrně z hlediska deformací. Tento přechod mezi jednotlivými stavy je sice velmi blízký rovnoměrné změně zatížení, nemusí být však vždy totožný. Proto v případech, kdy je třeba bezpodmínečně dodržet silovou zatěžovací cestu, je nutné příslušný úsek rozdělit na menší přírůstky. Zatěžovacím stavem, který je předepsaný, soustava prochází za všech okolností.

K dosažení jemnějšího dělení není nutné zadávat další zatěžovací stavy; jednoduchou možnost poskytuje program HPLS. Ve většině případů je však zbytečné sahat k podobným opatřením, protože zpomalují výpočet (zpravidla zbytečně). Přesnost integrace konstitutivních rovnic je zabezpečena jiným způsobem a v případě creepu těleso prochází rovnovážným stavem na konci každého časového přírůstku.

Zatěžovací stavy se definují následujícím způsobem:

1)

V souboru name.i2 se v rámci dávky AS 1 (tj. v prvním zatěžovacím stavu) běžným způsobem definuje materiál, nulová posunutí uzlů, ev. pružiny, výchozí teplota tělesa a páry kontaktních ploch. Dále uvedená zatížení (tlak, osamělé síly, atp.) se berou v úvahu jen v lineární úloze. Pro nelineárních úlohy se zatěžovací účinky popsané v AS 1 nedají využít a je nutné je zadat v AS 2 a výše. V přiřazení /R 0 $T_o$ $T_w$/ má význam jen počáteční teplota $T_o,$ odpovídající stavu bez napětí. Hodnota $T_w$ se neuplatní.

2)

V dávkách AS 2, AS 3, … se postupně vytvoří všechny zatěžovací stavy, které budou použity v nelineárním výpočtu. Jejich pořadí nehraje roli. Zatížení se nezadává přírůstkově, ale absolutně (vzhledem k nule). Jestliže se například vyskytne přiřazení /R 0 $T_o$ $T_w$/, vzniká zatěžovací stav, ve kterém má těleso teplotu $T_w.$ Hodnota $T_o$ se ignoruje.

3)

Vstupní data se zpracují programem RPD2/RPD3 a dále se provede výpočet matic tuhosti programem SRH2/SRH3 a eliminace soustavy programem FEFS. Možné je též vypočítat elastická napětí programem STR2/STR3.

4)

Sestaví se posloupnost zatěžovacích stavů $L_1,$ $L_2,$ …, kde čísla $L_i$ jsou pořadová čísla AS dávek. Např. $L_2=5$ znamená, že druhý zatěžovací stav byl definován v AS 5 (aktuální přírůstek oproti minulé konfiguraci je $L_2-L_1$). Pokud se dodatečně ukáže, že některý zatěžovací stav chybí, je nezbytné celou úlohu přepočítat od bodu 2.

5)

V souboru name.iP se na IP řádku zadá:

IP KREST NLC NCYC KMOD KCRP KLARG KCNT KTPR KURHS 0 $L_1$ $L_2$ $\dots$ $L_\mathtt{NLC}$

kde

$\mathtt{NLC}$ je počet členů zatěžovací posloupnosti,
$\mathtt{NCYC}$ je počet cyklů (opakování celé posloupnosti, výchozí hodnota je $1$),
$\mathtt{KMOD},$ $\mathtt{KCRP},$ $\mathtt{KLARG},$ $\mathtt{KCNT}$ a $\mathtt{KURHS}$ viz Specifikace výpočtového modelu.

Celkový počet řešených stavů je $\mathtt{NLC}\cdot\mathtt{NCYC}.$

6)

Jestliže úloha nezávisí fyzikálně na čase (např. elastoplasticita), nezadávají se žádné další údaje a $\mathtt{KCRP}=0.$ Pro creepovou úlohu je $\mathtt{KCRP}\ge1$ a v souboru name.iP se na RP řádku zadají časy v hodinách odpovídající všem zatěžovacím stavům:

RP 10*0 $t_1$ $t_2$ $\dots$ $t_\mathtt{NLC}$

Předpokládá se $t_\mathtt{NLC}>\ldots>t_2>t_1>0$ a $\mathtt{NCYC}$ je vždy $1$ (tj. je možné zapsat i $\mathtt{NCYC}=0$).
Creepový výpočet se automaticky kombinuje s elastoplasticitou. Vyloučení plastických konstitutivních vztahů je možné dosáhnout zadáním dostatečně vysoké meze kluzu v souboru name.i2.

Ve výpočtech dlouhodobého tečení s počátečním elastickým stavem materiálu se často postupuje tak, že se vytvoří jediný zatěžovací stav $L_e,$ popisující konstantní zátěž, a pak se zadá $\mathtt{NLC}=2.$ posloupnost $L_1=L_e,L_2=L_e$ a časy $t_1=0,t_2=t_\text{end}.$ Znamená to, že se těleso nejprve zatíží v nulovém čase na $L_e$ (a tedy elasticky) a potom probíhá tečení po dobu $t_\text{end}$ (koncový zatěžovací stav $L_e$ je stejný jako na začátku).

7)

Vstupní data se zpracují programem HPP2/HPP3 a úloha se spočítá programem HPLS. Pokud řešení proběhlo úspěšně, můžeme navázat v bodě 5 zadáním dodatečných zatěžovacích stavů, které však musely být definovány předem v rámci AS dávek v souboru name.i2 – návrat k bodu 2 znamená přepočítání celé úlohy. V souboru name.iP stačí zapsat na první pozici IP řádku $\mathtt{KREST}=2$:

IP 2 NLC NCYC KMOD $\dots$

U creepových úloh se přirozeně předpokládá, že první čas $t_1$ je větší nebo roven času, ve kterém bylo ukončeno předchozí řešení. V opačném případě program HPP2/HPP3 hlásí chybu. Restartu je možné využít pro řešení úloh s předpětím. Nejprve se vypočítá předpětí s $\mathtt{KREST}=1$ (např. zbytková pnutí po procesu chladnutí) a pak se odstartuje nová série zatížení (např. cyklického) s $\mathtt{KREST}=2.$ Tímto způsobem se dá pracovat s materiálovými vlastnostmi, které jsou změněny předchozí historií.

Příklad

Vyšetříme zbytková pnutí v tělese po ohřevu z $T_1$ na $T_2>T_1$ a následném zatížení silou $F.$

V souboru name.i2 bude:

…
;Definice materiálu, okrajových podmínek a základní teploty.
;Vzhledem k tomuto stavu budou odečítána posunutí a napětí.
AS 1 /… /R 0 <T₁> <T> ; T libovolná
;Síla F při teplotě 0°C.
AS 2 /<přiřazení F>
;Teplota T₂ bez síly.
AS 3 /R 0 <T> <T₂> ; T libovolná
;Teplota T₂ + síla F.
AS 4 /R 0 <T> <T₂> /<přiřazení F> ; T libovolná
;Výchozí stav (nutno zadat, aby bylo možno popsat odlehčení).
AS 5 /R 0 <T> <T₁> ; T libovolná
…

Sestavme nyní několik posloupností v souboru name.iP:
a) chybně:

IP 1 3 0 0 6*0 3 2 5

3 – ohřev z $T_1$ na $T_2$
2 – zatížení $F,$ ale současné ochlazení z $T_2$ na $0~^\circ\text{C}$
5 – odlehčení $F$ a současný ohřev z $0~^\circ\text{C}$ na $T_1$

b) chybně:

IP 1 3 0 0 6*0 4 3 5

4 – ohřev z $T_1$ na $T_2$ a současné zatížení $F$
3 – odlehčení $F$ při konstantní teplotě
5 – ochlazení z $T_2$ na $T_1$

c) správně:

IP 1 3 0 0 6*0 3 4 5

3 – ohřev z $T_1$ na $T_2$
4 – zatížení $F$ při konstantní teplotě $T_2$
5 – ochlazení z $T_2$ na $T_1$ a současné odlehčení $F$

Všimněme si, že ve všech případech je výsledný stav zatížení stejný, jako byl na začátku, ale zbytková pnutí se budou lišit. Zadání příkladu nejlépe odpovídá posloupnost c). Jinou možností by bylo pořadí 3 4 3 5, protože nebylo přesně specifikováno, zda má odlehčování probíhat současně nebo postupně.

Příklad

Tyčka je zatížena střídavým napětím $\sigma_a.$ Předtím materiál prodělal jednorázové zatížení a odlehčení napětím $\sigma_o.$ Vyšetříme průběh prvních dvaceti cyklů.

V souboru name.i2 bude:

…
AS 1 /…
AS 2 /<přiřazení σ = 0>
AS 3 /<přiřazení σ = σₐ>
AS 4 /<přiřazení σ = -σₐ>
AS 5 /<přiřazení σ = σₒ>

První zatížení a odlehčení se popíše v souboru name.iP:

IP 1 2 0 0 6*0 5 2

Po proběhnutí programu HPLS opakujeme řešení s $\mathtt{KREST}=2$:

IP 2 2 20 0 6*0 3 4

Opětovným spuštěním programu HPLS získáme průběh 20 cyklů.

Specifikace výpočtového modelu

Výpočtový model vyplývá z materiálových vlastností popsaných v souboru name.i2, příp. name.DAT, a ze vstupních parametrů v souboru name.iP:

IP KREST NLC NCYC KMOD KCRP KLARG KCNT KTPR KURHS 0 $L_1$ $L_2$ $\dots$ $L_\mathtt{NLC}$ RP 10*0 $t_1$ $t_2$ $\dots$ $t_\mathtt{NLC}$

Elastoplasticita

Model plasticity se volí klíčem $\mathtt{KMOD}$:

$=0$ … elasticita (výchozí)
$=1$ … von Misesův model
$=2$ … zobecněný asociovaný model
$=3$ … zobecněný neasociovaný model
$=4$ … Feigenbaum-Dafaliasův model

Zvolenému modelu musí odpovídat materiálové vlastnosti zadané v souboru name.i2.

Časová posloupnost $t_i$ se nezadává.

Creep

Model creepu se volí klíčem $\mathtt{KCRP}$:

$=0$ … bez creepu (výchozí)
$=1$ … PMD model
$=2xx$ … Bínův model
$=2$ … Nortonův model
$=3$ … Norton-Baileyův model
$=4$ … Time Hardening model
$=5$ … MPC Project Omega model
$=6$ … Klocův model

přičemž klíč $\mathtt{KMOD}=0.$ Při $\mathtt{KCRP}=1$ se v souboru name.i2 zadá závislost $\dot\varepsilon_c=\dot\varepsilon_c(\sigma_e,\varepsilon_c,T).$ Při $\mathtt{KCRP}>1$ se popis creepových vlastností materiálu čte ze souboru name.DAT.

Na RP řádku se přiřadí časy $t_i$ (v hodinách) všem zatěžovacím stavům $L_i,$ viz odstavec Posloupnost zatěžovacích stavů.

Elastoplasticita s creepem

Postupuje se podle předchozích odstavců. Zadá se klíč $\mathtt{KMOD}>0$ rozlišující plastický model, klíč $\mathtt{KCRP}>0$ rozlišující creepový model a zároveň se zadá časová posloupnost $t_i.$

Pokud $t_{i+1}=t_i,$ proběhne změna z $L_i$ do $L_{i+1}$ náhle s vyloučením creepu. Kombinovaný výpočet je vhodné spustit až po odladění separovaných úloh.

Geometricky nelineární úloha

Aktivuje se při $\mathtt{KLARG}>0.$

$\mathtt{KLARG}=0$ … geometricky lineární (výchozí)
$=1$ … totální Lagrangeovská formulace (velká posunutí, malé deformace)
$=2$ … aktualizovaná Lagrangeovská formulace (velká posunutí, malé deformace)
$=3$ … logaritmická formulace (velká posunutí, velké deformace)
$=4$ … korotační formulace
$=5$ … Mooney-Rivlinův hyperelastický model
$=6$ … Ogdenův hyperelastický model

$\mathtt{KURHS}=1$

Kontakt

Aktivuje se při $\mathtt{KCNT}=1.$

V souboru name.i2 je nutné zadat páry kontaktních ploch jako SV dávky s $\mathtt{KQT}=10$ a $11,$ $12$ a $13,$ atd.

V souboru name.iL je nutné nastavit hodnotu penalty $\mathtt{PENAL}.$ Pro styk ocelových materiálů se volí $\mathtt{PENAL}\approx10^{13}.$ Pokud úloha nekonverguje, doporučuje se hodnotu $\mathtt{PENAL}$ snížit nebo zatěžovací krok rozložit parametrem $\mathtt{NSUBI}$ na několik kroků. Pokud výsledky kontaktu nejsou optimální, doporučuje se hodnotu $\mathtt{PENAL}$ zvýšit. Také je možné začít s nižší hodnotou $\mathtt{PENAL}$ a tu postupně zvyšovat.

Tělesa v kontaktu musí být regulárně uložena. Volba tuhosti uložení „volného“ tělesa ve směru kontaktu musí být taková, aby neovlivňovala řešení.

Řídící parametry řešiče

Řídící parametry řešiče HPLS se zadávají v souboru name.iL:

IP KMET KOUT NSUBI NITER NINT KTPR RP UTOL RTOL XTOL PENAL

Význam řídících parametrů záleží na tom, jaký typ nelineární úlohy se řeší. Elastoplastická úloha vede k řešení soustavy nelineárních algebraických rovnic v každém zatěžovacím stavu. Současně probíhá integrace konstitutivních rovnic s deformačně řízeným krokem, na který nemá velikost přírůstku zatížení podstatný vliv. Pro elastoplastickou úlohu jsou tudíž důležité parametry vztahující se k iterační metodě řešení ($\mathtt{KMET},$ $\mathtt{NITER}$ a kritéria konvergence). Pro výpočet creepu je naopak použita časově řízená integrace, která sice nevyžaduje řešení soustavy nelineárních rovnic, klade však důraz na správnou volbu délky kroku. Důležité jsou proto parametry $\mathtt{NSUBI}$ a $\mathtt{NINT}.$

Význam parametrů:

$\mathtt{KMET}$ volí metodu řešení soustavy nelineárních algebraických rovnic. Pro čistě creepovou úlohu nemá tento parametr význam, neboť explicitní časová integrace vede k linearizovaným rovnicím.
- $=1$ … modifikovaná Newton-Raphsonova metoda
- $=2$ … metoda BFGS (výchozí)
$\mathtt{KOUT}$ aktivuje výstup vektorů řešení do binárního souboru name.PLS, kterým pak v programu STR2/STR3 odpovídá pořadové číslo $\mathtt{ILC}$ (viz name.i5). Pokud byla úloha restartována s $\mathtt{KREST}=2$ (viz odstavec Posloupnost zatěžovacích stavů), soubor name.PLS se přepíše, tj. řešení z předchozího běhu se musí zpracovat před restartem. Číslo v závorce udává počet zapsaných vektorů řešení:
- $=0$ … kontrola vstupních dat ($0$)
- $=1$ … po každém zatěžovacím stavu ($\mathtt{NCYC}\cdot\mathtt{NLC}$)
- $=2$ … po každém cyklu ($\mathtt{NCYC}$)
- $=3$ … jen výsledné řešení ($1$)
$\mathtt{NSUBI}$ specifikuje subinkrementaci zatěžovacích stavů (výchzí hodnota je $1$). Dílčí řešení se nezapisují do souboru name.PLS. Dělení má význam především pro creepovou úlohu, kdy je možné pevně předepsat délku časového kroku $\Delta t=(t_{i+1}-t_i)/\mathtt{NSUBI}$ ($t_i$ odpovídají zatěžovacím stavům $L_i,$ viz odstavec Posloupnost zatěžovacích stavů). Zároveň je vhodné nastavit $\mathtt{NINT}=1.$ Pro časově nezávislý problém se parametr $\mathtt{NSUBI}$ uplatní jen vyjímečně, v podstatě jen tehdy, když je třeba vynutit proporcionální silovou změnu mezi $L_i$ a $L_{i+1}$ (viz odstavec Posloupnost zatěžovacích stavů).
$\mathtt{NITER}$ specifikuje maximální počet iterací (výchozí hodnota je $10$).
$\mathtt{NINT}$ specifikuje dělení integračního intervalu (výchozí hodnota je $10$) a souvisí s přesností řešení. Řádová chyba v napětích je podstatně menší než $100/\mathtt{NINT}$ %, Výchozí hodnota dává ve většině případů malou chybu (cca 1–3 %) a nedoporučuje se ji proto měnit. Při řešení creepových úloh je časový krok programem nastavován automaticky. Položení $\mathtt{NINT}=1$ má za následek vyřazení krokovacího algoritmu z činnosti (připouští se 100% lokální chyba). V takovém případě je nutné předepsat časový krok volbou parametru $\mathtt{NSUBI}.$ U elastoplastických úloh může integrace s $\mathtt{NINT}=1$ významně urychlit výpočet, aniž by se příliš poškodilo řešení. Integrační metoda se pak redukuje na algoritmus prediktor-korektor (Euler forward – radial return), vhodné je však omezit velikost silových přírůstků pomocí $\mathtt{NSUBI}.$
Kritéria konvergence $\mathtt{UTOL},$ $\mathtt{RTOL}$ a $\mathtt{XTOL}$ se vztahují k iteračnímu řešiči. Řešení pokračuje tak dlouho, dokud nejsou splněna všechna kriteria současně, nebo se nepřekročí maximální dovolený počet iterací $\mathtt{NITER}.$ Pokud bylo řešení přerušeno v důsledku překročení $\mathtt{NITER},$ je nutné spustit program znovu. Výpočet je automaticky restartován od poslední iterace. Před opětovným spuštěním programu je možné změnit kriteria konvergence, metodu řešení $\mathtt{KMET}$ a integrační parametr $\mathtt{NINT}.$ Index $^{(i)}$ označuje $i$-tou iteraci, index $^{(0)}$ výchozí stav a $\sqrt{\mathtt{LSOL}}$ je odmocnina z celkového počtu neznámých (stupňů volnosti sítě):
- $||\Delta\mathbf{u}^{(i)} || < \mathtt{UTOL}\,||\mathbf{u}^{(i)}||$ (výchozí hodnota je $\mathtt{UTOL}=10^{-3}$)
- $||\mathbf{R}^{(i)}|| < \mathtt{RTOL}\,||\mathbf{R}^{(0)}||$ (výchozí hodnota je $\mathtt{RTOL}=10^{-3}$)
- $\sqrt{\mathtt{LSOL}}\max|\mathbf{R}^{(i)}|<\mathtt{XTOL}\,||\mathbf{R}^{(0)}||$ (výchozí hodnota je $\mathtt{XTOL}=10^{-2}$)
$\mathtt{PENAL}$ se zadává u kontaktních úloh a na jeho správné volbě závisí stabilita i přesnost řešení.

Následující soubor lze doporučit pro většinu úloh:

name.iL

; KMET KOUT NSUBI NITER NINT KTPR
IP  2    1    0     0     0    1
; UTOL RTOL XTOL
RP  0    0    0
EN
EN