Application Center - Maplesoft

App Preview:

The elastic prismatic beams on the Winkler foundation

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

[Inserted Image]

The elastic prismatic beams on the Winkler foundation

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

Chair of Mechanics of Materials,
Technical University of L?dz,

Al. Politechniki 6, 93-590 A?dz, POLAND,

email:
Marcin.Kaminski@p.lodz.pl

L?dz, December 2007

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

Homogeneous beam on homogeneous foundation under uniforlmy distributed loading

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

A deflection of the elastic prismatic beam on the elastic Winkler soil is analyzed here. The differential equation governing this problem has the following form:

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

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

where  

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

beta = 2^(1/2)*(k/EJ)^(1/4)/2

The entire length of this beam is partitioned into three different subsets, which reflects the changes in boundary conditions, various elastic constants of the subsoil as well as the beam stiffness variations. The boundary conditions bc1,bc2 reflect the situation at the left edge and bc3, bc4 - at the right end. The remaining boundary conditions are not neccessary here.

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

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

Now we determine the general solutions for those ordinary differential equations using ' dsolve' command as  

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

y(x) = -q/(4*beta^4*EJ)+_C1*exp(2^(1/2)*beta*x)+_C2*sin(2^(1/2)*beta*x)+_C3*cos(2^(1/2)*beta*x)+_C4*exp(-2^(1/2)*beta*x)

Let us introduce the boundary conditions for the left edge of the beam typical for free-end

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

bc1 := 1/4*(-q+4*_C1*beta^4*EJ+4*_C2*sin(0)*beta^4*EJ+4*_C3*cos(0)*beta^4*EJ+4*_C4*beta^4*EJ)/(beta^4*EJ) = 0

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

bc2 := 1/4*(-q+4*_C1*exp(2^(1/2)*beta*L)*beta^4*EJ+4*_C2*sin(2^(1/2)*beta*L)*beta^4*EJ+4*_C3*cos(2^(1/2)*beta*L)*beta^4*EJ+4*_C4*exp(-2^(1/2)*beta*L)*beta^4*EJ)/(beta^4*EJ) = 0

and, further, the boundary conditions for the right edge of the beam - the same as before

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

bc3 := -2*beta^2*(-_C1+_C2*sin(0)+_C3*cos(0)-_C4) = 0

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

bc4 := -2*beta^2*(-_C1*exp(2^(1/2)*beta*L)+_C2*sin(2^(1/2)*beta*L)+_C3*cos(2^(1/2)*beta*L)-_C4*exp(-2^(1/2)*beta*L)) = 0

Let us determine all integration constants for our problem as a solution for the boundary conditions applied  

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

The deflection functions have the following forms (after MAPLE simplifications):

> simplify(y(x)):

We need to introduce the specific values of the beam parameters like k, m, EJ, L before the graphs can be made. Before most of the graphs, the maximum and minimum of each functions (together with their locations along the beam) are found numerically thanks the 'Optimization' library

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

The graph of the beam deflections

> 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=`Deflection of the beam`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2});

The graph of the slope along the beam

> s:=diff(y,x): plot(s,x=0..L, title=`Slope of the beam`,axes=boxed);

The graph of the bending moment  

> 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=`Bending moment`,axes=boxed);

The graph of the shear force

> 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=`Shear force in the beam`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2});

[0.1315640565e-1, [x = .5000000000]]

[0.5017836702e-9, [x = 0.1192092896e-7]]

[Plot]

[Plot]

[.1263375217, [x = .5000000105]]

[0.6010649367e-8, [x = 0.1192092896e-7]]

[Plot]

[.5042098116, [x = 0.1192092896e-7]]

[-.5042098236, [x = 1.]]

[Plot]

As it is clear from those graphs, the deflection function and the bending moment variability are symmetric with respect to the axis x=L/2, whereas the slope as well as the shear force - must be and are antysymmetric. Since slope equals 0 for this particular value, the deflection is maximized at this point (the bending moment behaves analogously).

A homogeneous beam loaded with the concentrated force, concentrated moment and two different distributed loads

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

