Ugi?cie belek pryzmatycznych na spr?ystym podlo?u

Prof. Marcin Kaminski, D.Sc., Ph.D.

Chair of Mechanics of Materials,
Technical University of L?dz,
Al. Politechniki 6, 93-590 L?dz, POLAND,
email: Marcin.Kaminski@p.lodz.pl

L?dz, January 2008

Abstract: This script is entirely devoted to the static problem of a deflection of the linear elastic istotropic prismatic beams resting on the homogeneous linear elastic and isotropic foundation. It is demonstrated below how one can symbolically solve the problems of such beams consisting of one, two and three different intervals with spatially varying distributed load applied along the beam, different concentrated forces and moments. It is also possible to define some variability of the bending stiffness as well as the elastic foundation parameter as piecewise constant functions along the beam. This program is able to solve symbolically fourth order differential equations of the beam equilibrium, then - to determine the integration constants from the additional boundary and continuity conditions for the beam intervals. Finally, it computes maxima and minima of deflections, bending moments and shear forces and prepares the diagrams of all those functions. This script may be easily modified to extend the beam analysis towards the structural sensitivity studies and/or reliability analysis for the beams resting on the elastic subsoils.

Keywords: elastic beams, elastic foundation, ordinary differential equations, beam diagrams

 Belka jednorodna na jednorodnym podlo?u poddana obci??eniu stalemu rownomiernie rozlo?onemu

>	restart: with(plots): with(Optimization):

Warning, the name changecoords has been redefined

Analizujemy tutaj lini? ugi?cia dla belki pryzmatycznej o stalym przekroju le??cej na winklerowskim spr?ystym podlo?u. Material, z kt?rego wykonana jest belka jest r?wnie? liniowo-spr?ysty i izotropowy. R?wnanie r?niczkowe opisuj?ce poszukiwane ugi?cie ma posta? nast?puj?c?:  

>	Diff(y(x),x$4)-4*beta^4*y(x)-q/EJ=0;

gdzie   

>	beta=sqrt(sqrt(k/4/EJ));

Warunki brzegowe niezb?dne do obliczenia stalych calkowania to bc1,bc2 odpowiadaj?ce lewemu ko?cowi belki oraz bc3, bc4 dla prawego ko?ca belki. Inne warunki brzegowe, jak np. zerowanie si? k?ta ugi?cia na osi symetrii belki nie s? tutaj konieczne.  

>	de:=diff(y(x),x$4)-4*beta^4*y(x)-q/EJ=0;

Znajdujemy najpierw rozwi?zanie og?lne dla r?wnania r?niczkowego linii ugi?cia

>	dsolve(de,y(x)); assign(%);

Kolejno stosujemy warunki brzegowe dla lewego brzegu wla?ciwe dla swobodnego ko?ca  

>	bc1:=simplify(subs(x=0,y(x))=0);

>	bc2:=simplify(subs(x=L,y(x))=0);

a tak?e warunki brzegowe dla ko?ca prawego

>	bc3:=simplify(subs(x=0,diff(y(x),x$2))=0);

>	bc4:=simplify(subs(x=L,diff(y(x),x$2))=0);

Potem wyznaczamy warto?ci tych stalych poprzez rozwiazanie nastepuj?cego ukladu r?wna?:   

>	solve({bc1,bc2,bc3,bc4},{_C1,_C2,_C3,_C4}): assign(%):

Uzyskujemy w ten spos?b ostateczn? lini? ugi?cia dla rozpatrywanej belki.

>	simplify(y(x)):

W celu sporz?dzenia wykres?w wprowadzamy okre?lone warto?ci parametr?w gruntu oraz belki - k, m, EJ oraz L. Znajdujemy jednocze?nie warto?ci liczbowe odpowiadaj?ce ekstremom poszczeg?lnych funkcji.

>	k:=1: q:=1: EJ:=1: L:=1: beta:=sqrt(sqrt(k/4/EJ)):

Wykres ugi?cia belki jest nast?puj?cy:

>	y:=y(x): fmax:=Maximize(y,x=0..L): evalf(fmax); fmin:=Minimize(y,x=0..L): evalf(fmin); p1:=plot(-y,x=0..L, title=`Ugi?cie belki`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2},font=[TIMES, BOLD, 16]);

Kolejno sporz?dzamy wykres k?t?w ugi?cia dla tej belki

>	s:=diff(y,x): plot(s,x=0..L, title=`K?t ugi?cia belki`,axes=boxed,font=[TIMES, BOLD, 16]);

