Obr. 1: Tažený prut
Tažený prut z prostého betonu:
Nelineární dynamický systém s uvážením plasticity
a reologických jevů materiálu
Petr Frantík
1 Úvod
Tato práce předkládá modelování vlivu plasticity a reologických vlastností betonu na jednoduché konstrukci za účelem zjištění možností nástrojů teorie dynamických systémů pro simulaci těchto jevů. Použité náhrady vlastností betonu měly být samostatně vytvořeny a srovnány s náhradami použitými v běžné praxi. Tyto požadavky se podařilo splnit, ačkoliv byly zjištěny významné potíže při vytváření náhrady pracovního diagramu betonu. Význam těchto potíží, jak bude dále uvedeno, byl potlačen neuvažováním změny znaménka napětí v betonu (pro jednoduchost bylo zvoleno tahové napětí).
2.1 Konstrukce
Pokus modelovat společně všechny tyto vlastnosti dynamickým systémem vyžadoval zvolení jednoduché konstrukce. Jinými slovy konstrukce, která se dá bez podstatných ztrát zjednodušit na model s jedním stupněm volnosti a kde napětí na průřezu je konstantní. Takovouto konstrukcí je normálově namáhaný prut. Jeho schéma je na obr. 1. Tuto konstrukci lze zjednodušit na prut, kde normálové napětí rovnoběžné s osou prutu na takto namáhané konstrukci je konstantní po délce i po průřezu.
Obr. 1: Tažený prut
2.2 Model
Konstrukci na obr. 1 nahradíme modelem s jedním stupněm volnosti, jenž je schematicky znázorněn na obr. 2. Na obrázku 2 jsou zobrazeny tři parametry, z nichž každý odpovídá modelovému zjednodušení některé vlastnosti. Parametr m bude náhradní hmotnost hmotného bodu odpovídající kmitající hmotě. Zvolíme jej pro jednoduchost tak, že bude roven hmotnosti celého prutu. Parametr k je tuhost pružiny odpovídající tuhosti prutu v tahu. Jelikož budeme uvažovat nelineární plastické chování betonu, nebude tato tuhost konstantní. Parametr c je konstanta útlumu, který bude zajišťovat ztrátu energie vlivem jejího rozptýlení způsobeného třením, odporem vzduchu apod. Útlum budeme uvažovat lineární. Dále si rozebereme sestavení jednotlivých vztahů, které potřebujeme pro vytvoření dynamického systému.
Obr. 2: Dynamický model
Rovnice kmitání
Pro sestavení rovnic kmitání hmotný bod uvolníme a rozepíšeme si vše, co na něj působí. Všechny síly, které na hmotný bod působí jsou znázorněny na obrázku 3.
Obr. 3: Uvolněný hmotný bod
Síla Fk odpovídá působení pružiny, síla Fs odpovídá setrvačné síle, síla Ft tlumící síle, Fg je tíhová síla a Fb je budící síla odpovídající zatížení. Tyto síly lze rozepsat:
Kde A je plocha průřezu, s je normálové napětí rovnoběžné s osou prutu, x je výchylka hmotného bodu ve směru osy prutu, l je délka prutu, a je zrychlení hmotného bodu, m je hmotnost hmotného bodu, v je rychlost hmotného bodu, c je konstanta útlumu, g je tíhové zrychlení, F(t) je zatížení.
Pro takto uvolněný bod lze potom psát podmínku dynamické rovnováhy:
Protože soustavu diferenciálních rovnic budeme řešit numericky, potřebujeme upravit tuto diferenciální rovnici druhého řádu na soustavu dif. rovnic prvního řádu. Výsledné diferenciální rovnice kmitání tedy budou:
Pracovní diagram
Chování betonu pod zatížením je složité. Deformace roste s napětím nelineárně a poklesne-li napětí, pak deformace neklesá stejným způsobem jako stoupala (viz obr. 4: červená je zvyšování napětí, zelená je snižování napětí). Vznikají trvalé deformace. Takové chování popisujeme jako plasticitu.
Obr. 4: Schéma pracovního diagramu betonu
Pro tuto práci byl pracovní diagram zjednodušen, protože definice funkcí zajišťujících kompletní popis diagramu je složitá. Zjednodušení je takové, že nelineární zůstává pouze ohraničující funkce, která je funkcí, jež by byla správná, kdyby napětí neklesalo. V případě, že napětí poklesne, pak klesá i následně stoupá lineárně až po průsečík s ohraničující funkcí. Takový pracovní diagram je znázorněn na obrázku 5 (červená je zvyšování napětí na ohraničující funkci, modrá je snižování a zvyšování napětí pod ohraničující funkcí).
Obr. 5: Zjednodušení pracovního diagramu betonu
Pro popis ohraničující funkce byl vybrán neúplný polynom třetího stupně s lineárním a kubickým členem:
kde Ep je počáteční modul pružnosti a a je násobitel Ep pro kubický člen. Počáteční modul pružnosti byl zvolen 40 GPa a násobitel a byl zvolen -2.107. Lineární část diagramu je popsána návratovým modulem pružnosti Er. Návratový modul pružnosti byl zvolen 45 GPa.
Již v úvodu bylo zmíněno, že vznikly potíže se zjednodušeným pracovním diagramem při změně znaménka. Zjednodušení, které je použito totiž způsobí, že při dosažení nulového napětí na návratové lineární části dojde při myšleném zvýšení napětí k nejednoznačné situaci. Kdyby bylo připuštěno libovolně malé opačné napětí, musela by být ohraničující funkce posunuta z počátku, což by způsobylo onu nejednoznačnost. Tyto potíže by byly odstraněny s odstraněním zjednodušení.
Reologické jevy
Další komplikací pro výpočet betonových konstrukcí jsou reologické vlastnosti betonu. Pro výpočet jsou podstatné dvě. Smršťování betonu a dotvarování betonu. Smršťování betonu je způsobeno odpařováním vody z betonu. Poločas smrštění je přibližně 28 dní (za 28 dní dojde k polovičnímu smrštění) Velikost smrštění závisí na druhu betonu a na prostředí, v němž je beton uložen. Dotvarování betonu je způsobeno transportem vody v betonu vlivem zatěžování a jejím odpaření. Je to tedy podobný jev jako smršťování, ale je závislý na prostředí i na napětí.
Pro oba tyto jevy budeme uvažovat stejné funkce ve tvaru:
Poločas smrštění a dotvarování byl zvolen 28 dní. Tedy:
kde t28 je 28 dní v sekundách.
Konečné smrštění bylo zvoleno as = 1mm/m. Konečné dotvarování bylo zvoleno jako lineární funkce napětí:
kde parametr R je zvolen tak, aby konečné dovarování při napětí 3 MPa bylo 0.9mm/m. Tedy R = 3.10-10 Pa-1.
Pro délku nenapjatého prutu lze psát:
Kde L je počáteční délka nenapjatého prutu a bs = bd = b.
Diferenciální rovnice reologických jevů potom bude:
Na závěr si ukážeme srovnání průběhu zvolené funkce s funkcí použitou v normě. Norma používá funkci:
Na obrázku 6 je srovnání obou funkcí, přičemž parametry b = 2.7172 a a = 1. Je vidět, že normová funkce má nekonečnou hodnotu derivace v počátku a její hodnota pro 28 dní rozhodně není 0.5, což by odpovídalo poločasu smršťování.
Obr. 6: Srovnání zvolené funkce a funkce použité v normě
2.3 Simulace
Pro kontrolu bylo třeba provést simulace jednotlivých jevů. Byl vytvořen program v jazyce Pascal, který je uveden v části 5 a byly provedeny jednotlivé simulace.
Pracovní diagram
Obr. 7: Vypočtený tahový pracovní diagram betonu
Smršťování
Obr. 8: Vypočtený průběh smršťování
Dotvarování
Obr. 9: Vypočtený průběh dotvarování pro napětí 3 MPa
Zatěžovací pokus
Pro simulaci zatěžovacího pokusu byl vybrán model konstrukce s těmito parametry:
Zatěžovací plán byl použit takový, že se do bednění vybetonuje prut, po 10 dnech se odbední a po 28 dnech se na 4 dny zatíží silou 25 kN.
Obr. 10: Vypočtené napětí
Obr. 11: Vypočtené protažení prutu vlivem napětí
Obr. 12: Změna délky prutu vlivem reologických jevů
3 Závěr
Všechny požadované jevy se podařilo úspěšně simulovat pomocí diferenciálních rovnic dynamického systému.
Důležitým poznatkem byl problém s diskretizačním krokem numerických metod. Pro kmitání systému byl vyžadován krok maximálně jedna statitisícina sekundy, což pro požadovanou délku simulace bylo příliž málo (simulace by na AMD K6/300 MHz trvala příliž dlouho). Tento problém byl vyřešen odstavením řešení diferenciálních rovnic kmitání vždy tehdy, když byl systém ustálen a krok byl zvýšen na 10 sekund. Automaticky při změně rozhodujícího parametru (zatěžovací síla) byl pak krok snížen a diferenciální rovnice kmitání se započaly řešit.
Objevily se potíže při vytváření náhrady pracovního diagramu betonu, kterou by bylo vhodné vylepšit. Nesoulad mezi reologickými funkcemi neumím vysvětlit.
4 Literatura
5 Zdrojový kód programu
{$N+}
Uses Crt,Graph;
Const
{Metoda 0-RK, 1-EE, 2-SE, 3-ME, 4-O, 5-MO}
meti=0;
{Krok metod}
hm1=0.00001; {dynamicky 0.00001}
hm2=10; {staticky 10}
{}
dgx=0.00001; {0.00001}
cxgi=5; {5}
dgy=0.01; {0.01}
cygi=2; {2}
{}
{parametry}
Ep=40E+9; {40E+9 Pa}
alfa=-2E+7; {-2E+7}
Er=45E+9; {45E+9 Pa}
A=0.01; {0.01 m2}
m=150; {150 kg}
c=500000; {10000}
{}
Reo=3E-10; {ad=Reo*sig; max sig=30E+6: ad=0.009}
{}
time1=60*60*24; {v sekundach 1 den}
time10=time1*10; {v sekundach 10 dni}
time28=time1*28; {v sekundach 28 dni}
time32=time1*32; {v sekundach 32 dni}
{}
{auto deleni os}
osb=false;
Var
g:double; {9.81 m/s2}
{}
{}
grdi,gmi,erri,ki:integer;
{}
h,h2,hd2,hd6:double;
{Stavove promenne}
x,y,l:double;
x0,y0,l0:double;
xp,yp,lp:double;
F:double;
{}
{Casove promenne}
time:double;
{}
{parametry}
{reologie}
ad,bd,as,bs:double;
{}
{cleny funkce tahoveho diagramu}
sig:double;
epsm,epsp,epshr:double;
{}
kx,ky,px,py:real;
{}
il:longint;
{}
smodei:integer; {prepinac casoveho kroku 1-dynamic,2-static}
{}
fo:file of double;
zari:integer;
{Rozsirujici funkce}
{$I Tools}
{$I Input}
{****************************************************************************}
{Soustava diferencialnich rovnic systemu}
Function F1(y1,y2,y3:double):double;
begin
if smodei=1 then
begin
F1:=y2;
end
else
begin
F1:=0;
end;
end;
{******}
Function Fung(eps:double):double;
function ghr(x:double):double;
begin
ghr:=Ep*x*(1+alfa*x*x);
end;
function glin(x:double):double;
begin
glin:=Er*x;
end;
begin
if eps>=epshr then
begin
sig:=ghr(eps-epsp);
epshr:=eps;
epsm:=epshr-sig/Er;
end
else if epshr>eps then
begin
sig:=glin(eps-epsm);
end;
{ustaveni smrstovaci meze}
ad:=Reo*sig;
{}
fung:=sig;
end;
{******}
Function F2(y1,y2,y3:double):double;
begin
if smodei=1 then
begin
{F2:=-c/m*y2-A/m*fung(y1/y3)+g+F/m;}
F2:=-c/m*y2-A/m*fung(y1/l0)+g+F/m;
end
else
begin
F2:=0;
end;
end;
Function F3(y1,y2,y3:double):double;
begin
F3:=(ad/Moc(bd,time)*ln(bd)-as/Moc(bs,time)*ln(bs))*l0;
end;
{$I Library}
{Konec soustavy}
{****************************************************************************}
{Zobrazeni aktualnich informaci}
Procedure Info;
begin
Bar(1,1,638,30);
SetColor(white);
{stavove}
OutTextXY(10,10,'t: '+Fix(time,4,12)+'('+Num(smodei)+')'+
' l: '+Fix(l,4,6)+' F: '+Fix(F,0,5)+' sig: '+Fix(sig,0,6)+'/'+Fix(F/A,0,6));
OutTextXY(10,20,'x: '+Fix(x,8,12));
{pocatecni podminky}
{OutTextXY(558,10,'x0'+Fix(x0,0,1)+' y0'+Fix(y0,0,1)+' l0'+Fix(l0,0,2));}
{parametry}
OutTextXY(150,20,'Ep'+Fix(Ep/1E+9,0,2)+' alfa'+Fix(alfa/1E+6,0,2)+'E+6 Er'+Fix(Er/1E+9,0,2)+'E+9 A'+Fix(A,2,4)+
' g'+Fix(g,1,3)+' m'+Fix(m,1,3)+' c'+Fix(c,1,5)+' Reo'+Fix(Reo/1E-12,0,2)+'E-12');
end;
Procedure Ustav;
begin
if smodei=1 then h:=hm1 else h:=hm2;
h2:=h*2;
hd2:=h/2;
hd6:=h/6;
end;
{****************************************************************************}
{Zakladni funkce}
{Vynulovani casu systemu}
Procedure Restart;
begin
time:=0;
{Pocatecni podminky}
x:=x0;
y:=y0;
l:=l0;
xp:=x0;
yp:=y0;
lp:=l0;
{}
epsp:=0;
epsm:=0;
epshr:=0;
sig:=0;
end;
{Vypocet koeficientu snizujicich pocet parametru}
Procedure Substit;
begin
end;
{Inicializace promennych}
Procedure Init(kx0,ky0,px0,py0:real);
begin
Restart;
Substit;
{Zobrazovaci promenne}
kx:=kx0;
ky:=ky0;
px:=px0;
py:=py0;
end;
{Klavesova nabidka}
Function Menu(ki:integer):integer;
begin
if (ki>0) and (ki=/27) and (ki=/32) then
begin
if ki=8 then
begin
{}
end
else if ki=9 then
begin
Info;
ki:=290;
end
else if ki=13 then
begin
smodei:=3-smodei;
end
else if (ki=ord('i')) or (ki=ord('I')) then
begin
Info;
end
else if (ki=ord('p')) or (ki=ord('P')) then
begin
F:=F+1000;
Substit;
end
else if (ki=ord('o')) or (ki=ord('O')) then
begin
F:=F-1000;
Substit;
end
else
begin
if ki=1 then
begin
ky:=ky*1.1;
py:=py*1.1;
end
else if ki=2 then
begin
ky:=ky/1.1;
py:=py/1.1;
end
else if ki=3 then
begin
kx:=kx/1.1;
px:=px/1.1;
end
else if ki=4 then
begin
kx:=kx*1.1;
px:=px*1.1;
end
else if (ki=ord('a')) or (ki=ord('A')) then
begin
px:=px-10;
end
else if (ki=ord('d')) or (ki=ord('D')) then
begin
px:=px+10;
end
else if (ki=ord('w')) or (ki=ord('W')) then
begin
py:=py-10;
end
else if (ki=ord('s')) or (ki=ord('S')) then
begin
py:=py+10;
end;
ki:=290;
end;
end;
Menu:=ki;
end;
{Konec zakladnich funkci}
{****************************************************************************}
procedure Save;
var hd:double;
begin
if time=0 then
begin
Assign(fo,'G:\Out.num');
Rewrite(fo);
zari:=0;
end;
if zari=0 then
begin
hd:=time/time1;
Write(fo,hd);
Write(fo,sig);
Write(fo,x);
Write(fo,l);
zari:=1000;
end
else zari:=zari-1;
end;
{****************************************************************************}
Begin
{mod kroku}
smodei:=2; {1}
{Pocatecni podminky}
x0:=0; {musi byt nula (napjatost)}
y0:=0;
l0:=10;
{}
F:=0;
g:=0;
{Parametry reologie}
ad:=0; {konecne dotvarovani}
bd:=moc(2,1/time28); {rychlost dotvarovani 50% za 28 dni}
{}
as:=0.001; {0.001 m/m konecne smrsteni}
bs:=moc(2,1/time28); {rychlost smrstovani 50% za 28 dni}
{}
{textovy uvod}
ClrScr;
TextColor(green);
WriteLn('Program TahPrut');
TextColor(lightgray);
WriteLn('Autor '+Autor+', rok emise 2001.');
WriteLn('(Program pro simulaci disipativniho dynamickeho systemu)');
WriteLn;
Writeln(' Esc - ukonceni');
WriteLn(' Enter - ulozeni');
WriteLn(' Sipky - zoom');
WriteLn(' A,W,S,D - posun os');
WriteLn(' I - informacni okenko');
WriteLn(' P - dana zmena daneho parametru');
WriteLn(' Mezera - restart simulace');
WriteLn('BackSpace - vynulovani znacky lokalni extremni vychylky');
WriteLn(' Tab - ulozeni a vymazani obrazovky');
WriteLn;
TextColor(white);
WriteLn('(Stiskni klavesu)');
WriteLn;
TextColor(lightgray);
if Get=/27 then
begin
{Inicializace grafiky}
grdi:=Detect;
InitGraph(grdi,gmi,'');
erri:=GraphResult;
if erri=GrOK then
begin
{Pauza pro grafiku}
Delay(1000);
SetFillStyle(1,7);
{}
Ustav;
Init(Moc(1.1,150),Moc(1.1,80),0,0);
{cyklus pripadu}
repeat
Restart;
gCls;
MoveTo(Round(320+px),240);
{cyklus simulace}
repeat
ShowTraj(x,y,xp,yp);
{}
Save;
{}
step(meti,x,y,l,xp,yp,lp);
{}
{Zatezovaci plan}
if (time>=time10) and (time10+hm2>=time) then
begin
g:=9.81;
smodei:=1;
end
else if (time>=time28) and (time28+hm2>=) then
begin
F:=25000;
smodei:=1;
end
else if (time>=time32) and (time32+hm2>=time) then
begin
F:=0;
smodei:=1;
end;
{}
ki:=Menu(key);
Ustav;
if ki=290 then gCls;
until (ki=32) or (ki=27);
{konec cyklu simulace}
until ki=27;
{konec cyklu pripadu}
{Uzavreni grafiky}
CloseGraph;
Delay(1000);
end
else
begin
WriteLn('Chyba grafiky: ',GraphErrorMsg(erri));
end;
end;
ClrScr;
end.