The entire length of the beam is partitioned into two intervals (with different differential equations), which reflects the changes in boundary conditions, various elastic constants of the subsoil as well as the beam stiffness variations. The essential change with respect to the previous analysis is application of the spatially varying load q=q(x) for both intervals - q1(x) and q2(x) separately. The differential equations valid for those partitions are de1 and de2 and they return the constants _C1,..., _C8 (4 constants for each interval). The boundary conditions bc1,bc2 reflect the situation at the left edge and bc3, bc4 - at the right end, whi;e bc5, bc6, bc7, bc8 are the remaining boundary conditions equivalent to the continuity of this structure for x=L1.

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

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;

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

We need to define first the loadings q1(x) and q2(x) according to the static scheme of this beam. We have

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

Now we determine the general solutions for those ordinary differential equations using ' dsolve' command as  

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

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

y[1](x) = -5*(L1+2*x)/(4*beta1^4*L1*EJ1)+_C1*sin(2^(1/2)*beta1*x)+_C2*exp(2^(1/2)*beta1*x)+_C3*cos(2^(1/2)*beta1*x)+_C4*exp(-2^(1/2)*beta1*x)

y[2](x) = -5/(4*beta2^4*EJ2)+_C1*sin(2^(1/2)*beta2*x)+_C2*exp(2^(1/2)*beta2*x)+_C3*cos(2^(1/2)*beta2*x)+_C4*exp(-2^(1/2)*beta2*x)

We introduce the constants  _C5 .. _C8 instead of  _C1 .. _C4 to define y2(x)  

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

y[2](x) := -5/(4*beta2^4*EJ2)+_C5*sin(2^(1/2)*beta2*x)+_C6*exp(2^(1/2)*beta2*x)+_C7*cos(2^(1/2)*beta2*x)+_C8*exp(-2^(1/2)*beta2*x)

Let us introduce the boundary conditions for the left edge of the beam typical for free-end

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

bc1 := 2*beta1^2*(-_C1*sin(0)+_C2-_C3*cos(0)+_C4) = 0

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

bc2 := 2*2^(1/2)*beta1^3*(-_C1+_C2-_C4) = 0

and, further, the boundary conditions for the right edge of the beam - the same as before

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

bc3 := 2*beta2^2*(-_C5*sin(2^(1/2)*beta2*(L1+L2))+_C6*exp(2^(1/2)*beta2*(L1+L2))-_C7*cos(2^(1/2)*beta2*(L1+L2))+_C8*exp(-2^(1/2)*beta2*(L1+L2))) = 0

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

bc4 := 2*2^(1/2)*beta2^3*(-_C5*cos(2^(1/2)*beta2*(L1+L2))+_C6*exp(2^(1/2)*beta2*(L1+L2))+_C7*sin(2^(1/2)*beta2*(L1+L2))-_C8*exp(-2^(1/2)*beta2*(L1+L2))) = 0

Let us introduce all the continuity conditions connecting the deflection functions y1 and y2 for x=L as well as y2 and y3 for x=2L

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