Nast?pnie mamy wykres momentu zginaj?cego

>	m:=-EJ*diff(y,x$2): Mmax:=Maximize(m,x=0..L): evalf(Mmax); Mmin:=Minimize(m,x=0..L): evalf(Mmin); plot(m,x=0..L, title=`Moment zginaj?cy`,axes=boxed,font=[TIMES, BOLD, 16]);

oraz wykres sily tn?cej na dlugo?ci belki

>	t:=-EJ*diff(y,x$3): Tmax:=Maximize(t,x=0..L): evalf(Tmax); Tmin:=Minimize(t,x=0..L): evalf(Tmin); p1:=plot(t,x=0..L, title=`Sila tn?ca`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2},font=[TIMES, BOLD, 16]);

Jak wynika z tych wykres?w, funkcja ugi?cia i momentu zginaj?cego s? symetryczne wzgl?dem osi

x=L/2, podczas gdy wykres k?ta ugi?cia oraz sily tn?cej jako pochodne rz?du nieparzystego s? antysymetryczne wzgl?dem tej samej prostej. Warto?? maksymaln? ugi?cia otrzymujemy oczywi?cie dla warto?ci x, dla kt?rej k?t jest r?wny 0 i zmienia znak w otoczeniu tego punktu (moment zginaj?cy i sia tn?ca wykazuj? podobn? zale?no??).

Jednorodna belka obci??ona sil? skupion? oraz skupionym momentem oraz dwoma rodzajami obcia?enia r?wnomiernie rozlo?onego

>	restart: with(plots): with(Optimization):

Warning, the name changecoords has been redefined

Ze wzgl?du na zmienno?? obci??enia na dlugo?ci belki dzielimy j? na dwa rozlaczne przedzialy (gdzie obowi?zuj? r?ne r?wnania r?niczkowe ugi?cia), co odzwierciedla zmiany w sztywno?ci gi?tnej, parametrach podlo?a a tak?e r?ne warunki brzegowe. Zasadnicz? r?nic? w stosunku do poprzednio rozwi?zanego zadania jest obecno?? obci??e? r?wnomiernie rozlo?onych na dlugo?ci przedzial?w - trapezowego q1(x) (staego + trojkatnego) na pierwszym odcinku oraz stalego q2(x) na odcinku drugim. R?wnania r?niczkowe ugi?cia oznaczamy jako de1 oraz de2, w kt?rych wyst?puj? stale  _C1,..., _C8 (po 4 stale dla ka?dego odcinka). Warunki brzegowe bc1,bc2 odzwierciedlaj? sytuacj? na lewym ko?cu schematu statycznego, natomiast bc3, bc4 - na prawym ko?cu ukladu. W warunkach bc5, bc6, bc7, bc8 rozpoznajemy warunki zszycia dla punktu wsp?lnego obydwu linii ugi?cia, czyli dla x=L1.

>	de1:=diff(y[1](x),x$4)-4*beta1^4*y[1](x)-q1(x)/EJ1=0;

>	de2:=diff(y[2](x),x$4)-4*beta2^4*y[2](x)-q2(x)/EJ2=0;

Definiujemy obci??enia q1(x) oraz q2(x) zgodnie z wprowadzonym schematem statycznym. Mamy

>	q0:=5: q1(x):=q0*(1+2*x/L1): q2(x):=q0:

Kolejno otrzymujemy rozwi?zania og?lne dla tych rowna? ro?niczkowych jako  

>	dsolve(de1,y[1](x)); assign(%);

>	dsolve(de2,y[2](x)); assign(%);

Wprowadzamy nowe stale calkowania dla drugiej linii ugi?cia:  _C5 .. _C8  zamiast  _C1 .. _C4

>	y[2](x):=subs(_C1=_C5,_C2=_C6,_C3=_C7,_C4=_C8,y[2](x));

Stosujemy kolejno warunki brzegowe typowe dla swobodnego ko?ca na lewym brzegu  

>	bc1:=simplify(subs(x=0,diff(y[1](x),x$2))=0);

>	bc2:=simplify(subs(x=0,diff(y[1](x),x$3))=0);

i kolejno takie same warunki dla prawego brzegu

>	bc3:=simplify(subs(x=L1+L2,diff(y[2](x),x$2))=0);

>	bc4:=simplify(subs(x=L1+L2,diff(y[2](x),x$3))=0);

