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:

  • Tíhové zrychlení g = 9.81 N.m-2
  • Délka prutu l = 10 m
  • Plocha průřezu A = 0.01 m2
  • Hmotnost hmotného bodu m = 150 kg
  • Parametr útlumu c = 5.105
  • N.s.m-1
  • Parametr R = 3.10-10 Pa-1
  • Počáteční modul pružnosti Ep = 40 GPa
  • Násobitel a  = -2.107 GPa
  • Návratový modul pružnosti Er = 45 GPa
  • Konečné smrštění as = 0.001
  • Parametr rychlosti smršťování bs = 2(1/t28)
  • Parametr rychlosti dotvarování bd = 2(1/t28)
  • 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

    1. RNDr. Jiří Macur, CSc., Úvod do teorie dynamických systémů a jejich simulace, 1. vyd., VUT v Brně, PC-DIR, spol. s r.o. - Nakladatelství Brno, prosinec 1995
    2. Iain G. Main, Kmity a vlny ve fyzice, 1. vyd., Academia Praha 1990, TISK, n. p., Brno

    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.