bc5 := 1/4*(-15+4*_C1*sin(2^(1/2)*beta1*L1)*beta1^4*EJ1+4*_C2*exp(2^(1/2)*beta1*L1)*beta1^4*EJ1+4*_C3*cos(2^(1/2)*beta1*L1)*beta1^4*EJ1+4*_C4*exp(-2^(1/2)*beta1*L1)*beta1^4*EJ1)/(beta1^4*EJ1) = 1/4*(-...bc5 := 1/4*(-15+4*_C1*sin(2^(1/2)*beta1*L1)*beta1^4*EJ1+4*_C2*exp(2^(1/2)*beta1*L1)*beta1^4*EJ1+4*_C3*cos(2^(1/2)*beta1*L1)*beta1^4*EJ1+4*_C4*exp(-2^(1/2)*beta1*L1)*beta1^4*EJ1)/(beta1^4*EJ1) = 1/4*(-...

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

bc6 := 1/2*(-5+2*_C1*cos(2^(1/2)*beta1*L1)*2^(1/2)*beta1^5*L1*EJ1+2*_C2*2^(1/2)*beta1^5*exp(2^(1/2)*beta1*L1)*L1*EJ1-2*_C3*sin(2^(1/2)*beta1*L1)*2^(1/2)*beta1^5*L1*EJ1-2*_C4*2^(1/2)*beta1^5*exp(-2^(1/...bc6 := 1/2*(-5+2*_C1*cos(2^(1/2)*beta1*L1)*2^(1/2)*beta1^5*L1*EJ1+2*_C2*2^(1/2)*beta1^5*exp(2^(1/2)*beta1*L1)*L1*EJ1-2*_C3*sin(2^(1/2)*beta1*L1)*2^(1/2)*beta1^5*L1*EJ1-2*_C4*2^(1/2)*beta1^5*exp(-2^(1/...

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

bc7 := 2*EJ1*beta1^2*(-_C1*sin(2^(1/2)*beta1*L1)+_C2*exp(2^(1/2)*beta1*L1)-_C3*cos(2^(1/2)*beta1*L1)+_C4*exp(-2^(1/2)*beta1*L1)) = EJ1*(M-2*EJ2*_C5*sin(2^(1/2)*beta2*L1)*beta2^2+2*EJ2*_C6*beta2^2*exp(...bc7 := 2*EJ1*beta1^2*(-_C1*sin(2^(1/2)*beta1*L1)+_C2*exp(2^(1/2)*beta1*L1)-_C3*cos(2^(1/2)*beta1*L1)+_C4*exp(-2^(1/2)*beta1*L1)) = EJ1*(M-2*EJ2*_C5*sin(2^(1/2)*beta2*L1)*beta2^2+2*EJ2*_C6*beta2^2*exp(...

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

bc8 := 2*EJ1*2^(1/2)*beta1^3*(-_C1*cos(2^(1/2)*beta1*L1)+_C2*exp(2^(1/2)*beta1*L1)+_C3*sin(2^(1/2)*beta1*L1)-_C4*exp(-2^(1/2)*beta1*L1)) = EJ1*(P-2*EJ2*_C5*cos(2^(1/2)*beta2*L1)*2^(1/2)*beta2^3+2*EJ2*...bc8 := 2*EJ1*2^(1/2)*beta1^3*(-_C1*cos(2^(1/2)*beta1*L1)+_C2*exp(2^(1/2)*beta1*L1)+_C3*sin(2^(1/2)*beta1*L1)-_C4*exp(-2^(1/2)*beta1*L1)) = EJ1*(P-2*EJ2*_C5*cos(2^(1/2)*beta2*L1)*2^(1/2)*beta2^3+2*EJ2*...

Let us determine all integration constants for our problem as a solution for the boundary conditions applied  

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

The deflection functions have the following forms (after MAPLE simplifications):

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

We need to introduce the specific values of the beam parameters like k, m, EJ, L before the graphs can be made.

> 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:

The graph of the beam deflections

> 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=`Deflection of the beam`,axes=boxed);

The graph of the slope along the beam

> 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=`Slope of the beam`,axes=boxed);

The graph of the bending moment must show a jump for x=L1 equal to concentrated moment with the magnitude 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=`Bending moment`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2});

The graph of the shear force shows analogous tendency - look at the coordinate x=L1, where the concentrated force P has been applied

> 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=`Shear force in the beam`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2});

[-0., [x = 2.]]

[-8.049701726, [x = 0.1192092896e-7]]

[Plot]

[Plot]

[.6742511914, [x = 1.000000008]]

[-.3257487901, [x = .9999999972]]

[Plot]

[.4183876821, [x = .2744096507]]

[-2.498208513, [x = .9999999936]]

[Plot]

Civil engineering case study

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

A deflection of the elastic prismatic beam on the elastic Winkler soil is analyzed here. The differential equation governing this problem has the following form:

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

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

where  

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

beta = 2^(1/2)*(k/EJ)^(1/4)/2

The entire length of this beam is partitioned into three different subsets, which reflects the changes in boundary conditions, various elastic constants of the subsoil as well as the beam stiffness variations. The above equation reduces to classical Euler-Bernoulli beam deflection when one drop k=0 for the one of  the intervals along this structure. The differential equations valid for those partitions are de1, de2 and de3. They return the constants _C1,..., _C12, while the boundary conditions bc1,bc2 reflect the situation at the left edge and bc3, bc4 - at the right end. The remaining boundary conditions are equivalent to the continuity conditions for x=L1 and x=L2. As it is apparent from the additional static scheme, the second interval does not contain elastic foundation, so that the corresponding equilibrium differential equation has classical polynomial solution obtained from the straightforward integration. Furthermore, the beam is once more symmetric, so that the deflections and the bending moment functions are symmetric, whereas the slope and the shear force - antysymmetric.

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

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;

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;

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

Now we determine the general solutions for those ordinary differential equations using ' dsolve' command as  

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

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

y[1](x) = -q1/(4*beta1^4*EJ1)+_C1*sin(2^(1/2)*beta1*x)+_C2*exp(2^(1/2)*beta1*x)+_C3*cos(2^(1/2)*beta1*x)+_C4*exp(-2^(1/2)*beta1*x)

y[2](x) = -q2/(4*beta2^4*EJ2)+_C1*exp(2^(1/2)*beta2*x)+_C2*sin(2^(1/2)*beta2*x)+_C3*cos(2^(1/2)*beta2*x)+_C4*exp(-2^(1/2)*beta2*x)

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

y[3](x) = -q3/(4*beta3^4*EJ3)+_C1*cos(2^(1/2)*beta3*x)+_C2*sin(2^(1/2)*beta3*x)+_C3*exp(2^(1/2)*beta3*x)+_C4*exp(-2^(1/2)*beta3*x)

We introduce the constants  _C5 .. _C8 instead of  _C1 .. _C4 to define y2(x) as well as the constants  _C9 .. _C12  instead of  _C1 .. _C4 for the needs of y3(x)

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

y[2](x) := -q2/(4*beta2^4*EJ2)+_C5*exp(2^(1/2)*beta2*x)+_C6*sin(2^(1/2)*beta2*x)+_C7*cos(2^(1/2)*beta2*x)+_C8*exp(-2^(1/2)*beta2*x)

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

y[3](x) := -q3/(4*beta3^4*EJ3)+_C9*cos(2^(1/2)*beta3*x)+_C10*sin(2^(1/2)*beta3*x)+_C11*exp(2^(1/2)*beta3*x)+_C12*exp(-2^(1/2)*beta3*x)

Let us introduce the boundary conditions for the left edge of the beam typical for the free support in the direction perpendicular to its length

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

bc1 := 1/4*(-q1+4*_C1*sin(0)*beta1^4*EJ1+4*_C2*beta1^4*EJ1+4*_C3*cos(0)*beta1^4*EJ1+4*_C4*beta1^4*EJ1)/(beta1^4*EJ1) = 0

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

bc2 := -2*beta1^2*(_C1*sin(0)-_C2+_C3*cos(0)-_C4) = 0

and, further, the boundary conditions for the right edge of the beam - the same as before (free end)

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

bc3 := -2*beta3^2*(_C9*cos(2^(1/2)*beta3*L)+_C10*sin(2^(1/2)*beta3*L)-_C11*exp(2^(1/2)*beta3*L)-_C12*exp(-2^(1/2)*beta3*L)) = 0

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

bc4 := -2*2^(1/2)*beta3^3*(-_C9*sin(2^(1/2)*beta3*L)+_C10*cos(2^(1/2)*beta3*L)-_C11*exp(2^(1/2)*beta3*L)+_C12*exp(-2^(1/2)*beta3*L)) = 0

Let us introduce all the continuity conditions connecting the deflection functions y1 and y2 for x=L as well as y2 and y3 for x=2L

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

bc5 := 1/4*(-q1*beta2^4*EJ2+4*_C1*sin(2^(1/2)*beta1*L1)*beta1^4*EJ1*beta2^4*EJ2+4*_C2*exp(2^(1/2)*beta1*L1)*beta1^4*EJ1*beta2^4*EJ2+4*_C3*cos(2^(1/2)*beta1*L1)*beta1^4*EJ1*beta2^4*EJ2+4*_C4*exp(-2^(1/...bc5 := 1/4*(-q1*beta2^4*EJ2+4*_C1*sin(2^(1/2)*beta1*L1)*beta1^4*EJ1*beta2^4*EJ2+4*_C2*exp(2^(1/2)*beta1*L1)*beta1^4*EJ1*beta2^4*EJ2+4*_C3*cos(2^(1/2)*beta1*L1)*beta1^4*EJ1*beta2^4*EJ2+4*_C4*exp(-2^(1/...

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

bc6 := -2^(1/2)*(-_C1*cos(2^(1/2)*beta1*L1)*beta1-_C2*beta1*exp(2^(1/2)*beta1*L1)+_C3*sin(2^(1/2)*beta1*L1)*beta1+_C4*beta1*exp(-2^(1/2)*beta1*L1)+_C5*beta2*exp(2^(1/2)*beta2*L1)+_C6*cos(2^(1/2)*beta2...bc6 := -2^(1/2)*(-_C1*cos(2^(1/2)*beta1*L1)*beta1-_C2*beta1*exp(2^(1/2)*beta1*L1)+_C3*sin(2^(1/2)*beta1*L1)*beta1+_C4*beta1*exp(-2^(1/2)*beta1*L1)+_C5*beta2*exp(2^(1/2)*beta2*L1)+_C6*cos(2^(1/2)*beta2...

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

bc7 := -2*EJ1*_C1*sin(2^(1/2)*beta1*L1)*beta1^2+2*EJ1*_C2*beta1^2*exp(2^(1/2)*beta1*L1)-2*EJ1*_C3*cos(2^(1/2)*beta1*L1)*beta1^2+2*EJ1*_C4*beta1^2*exp(-2^(1/2)*beta1*L1)-2*EJ2*_C5*beta2^2*exp(2^(1/2)*b...bc7 := -2*EJ1*_C1*sin(2^(1/2)*beta1*L1)*beta1^2+2*EJ1*_C2*beta1^2*exp(2^(1/2)*beta1*L1)-2*EJ1*_C3*cos(2^(1/2)*beta1*L1)*beta1^2+2*EJ1*_C4*beta1^2*exp(-2^(1/2)*beta1*L1)-2*EJ2*_C5*beta2^2*exp(2^(1/2)*b...

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

bc8 := -2*2^(1/2)*(EJ1*_C1*cos(2^(1/2)*beta1*L1)*beta1^3-EJ1*_C2*beta1^3*exp(2^(1/2)*beta1*L1)-EJ1*_C3*sin(2^(1/2)*beta1*L1)*beta1^3+EJ1*_C4*beta1^3*exp(-2^(1/2)*beta1*L1)+EJ2*_C5*beta2^3*exp(2^(1/2)*...bc8 := -2*2^(1/2)*(EJ1*_C1*cos(2^(1/2)*beta1*L1)*beta1^3-EJ1*_C2*beta1^3*exp(2^(1/2)*beta1*L1)-EJ1*_C3*sin(2^(1/2)*beta1*L1)*beta1^3+EJ1*_C4*beta1^3*exp(-2^(1/2)*beta1*L1)+EJ2*_C5*beta2^3*exp(2^(1/2)*...

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

bc9 := 1/4*(-q2*beta3^4*EJ3+4*_C5*exp(2^(1/2)*beta2*L21)*beta2^4*EJ2*beta3^4*EJ3+4*_C6*sin(2^(1/2)*beta2*L21)*beta2^4*EJ2*beta3^4*EJ3+4*_C7*cos(2^(1/2)*beta2*L21)*beta2^4*EJ2*beta3^4*EJ3+4*_C8*exp(-2^...bc9 := 1/4*(-q2*beta3^4*EJ3+4*_C5*exp(2^(1/2)*beta2*L21)*beta2^4*EJ2*beta3^4*EJ3+4*_C6*sin(2^(1/2)*beta2*L21)*beta2^4*EJ2*beta3^4*EJ3+4*_C7*cos(2^(1/2)*beta2*L21)*beta2^4*EJ2*beta3^4*EJ3+4*_C8*exp(-2^...bc9 := 1/4*(-q2*beta3^4*EJ3+4*_C5*exp(2^(1/2)*beta2*L21)*beta2^4*EJ2*beta3^4*EJ3+4*_C6*sin(2^(1/2)*beta2*L21)*beta2^4*EJ2*beta3^4*EJ3+4*_C7*cos(2^(1/2)*beta2*L21)*beta2^4*EJ2*beta3^4*EJ3+4*_C8*exp(-2^...

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

bc10 := -2^(1/2)*(-_C5*beta2*exp(2^(1/2)*beta2*L21)-_C6*cos(2^(1/2)*beta2*L21)*beta2+_C7*sin(2^(1/2)*beta2*L21)*beta2+_C8*beta2*exp(-2^(1/2)*beta2*L21)-_C9*sin(2^(1/2)*beta3*L21)*beta3+_C10*cos(2^(1/2...bc10 := -2^(1/2)*(-_C5*beta2*exp(2^(1/2)*beta2*L21)-_C6*cos(2^(1/2)*beta2*L21)*beta2+_C7*sin(2^(1/2)*beta2*L21)*beta2+_C8*beta2*exp(-2^(1/2)*beta2*L21)-_C9*sin(2^(1/2)*beta3*L21)*beta3+_C10*cos(2^(1/2...

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

bc11 := 2*EJ2*_C5*beta2^2*exp(2^(1/2)*beta2*L21)-2*EJ2*_C6*sin(2^(1/2)*beta2*L21)*beta2^2-2*EJ2*_C7*cos(2^(1/2)*beta2*L21)*beta2^2+2*EJ2*_C8*beta2^2*exp(-2^(1/2)*beta2*L21)+2*EJ3*_C9*cos(2^(1/2)*beta3...bc11 := 2*EJ2*_C5*beta2^2*exp(2^(1/2)*beta2*L21)-2*EJ2*_C6*sin(2^(1/2)*beta2*L21)*beta2^2-2*EJ2*_C7*cos(2^(1/2)*beta2*L21)*beta2^2+2*EJ2*_C8*beta2^2*exp(-2^(1/2)*beta2*L21)+2*EJ3*_C9*cos(2^(1/2)*beta3...

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

bc12 := -2*2^(1/2)*(-EJ2*_C5*beta2^3*exp(2^(1/2)*beta2*L21)+EJ2*_C6*cos(2^(1/2)*beta2*L21)*beta2^3-EJ2*_C7*sin(2^(1/2)*beta2*L21)*beta2^3+EJ2*_C8*beta2^3*exp(-2^(1/2)*beta2*L21)+EJ3*_C9*sin(2^(1/2)*be...bc12 := -2*2^(1/2)*(-EJ2*_C5*beta2^3*exp(2^(1/2)*beta2*L21)+EJ2*_C6*cos(2^(1/2)*beta2*L21)*beta2^3-EJ2*_C7*sin(2^(1/2)*beta2*L21)*beta2^3+EJ2*_C8*beta2^3*exp(-2^(1/2)*beta2*L21)+EJ3*_C9*sin(2^(1/2)*be...

We need to introduce the specific values of the beam parameters like k, m, EJ, L in each interval separately before the final graphs can be made.

> 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):

Let us determine all integration constants for our problem as a solution for the boundary conditions applied  

> 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(%):

We simplify the deflections functions in each interval of the beam. There holds  

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

The graph of the beam deflections - everyone must remember that the real deflection is directed downwards - in opposite direction to that visualized on this graph

> 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=`Deflection of the beam`,axes=boxed);  

The graph of the slope along the beam

> 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=`Slope of the beam`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2});

The graph of the bending moment  

> 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=`Bending moment`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2});

Finally, we make the graph of the shear force

> 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=`Shear force in the beam`,axes=boxed): p2:=plot(0,x=0..L,axes=boxed,colour=black): display({p1,p2});

[10.73514182, [x = 14.99999999]]

[0.2412310907e-9, [x = 0.1238729139e-7]]

[Plot]

[Plot]

[75.53341296, [x = 9.999522726]]

[-96.82806004, [x = 5.000000001]]

[Plot]

[53.58527261, [x = 5.000000008]]

[-26.07249791, [x = 12.67795658]]

[Plot]

Let us note as before, that an application of the concentrated loadings must result in the jumps on the additional diagrams.

Let us mention also that one can apply here the other supports for the beam resting on the elastic foundation like telescope or clamped edge. If some of them are inserted for 0<x<L then two of the equations bc5-bc12 must contain some extra conditions like zeroing of the displacements and moments (for a simple support). When some support is added for x=0 or x=L, the conditions bc1-bc4 must be corrected accordingly by setting y(x)=0 and diff(y(x),x)=0 in the case of the clamped adge.

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.

[Inserted Image]