Nast?pnie wprowadzamy warunki ci?glo?ci poszczeg?lnych funkcji dla x=L, w kt?rych uwzgl?niamy skokow? zmian? momentu oraz sily tn?cej w dw?ch ostatnich warunkach

>	bc5:=simplify(subs(x=L1,y[1](x))=subs(x=L1,y[2](x)));

>	bc6:=simplify(subs(x=L1,diff(y[1](x),x))=subs(x=L1,diff(y[2](x),x)));

>	bc7:=EJ1*simplify(subs(x=L1,diff(y[1](x),x$2))=M+EJ2*subs(x=L1,diff(y[2](x),x$2)));

>	bc8:=EJ1*simplify(subs(x=L1,diff(y[1](x),x$3))=P+EJ2*subs(x=L1,diff(y[2](x),x$3)));

Ostatecznie wyznaczamy warto?ci wszystkich stalych calkowania jako rozwi?zanie dla nast?puj?cego ukladu r?wna?  

>	solve({bc1,bc2,bc3,bc4,bc5,bc6,bc7,bc8},{_C1,_C2,_C3,_C4,_C5,_C6,_C7,_C8}): assign(%):

Upraszczamy posta? funkcji ugi?cia  

>	simplify(y[1](x)): simplify(y[2](x)):

oraz definiujemy warto?ci parametr?w k, m, EJ, L w celu wykonania odpowiednich rysunk?w

>	k1:=1: q1:=1: EJ1:=1: L1:=1: k2:=1: q2:=1: EJ2:=1: L2:=1: beta1:=sqrt(sqrt(k1/4/EJ1)): beta2:=sqrt(sqrt(k2/4/EJ2)): P:=1: M:=1: L:=L1+L2:

Otrzymujemy wykres linii ugi?cia

>	y1:=y[1](x): y2:=y[2](x): deflection:=piecewise(x<L1,y1,x<L,y2): fmax:=Maximize(deflection,x=0..L): evalf(fmax); fmin:=Minimize(deflection,x=0..L): evalf(fmin); plot(deflection,x=0..L, title=`Linia ugi?cia belki`,axes=boxed,axes=boxed,font=[TIMES, BOLD, 16]);

kolejno znajdujemy wykres k?t?w ugi?cia dla rozpatrywanej belki

>	y1:=diff(y[1](x),x): y2:=diff(y[2](x),x): Slope:=piecewise(x<L1,y1,x<L,y2):

>	plot(Slope,x=0..L, title=`K?t ugi?cia belki`,axes=boxed,axes=boxed,font=[TIMES, BOLD, 16]);

Wykres momentu zginaj?cego musi mie? nieci?glo?? w miejscu polaczenia y1(x) oraz y2(x) o wielko?ci M  

>

y1:=-EJ1*diff(y[1](x),x$2): y2:=-EJ2*diff(y[2](x),x$2): moment:=piecewise(x<L1,y1,x<L,y2): Mmax:=Maximize(moment,x=0..L): evalf(Mmax); Mmin:=Minimize(moment,x=0..L): evalf(Mmin); p1:=plot(moment,x=0..L, title=`Moment zginajacy`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2},axes=boxed,font=[TIMES, BOLD, 16]);

Wykres sily tn?cej musi wykazywa? podobn? wla?ciwo?? - skok o wielko?? P dla wsp?lrz?dnej x=L1, gdzie przylo?ono obci??enie zewn?trzne

>	y1:=-EJ1*diff(y[1](x),x$3): y2:=-EJ2*diff(y[2](x),x$3): shear:=piecewise(x<L1,y1,x<L,y2): Tmax:=Maximize(shear,x=0..L): evalf(Tmax); Tmin:=Minimize(shear,x=0..L): evalf(Tmin);

>	p1:=plot(shear,x=0..L, title=`Sila tn?ca`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2},axes=boxed,font=[TIMES, BOLD, 16]);

Belka wieloprz?slowa na spr?ystym podlo?u  

>	restart: with(plots): with(Optimization): with(LinearAlgebra):

Warning, the name changecoords has been redefined

Na zako?czenie rozpatrzymy lini? ugi?cia belki tr?jprz?slowej pokazanej na powy?szym rysunku, gdzie dla ka?dego odcinka obowi?zuje analogiczne r?wnanie linii ugi?cia:

>	Diff(y(x),x$4)-4*beta^4*y(x)-q/EJ=0;

i gdzie oczywi?cie jest

>	beta=sqrt(sqrt(k/4/EJ));

W tym przypadku dokonujemy podzialu calej dlugo?ci belki na trzy rozl?czne podprzedzialy w celu uwzgl?dnienia zmian sztywno?ci, obci??en i warunkow gruntowych. Zauwazmy, ze jak poprzednio, r?wnanie to redukuje si? do klasycznego r?wnania linii ugi?cia belki zgodnie z teori? rz?du pierwszego poprzez podstawienie k=0. R?wnania r?niczkowe odpowiadaj?ce temu podzialowi to odpowiednio de1, de2 oraz de3. Ich czterokrotne calkowanie powoduje pojawienie sie stalych calkowania C1,..., _C12, kt?rych warto?ci znajdujemy z warunkow brzegowych bc1,bc2 dla lewego ko?ca struktury, warunk?w bc3 oraz bc4 - dla prawego ko?ca belki oraz warunk?w zszycia odpowiadaj?cych wsp?lrz?dnym x=L1 oraz x=L2.  

>	de1:=diff(y[1](x),x$4)-4*beta1^4*y[1](x)-q1/EJ1=0;

>	de2:=diff(y[2](x),x$4)-4*beta2^4*y[2](x)-q2/EJ2=0;

>	de3:=diff(y[3](x),x$4)-4*beta3^4*y[3](x)-q3/EJ3=0;

Kolejno wyznaczamy rozwi?zania og?lne dla tych r?wna? jako  

>	dsolve(de1,y[1](x)); assign(%);

>	dsolve(de2,y[2](x)); assign(%);

>	dsolve(de3,y[3](x)); assign(%);

Wprowadzamy stale  _C5 .. _C8 zamiast  _C1 .. _C4 dla okre?lenia funkcji y2(x) oraz stale  _C9 .. _C12 zamiast  _C1 .. _C4 - w celu zdefiniowania linii ugi?cia y3(x)

>	y[2](x):=subs(_C1=_C5,_C2=_C6,_C3=_C7,_C4=_C8,y[2](x));

>	y[3](x):=subs(_C1=_C9,_C2=_C10,_C3=_C11,_C4=_C12,y[3](x));

Wprowadzamy warunki brzegowe typowe dla swobodnego ko?ca na lewym ko?cu belki  

>	bc1:=simplify(subs(x=0,diff(y[1](x),x$2))=0);

>	bc2:=simplify(subs(x=0,diff(y[1](x),x$3))=0);

oraz warunki na prawym brzegu typowe dla swobodnego ko?ca

>	bc3:=simplify(subs(x=L,diff(y[3](x),x$2))=0);

>	bc4:=simplify(subs(x=L,diff(y[3](x),x$3))=0);

Kolejno wprowadzamy warunki ci?glo?ci dla poszczeg?lnych pochodnych funkcji ugi?? y1(x) i y2(x) dla wsp?lrz?dnych x=L oraz zszycia funkcji y2 oraz y3 w x=2L

>	bc5:=simplify(subs(x=L1,y[1](x))-subs(x=L1,y[2](x))=0);

>	bc6:=simplify(subs(x=L1,diff(y[1](x),x))-subs(x=L1,diff(y[2](x),x))=0);

>	bc7:=simplify(subs(x=L1,EJ1*diff(y[1](x),x$2))-subs(x=L1,EJ2*diff(y[2](x),x$2))=0);

>	bc8:=simplify(subs(x=L1,EJ1*diff(y[1](x),x$3))-subs(x=L1,EJ2*diff(y[2](x),x$3))=P);

>	bc9:=simplify(subs(x=L21,y[2](x))-subs(x=L21,y[3](x))=0);

>	bc10:=simplify(subs(x=L21,diff(y[2](x),x))-subs(x=L21,diff(y[3](x),x))=0);

>	bc11:=simplify(subs(x=L21,EJ2*diff(y[2](x),x$2))-subs(x=L21,EJ3*diff(y[3](x),x$2))=M);

>	bc12:=simplify(subs(x=L21,EJ2*diff(y[2](x),x$3))-subs(x=L21,EJ3*diff(y[3](x),x$3))=0);

Definiujemy warto?ci parametr?w k, m, EJ, L oddzielnie w ka?dym z wprowadzonych przedzial?w w celu wykonania odpowiednich wykres?w  

>	k1:=10: q1:=10: EJ1:=2000: L1:=5: beta1:=sqrt(sqrt(k1/4/EJ1)): k2:=10: q2:=0: EJ2:=210: L2:=6: beta2:=sqrt(sqrt(k2/4/EJ2)): k3:=10: q3:=0: EJ3:=300: L3:=4: beta3:=sqrt(sqrt(k3/4/EJ3)): P:=100.0: M:=10.0: L:=L1+L2+L3: L21:=L1+L2:

>	evalf(bc1): evalf(bc2): evalf(bc3): evalf(bc4): evalf(bc5): evalf(bc6): evalf(bc7): evalf(bc8): evalf(bc9): evalf(bc10): evalf(bc11): evalf(bc12):

Rozwi?zujemy ukad r?wna? wynikaj?cy z warunk?w brzegowych w celu wyznaczenia stalych calkowania i mamy   

>	fsolve({bc1,bc2,bc3,bc4,bc5,bc6,bc7,bc8,bc9,bc10,bc11,bc12},{_C1,_C2,_C3,_C4,_C5,_C6,_C7,_C8,_C9,_C10,_C11,_C12}): assign(%):

Upraszczamy jednocze?nie r?wnania linii ugi?cia w poszczeg?lnych przedzialach   

>	simplify(y[1](x)): simplify(y[2](x)): simplify(y[3](x)):

Kolejno otrzymujemy wykres linii ugi?cia

>	y1:=y[1](x): y2:=y[2](x): y3:=y[3](x): deflection:=piecewise(x<L,y1,x<L1+L2,y2,x<L,y3): fmax:=Maximize(deflection,x=0..L): evalf(fmax); fmin:=Minimize(deflection,x=0..L): evalf(fmin); plot(deflection,x=0..L, title=`Linia ugiecia belki`,axes=boxed,font=[TIMES, BOLD, 16]);

oraz wykres funkcji k?ta ugi?cia  

>	y1:=diff(y[1](x),x): y2:=diff(y[2](x),x): y3:=diff(y[3](x),x): Slope:=piecewise(x<L1,y1,x<L1+L2,y2,x<L,y3):

>	p1:=plot(Slope,x=0..L, title=`Kat ugiecia belki`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2},font=[TIMES, BOLD, 16]);

wykres momentu zginaj?cego   

>

y1:=-EJ1*diff(y[1](x),x$2): y2:=-EJ2*diff(y[2](x),x$2): y3:=-EJ3*diff(y[3](x),x$2): moment:=piecewise(x<L1,y1,x<L1+L2,y2,x<L,y3): Mmax:=Maximize(moment,x=0..L): evalf(Mmax); Mmin:=Minimize(moment,x=0..L): evalf(Mmin); p1:=plot(moment,x=0..L, title=`Wykres momentu zginajacego`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2},font=[TIMES, BOLD, 16]);

oraz sily tn?cej  

>	y1:=-EJ1*diff(y[1](x),x$3): y2:=-EJ2*diff(y[2](x),x$3): y3:=-EJ3*diff(y[3](x),x$3): shear:=piecewise(x<L1,y1,x<L1+L2,y2,x<L,y3): Tmax:=Maximize(shear,x=0..L): evalf(Tmax); Tmin:=Minimize(shear,x=0..L): evalf(Tmin);

>	p1:=plot(shear,x=0..L, title=`Wykres sily tnacej`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2},font=[TIMES, BOLD, 16]);

Zauwa?my, ?e tak jak w poprzednim przypadku, zastosowanie obci??e? skupionych musi wplywa? na skokow? zmienno?? odpowiedniego wykresu. Dodajmy, ?e mo?emy zastosowa? tutaj r?wnie? inne warunki podparcia wla?ciwe dla zwyklych belek. I tak np. mo?emy zastosowa? tutaj teleskop albo zamocowanie (utwierdzenie) - je?li jakakolwiek inna podpora pojawi si? dla 0<x<L, w?wczas dwa z r?wna? bc5-bc12 musi zawiera? dodatkowo warunki zerowania si? przemieszcze? i moment?w (jak dla przegubowego podparcia). Gdy stosujemy podpor? dla x=0 lub x=L, w?wczas warunki brzegowe bc1-bc4 musz? by? odpowiednio zmienione poprzez np. wprowadzenie y(x)=0 oraz  diff(y(x),x)=0 w przypadku utwierdzenia jednego z ko?c?w belki.

Reference

S. Timoshenko, J.N. Goodier, Theory of Elasticity, second edition, McGraw-Hill Book Company Ltd., New York-Toronto-London, 1951.

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.