Open Educational Resources

Approximate Methods: Rayleigh Ritz Method

The Rayleigh Ritz method is a classical approximate method to find the displacement function of an object such that the it is in equilibrium with the externally applied loads. It is regarded as an ancestor of the widely used Finite Element Method (FEM). The Rayleigh Ritz method relies on the principle of minimum potential energy for conservative systems. The method involves assuming a form or a shape for the unknown displacement functions, and thus, the displacement functions would have a few unknown parameters. These assumed shape functions are termed Trial Functions. Afterwards, the potential energy function of the system is written in terms of those few parameters, and the values of those parameters that would minimize the potential energy function of the system are calculated. Mathematically speaking, we assume that the unknown displacement function u is a member of a certain space of functions; for example, we can assume that u\in V where V is the set of all the possible vector valued linear functions defined on the body represented by the set \Omega_0\subset \mathbb{R}^3:V=\{u:\Omega_0\rightarrow \mathbb{R}^3|\forall x\in\Omega_0:u(x)=Bx, B\in\mathbb{M}^3\}. With this assumption, we restrict the possible solutions to those in the form u=Bx, and thus, the unknown parameters are restricted to the matrix B, which in this particular example, has nine components. Finally, the nine components that would minimize the potential-energy function are obtained. The method will be illustrated using various one-dimensional examples.

Bars under Axial Loads

The Rayleigh Ritz method seeks to find an approximate solution to minimize the potential energy of the system. For plane bars under axial loading, the unknown displacement function is u_1. In this case, the potential energy of the system has the following form (See beams under axial loads and energy expressions):

    \[ PE=\mbox{Total Strain Energy} - W=\int\!\frac{EA}{2}\left(\frac{\mathrm{d}u_1}{\mathrm{d}X_1}\right)^2 \,\mathrm{d}X_1 - W \]

where E is Young’s modulus, A is the cross sectional area, and W is the work done by all the external forces acting on the bar (concentrated forces and distributed loads). To obtain an approximate solution for u_1, a trial function of a finite number of unknowns is first selected. Then, this trial function has to satisfy the “essential” boundary conditions, which are those boundary conditions that would satisfy the physical stability of the system. In the case of bars under axial loads, the trial function u_1 has to satisfy the displacement boundary conditions. The following examples show the method for different external loads, and bar configurations.

Example 1: A Bar under Axial Forces

Let a bar with a length L and a constant cross sectional area A be aligned with the coordinate axis X_1. Assume that the bar is subjected to a horizontal body force (in the direction of the axis X_1) p=CX_1 units of force/unit length where C is a known constant. Also, assume that the bar has a constant force of value P applied at the end X_1 = L. If the bar is fixed at the end X_1 = 0, find the displacement of the bar by directly solving the differential equation of equilibrium. Also, find the displacement function using the Rayleigh Ritz method assuming a polynomial function for the displacement of the degrees (0, 1, 2, 3, and 4). Assume that the bar is linear elastic with Young’s modulus E and that the small strain tensor is the appropriate measure of strain. Ignore the effect of Poisson’s ratio.
RR1

Solution

Exact Solution:
First, the exact solution that would satisfy the equilibrium equations can be obtained. The equilibrium equation as shown in the beams under axial loading section when E and A are constant is:

    \[ \frac{\mathrm{d}^2u_1}{\mathrm{d}X_1^2}=-\frac{CX_1}{EA} \]

Therefore, the exact solution has the form:

    \[ u_1=-\frac{CX_1^3}{6EA}+C_1X_1 + C_2 \]

where C_1 and C_2 can be obtained from the boundary conditions:

    \[ \begin{split} & @X_1=0:u_1=0\Rightarrow C_2=0\\ & @X_1=L:\sigma_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=\frac{P}{A}\Rightarrow C_1=\frac{P}{AE}+\frac{CL^2}{2EA} \end{split} \]

Therefore, the exact solution for the displacement u_1 and the stress \sigma_{11} are:

    \[\begin{split} u_1 & =-\frac{CX_1^3}{6EA}+\left(\frac{P}{AE}+\frac{CL^2}{2EA}\right)X_1 \\ \sigma_{11}& = E\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=-\frac{CX_1^2}{2A}+\left(\frac{P}{A}+\frac{CL^2}{2A}\right) \end{split} \]

View Mathematica Code
Clear[u, x]
a = DSolve[{u''[x] == -c*x/(Ee*A), u'[L] == P/(Ee*A), u[0] == 0}, u, x]
u = u[x] /. a[[1]]
FullSimplify[u]
sigma = Ee*D[u, x];
FullSimplify[sigma]

The Rayleigh Ritz Method:
The first step in the Rayleigh Ritz method finds the minimizer of the potential energy of the system which can be written as:

    \[ PE=\int\!\overline{U}\,dX - W=\int_0^L\!\frac{\sigma_{11}\varepsilon_{11}}{2}A \,\mathrm{d}X_1 - Pu_1|_{X_1=L}-\int_0^L\!pu_1\,\mathrm{d}X_1 \]

Notice that the potential energy lost by the action of the end force P is equal to the product of P and the displacement u_1 evaluated at X_1= L while the potential energy lost by the action of the distributed body forces is an integral since it acts on each point along the beam length. Further, we will use the constitutive equation to rewrite the potential energy function in terms of the function u_1:

    \[ PE=\int_0^L\!\frac{EA}{2}\left(\frac{\mathrm{d}u_1}{\mathrm{d}X_1}\right)^2 \,\mathrm{d}X_1 - Pu_1|_{X_1=L}-\int_0^L\!CX_1u_1\,\mathrm{d}X_1 \]

To find an approximate solution, an assumption for the shape or the form of u_1 has to be introduced.
Polynomial of the Zero Degree:
First, we can assume that u_1 is a constant function (a polynomial of the zero degree):

    \[ u_1 = a_0 \]

However, to satisfy the “essential” boundary condition that u_1 = 0 when X_1 = 0, the constant a_0 has to equal to zero; this leads to the trivial solution u_1= 0, which will automatically be rejected.
Polynomial of the First Degree:
As a more refined approximation, we can assume that u_1 is a linear function (a polynomial of the first degree):

    \[ u_1 = a_0+a_1 X_1 \]

To satisfy the essential boundary condition @X_1=0:u_1=0 the coefficient a_0 is equal to zero. Therefore:

    \[ u_1 = a_1 X_1 \]

The strain component \varepsilon_{11} can be computed as follows:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=a_1 \]

Substituting into the equation for the potential energy of the system:

    \[ PE=\int_0^L\!\frac{EA}{2}a_1^2 \,\mathrm{d}X_1 - Pa_1L-\int_0^L\!CX_1a_1X_1\,\mathrm{d}X_1=\frac{EAL}{2}a_1^2-Pa_1L-\frac{Ca_1L^3}{3} \]

The only variable that can be controlled to minimize the potential energy of the system is a_1. Therefore:

    \[ \frac{\mathrm{d}PE}{\mathrm{d}a_1}=0\Rightarrow EAa_1L-PL-\frac{CL^3}{3}=0\Rightarrow a_1=\frac{P+\frac{CL^2}{3}}{EA} \]

Therefore, the best linear function according to the Rayleigh Ritz method is:

    \[ u_1=\left(\frac{P+\frac{CL^2}{3}}{EA}\right)X_1 \]

A linear displacement would produce a constant stress:

    \[ \sigma_{11}=\left(\frac{P+\frac{CL^2}{3}}{A}\right) \]

Unfortunately, this solution is far from being accurate, since it does not satisfy the differential equation of equilibrium (Why?). However, within the only possible or allowed linear shapes, the obtained solution is the minimizer of the potential energy.

Polynomial of the Second Degree:
As a more refined approximation, we can assume that u_1 is a polynomial of the second degree:

    \[ u_1 = a_0+a_1 X_1+a_2X_1^2 \]

To satisfy the essential boundary condition @X_1=0:u_1=0 the coefficient a_0 is equal to zero. Therefore:

    \[ u_1 = a_1 X_1+a_2X_1^2 \]

The strain component \varepsilon_{11} can be computed as follows:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=a_1+2a_2X_1 \]

Substituting into the equation for the potential energy of the system:

    \[\begin{split} PE&=\int_0^L\!\frac{EA}{2}(a_1+2a_2X_1)^2 \,\mathrm{d}X_1 - P(a_1L+a_2L^2)-\int_0^L\!C(a_1X_1^2+a_2X_1^3)\,\mathrm{d}X_1\\ & = \frac{2EAL^3}{3}a_2^2+EAL^2a_2a_1+\frac{EAL}{2}a_1^2-PL^2a_2-PLa_1-\frac{CL^4}{4}a_2-\frac{CL^3}{3}a_1 \end{split} \]

a_1 and a_2 are the variables that can be controlled to minimize the potential energy of the system. Therefore:

    \[ \frac{\partial PE}{\partial a_1}=0\Rightarrow EAL^2a_2+EALa_1-PL-\frac{CL^3}{3}=0 \]

and

    \[ \frac{\partial PE}{\partial a_2}=0\Rightarrow \frac{4EAL^3}{3}a_2+EAL^2a_1-PL^2-\frac{CL^4}{4}=0 \]

Solving the above two equations yields:

    \[ a_1=\frac{7CL^2+12P}{12EA}\qquad a_2=-\frac{CL}{4EA} \]

Therefore, the best second degree polynomial according to the Rayleigh Ritz method is:

    \[ u_1=-\frac{cL}{4EA}X_1^2+\frac{(7CL^2+12P)}{12EA}X_1 \]

A parabolic displacement would produce a linear stress:

    \[ \sigma_{11}=-\frac{cL}{2A}X_1+\frac{(7CL^2+12P)}{12A} \]

View Mathematica Code
Clear[u,x,a2,a1,Ee,A,c,P,L]  
u=a2*x^2+a1*x;
PE=Integrate[(1/2)*Ee*A*D[u,x]^2,{x,0,L}]-P*(u/.x->L)-Integrate[c*x*u,{x,0,L}];
Eq1=D[PE,a2]  
Eq2=D[PE,a1]  
s=Solve[{Eq1==0,Eq2==0},{a2,a1}]  
u=u/.s[[1]]  
stress=FullSimplify[Ee*D[u,x]]  

Polynomial of the Third Degree:
Similarly, choosing a polynomial of the third degree as an approximate displacement function means:

    \[ u_1=a_0+a_1X_1+a_2X_1^2+a_3X_1^3 \]

with the following approximation for the axial strain

    \[ \varepsilon_{11}=a_1+2a_2X_1+3a_3X_1^2 \]

The essential boundary condition lead to a_0=0. The potential energy of the system has the form:

    \[ PE=\int_0^L\!\frac{EA}{2}(a_1+2a_2X_1+3a_3X_1^2)^2 \,\mathrm{d}X_1 - P(a_1L+a_2L^2+a_3L^3)-\int_0^L\!C(a_1X_1^2+a_2X_1^3+a_3X_1^4)\,\mathrm{d}X_1 \]

Clearly, the method generates a large set of equations that are difficult to solve by hand. Using Mathematica we can find the values of the constants a_1, a_2, and a_3 that would minimize the potential energy of the system:
View Mathematica Code:

Clear[u,x,a2,a1,a3,Ee,A,c,P,L]  
u=a3*x^3+a2*x^2+a1*x;
PE=Integrate[(1/2)*Ee*A*D[u,x]^2,{x,0,L}]-P*(u/.x->L)-Integrate[c*x*u,{x,0,L}];
Eq1=D[PE,a2]  
Eq2=D[PE,a1]  
Eq3=D[PE,a3]  
s=Solve[{Eq1==0,Eq2==0,Eq3==0},{a2,a1,a3}]  
u=u/.s[[1]]  
stress=FullSimplify[Ee*D[u,x]]  

Since the assumed trial function is a polynomial of the third degree, which is similar to the exact solution, the Rayleigh Ritz method is able to produce the exact solution for the displacement and the stresses:

    \[\begin{split} u_1 & =-\frac{CX_1^3}{6EA}+\left(\frac{P}{AE}+\frac{CL^2}{2EA}\right)X_1 \\ \sigma_{11}& = E\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=-\frac{CX_1^2}{2A}+\left(\frac{P}{A}+\frac{CL^2}{2A}\right) \end{split} \]

 

Example 2: A Bar under Axial Forces with Varying Cross Sectional Area

Let a bar with a length 2m be aligned with the coordinate axis X_1, and let the width of the bar be equal to 250mm while the height varies linearly such that the height is equal to 500mm when X_1=0 and is equal to 250mm when X_1=2m. Assume that the bar has a constant force of value 200N applied at the end X_1 =2m and that the bar is fixed at the end X_1=0. Find the displacement of the bar by directly solving the differential equation of equilibrium and then using the Rayleigh Ritz method assuming a polynomials of the degrees (0, 1, 2 and 3). Assume the bar to be linear elastic with Young’s modulus E=100,000Pa and assume that the small strain tensor is the appropriate measure of strain. Ignore the effect of Poisson’s ratio. Compare the solution obtained using the Rayleigh Ritz method with the exact solution.
RR2

Solution

Exact Solution:
First, the exact solution that would satisfy the equilibrium equations can be obtained. The equilibrium equation as shown in the beams under axial loading section when E is constant and p=0 is:

    \[ EA\frac{\mathrm{d}^2u_1}{\mathrm{d}X_1^2}+E\frac{\mathrm{d}A}{\mathrm{d}X_1}\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=0 \]

The area A varies linearly with X_1 according to the equation:

    \[ A=0.25\left(0.5-\frac{0.25}{2}X_1\right) \]

The two boundary conditions are:

    \[\begin{split} @X_1=0 & :u_1=0\\ @X_1=L & :\sigma_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=\frac{P}{A}\bigg|_{X_1=L} \end{split} \]

Therefore, the exact solution for the displacement u_1 and the stress \sigma_{11} are:

    \[\begin{split} u_1 & =\frac{8}{125}\left(\ln{4}-\ln\left(4-X_1\right)\right)\\ \sigma_{11} & =\frac{6400}{4-X_1} \end{split} \]

Note that the exact solution for the stress could have been directly obtained by dividing the force P by the variable cross sectional area. The displacement function would then be obtained by integrating the expression \sigma_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1}.
View Mathematica Code

Clear[u,x]  
Ax=25/100(5/10-25/200*x);
Ee=100000;
L=2;
P=200;
s=DSolve[{Ee*u''[x]*Ax+Ee*u'[x]*D[Ax,x]==0,u'[L]==200/Ee/(Ax/.x->L),u[0]==0},u,x]  
u=u[x]/.s[[1]]  
stress=Ee*D[u,x]

The Rayleigh Ritz Method:
The first step in the Rayleigh Ritz method finds the minimizer of the potential energy of the system which can be written as:

    \[ PE=\int\!\overline{U}\,dX - W=\int_0^L\!\frac{\sigma_{11}\varepsilon_{11}}{2}A \,\mathrm{d}X_1 - Pu_1|_{X_1=L} \]

Notice that the potential energy lost by the action of the end force P is equal to the product of P and the displacement u_1 evaluated at X_1= L. Further, we will use the constitutive equation to rewrite the potential energy function in terms of the function u_1:

    \[ PE=\int_0^L\!\frac{EA}{2}\left(\frac{\mathrm{d}u_1}{\mathrm{d}X_1}\right)^2 \,\mathrm{d}X_1 - Pu_1|_{X_1=L} \]

To find an approximate solution, an assumption for the shape or the form of u_1 has to be introduced. Similar to the previous example, a polynomial of the zero degree cannot be used as it will automatically render the displacement u_1=0.
Polynomial of the First Degree:
As a more refined approximation, we can assume that u_1 is a linear function (a polynomial of the first degree):

    \[ u_1 = a_0+a_1 X_1 \]

To satisfy the essential boundary condition @X_1=0:u_1=0 the coefficient a_0 is equal to zero. Therefore:

    \[ u_1 = a_1 X_1 \]

The strain component \varepsilon_{11} can be computed as follows:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=a_1 \]

Substituting into the equation for the potential energy of the system:

    \[ PE=\int_0^L\!\frac{EA}{2}a_1^2 \,\mathrm{d}X_1 - Pa_1L=0.125Ea_1^2\int_0^L\!(0.5-0.125X_1) \,\mathrm{d}X_1-Pa_1L \]

The only variable that can be controlled to minimize the potential energy of the system is a_1. Therefore:

    \[ \frac{\mathrm{d}PE}{\mathrm{d}a_1}=0\Rightarrow 18750a_1=400 \]

Therefore, the best linear function (units of m) according to the Rayleigh Ritz method is:

    \[ u_1=\frac{8}{375}X_1 \]

A linear displacement would produce a constant stress (units of N/m^2):

    \[ \sigma_{11}=\frac{6400}{3} \]

Unfortunately, this solution is far from being accurate, since it does not satisfy the differential equation of equilibrium (Why?). However, within the only possible or allowed linear shapes, the obtained solution is the minimizer of the potential energy.

Polynomial of the Second Degree:
As a more refined approximation, we can assume that u_1 is a polynomial of the second degree:

    \[ u_1 = a_0+a_1 X_1+a_2X_1^2 \]

To satisfy the essential boundary condition @X_1=0:u_1=0 the coefficient a_0 is equal to zero. Therefore:

    \[ u_1 = a_1 X_1+a_2X_1^2 \]

The strain component \varepsilon_{11} can be computed as follows:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=a_1+2a_2X_1 \]

Substituting into the equation for the potential energy of the system:

    \[\begin{split} PE&=\int_0^L\!\frac{EA}{2}(a_1+2a_2X_1)^2 \,\mathrm{d}X_1 - P(a_1L+a_2L^2)\\ & = 9375a_1^2+\frac{100000}{3}a_1a_2+\frac{125000}{3}a_2^2-200(2a_1+4a_2) \end{split} \]

a_1 and a_2 are the variables that can be controlled to minimize the potential energy of the system. Therefore:

    \[ \frac{\partial PE}{\partial a_1}=0\Rightarrow -400+18750a_1+\frac{100000}{3}a_2=0 \]

and

    \[ \frac{\partial PE}{\partial a_2}=0\Rightarrow -800+\frac{100000}{3}a_1+\frac{250000}{3}a_2=0 \]

Solving the above two equations yields:

    \[ a_1=\frac{24}{1625}\qquad a_2=\frac{6}{1625} \]

Therefore, the best second degree polynomial according to the Rayleigh Ritz method is:

    \[ u_1=\frac{24}{1625}X_1+\frac{6}{1625}X_1^2 \]

A parabolic displacement would produce a linear stress:

    \[ \sigma_{11}=\frac{9600}{13}(2+X_1) \]

Polynomial of the Third Degree:
As a more refined approximation, we can assume that u_1 is a polynomial of the third degree:

    \[ u_1 = a_0+a_1 X_1+a_2X_1^2+a_3X_1^3 \]

To satisfy the essential boundary condition @X_1=0:u_1=0 the coefficient a_0 is equal to zero. Therefore:

    \[ u_1 = a_1 X_1+a_2X_1^2+a_3X_1^3 \]

The strain component \varepsilon_{11} can be computed as follows:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=a_1+2a_2X_1+3a_3X_1^2 \]

Substituting into the equation for the potential energy of the system:

    \[\begin{split} PE&=\int_0^L\!\frac{EA}{2}(a_1+2a_2X_1+3a_3X_1^2)^2 \,\mathrm{d}X_1 - P(a_1L+a_2L^2+a_3L^3)\\ & = 9375a_1^2+\frac{100000}{3}a_1a_2+\frac{125000}{3}a_2^2+62500a_1a_3+180000a_2a_3+210000a_3^2\\ &-200(2a_1+4a_2+8a_3) \end{split} \]

a_1, a_2, and a_3 are the variables that can be controlled to minimize the potential energy of the system. Therefore:

    \[ \frac{\partial PE}{\partial a_1}=0\Rightarrow -400+18750a_1+\frac{100000}{3}a_2+62500a_3=0 \]

Also:

    \[ \frac{\partial PE}{\partial a_2}=0\Rightarrow -800+\frac{100000}{3}a_1+\frac{250000}{3}a_2+180000a_3=0 \]

And:

    \[ \frac{\partial PE}{\partial a_3}=0\Rightarrow -1600+62500a_1+180000a_2+420000a_3=0 \]

Solving the above three equations yields:

    \[ a_1=\frac{128}{7875}\qquad a_2=\frac{2}{1575} \qquad a_3=\frac{4}{4725} \]

Therefore, the best third degree polynomial according to the Rayleigh Ritz method is:

    \[ u_1=\frac{128}{7875}X_1+\frac{2}{1575}X_1^2+\frac{4}{4725}X_1^3 \]

A cubic displacement would produce a parabolic stress distribution:

    \[ \sigma_{11}=\frac{3200}{63}(32+5X_1(1+X_1)) \]

View Mathematica Code
Clear[u, u1, u2, u3, u4, stress, stress1, stress2, stress3, x]
(*Exact Solution *)
Ax=25/100(5/10-25/200*x);
Ee=100000;
L=2;
P=200;
s=DSolve[{Ee*u''[x]*Ax+Ee*u'[x]*D[Ax,x]==0,u'[L]==200/Ee/(Ax/.x->L),u[0]==0},u,x]
u=u[x]/.s[[1]];
stress=Ee*D[u,x]
(*First Degree *)
u1=a1*x;
PE=Integrate[1/2*Ee*D[u1,x]^2*Ax,{x,0,L}]-P*u1/.x->L;
Eq1=D[PE,a1];
s=Solve[Eq1==0,a1];
u1=u1/.s[[1]];
stress1=Ee*D[u1,x];
(*Second Degree *)
u2=a1*x+a2*x^2;
PE=Integrate[1/2*Ee*D[u2,x]^2*Ax,{x,0,L}]-P*u2/.x-> L;
Eq1=D[PE,a1];
Eq2=D[PE,a2];
s=Solve[{Eq1==0,Eq2==0},{a1,a2}];
u2=u2/.s[[1]];
stress2=Ee*D[u2,x];
(*Third Degree *)
u3=a1*x+a2*x^2+a3*x^3;
PE=Integrate[1/2*Ee*D[u3,x]^2*Ax,{x,0,L}]-P*u3/.x->L;
Eq1=D[PE,a1];
Eq2=D[PE,a2];
Eq3=D[PE,a3];
s=Solve[{Eq1==0,Eq2==0,Eq3==0},{a1,a2,a3}];
u3=u3/.s[[1]];
stress3=Ee*D[u3,x];
(*comparison*)
Plot[{u, u1, u2, u3}, {x, 0, L}, PlotLegends -> {Style["Exact", Bold, 12], Style["u1", Bold, 12], Style["u2", Bold, 12], Style["u3", Bold, 12]}, AxesLabel -> {"X1 m", "u(m)"}, BaseStyle -> Directive[Bold, 12]]  
Plot[{stress, stress1, stress2, stress3}, {x, 0, L}, AxesLabel -> {"X1 m", "stress(Pa)"}, BaseStyle -> Directive[Bold, 12], PlotLegends -> {Style["Exact", Bold, 12], Style["u1", Bold, 12],   Style["u2", Bold, 12], Style["u3", Bold, 12]}]

The exact solution obtained is a logarithmic function, while the approximations used for the Rayleigh Ritz method are polynomial functions. The Taylor series of a logarithmic function can be represented by a linear combination of polynomial functions, and so, the higher the degree of approximation, the closer the approximation to the exact solution. As seen in the plots below, good accuracy is obtained for a polynomial of the second degree for the displacement function; however, the stresses in that case have higher errors. A better approximation for the stresses is obtained using a polynomial of the third degree. In general, the error in the approximate solution (i.e., the error in u) is less than the error in its derivatives which in this case is the value of the strain component \varepsilon_{11}. This means as well that in order to achieve higher accuracy in the stress \sigma_{11}=E\frac{\mathrm{d}u_1}{\mathrm{d}X_1}, higher polynomials are needed.
RR2Sol

 

Example 3: A Bar under Varying Body Forces with Fixed Ends

Let a bar with a length 2mbe aligned with the coordinate axis X_1, and let the cross sectional area A of the bar be equal to 250\times 250mm^2. Assume that the bar is fixed at both ends X_1= 0 and X_1=2m. Assume that the bar is subjected to a varying body force that is equal to p=5X_1^2 N/m in the direction of the coordinate axis X_1. Find the displacement of the bar by directly solving the differential equation of equilibrium and then using the Rayleigh Ritz method assuming polynomials of the degrees (0, 1, 2, 3, and 4). Assume that the bar material is linear elastic with Young’s modulus E=100,000 Pa and that the small strain tensor is the appropriate measure of strain. Ignore the effect of Poisson’s ratio. Compare the solution obtained using the Rayleigh Ritz method using three- and four-degree polynomials as approximations with the exact solution.
RR3

Solution

Exact Solution:
Similar to the previous examples, the equilibrium equation is:

    \[ \frac{\mathrm{d}^2u_1}{\mathrm{d}X_1^2}=-\frac{p}{EA}=-\frac{5X_1^2}{0.25^2E} \]

Therefore, the solution for u_1 has the form:

    \[ u_1=-\frac{5X_1^4}{12(0.25)^2E}+C_1X_1+C_2 \]

The boundary conditions for this example are:

    \[\begin{split} @X_1=0&:u_1=0\Rightarrow C_2=0\\ @X_1=L&:u_1=0\Rightarrow C_1=\frac{5L^3}{12(0.25)^2E} \end{split} \]

Therefore the exact solutions for the displacement (in m) and the stress (in N/m^2) are:

    \[\begin{split} u_1&=\frac{8X_1-X_1^4}{15000}\\ \sigma_{11}&=\frac{20}{3}(8-4X_1^3) \end{split} \]

The Rayleigh Ritz Method:
The first step in the Rayleigh Ritz method finds the minimizer of the potential energy of the system which can be written as:

    \[ PE=\int\!\overline{U}\,dX - W=\int_0^L\!\frac{\sigma_{11}\varepsilon_{11}}{2}A \,\mathrm{d}X_1 - \int_0^L\!pu_1\,\mathrm{d}X_1 \]

Notice that the potential energy lost by the action of the distributed body forces is an integral since it acts on each point along the beam length. Further, we will use the constitutive equation to rewrite the potential energy function in terms of the function u_1:

    \[ PE=\int_0^L\!\frac{EA}{2}\left(\frac{\mathrm{d}u_1}{\mathrm{d}X_1}\right)^2 \,\mathrm{d}X_1 - \int_0^L\!5X_1^2u_1\,\mathrm{d}X_1 \]

To find an approximate solution, an assumption for the shape or the form of u_1 has to be introduced.
Polynomials of the Zero and First Degrees:
Assuming that that u_1 is a constant function (a polynomial of the zero degree) or a linear function (a polynomial of the first degree) and incorporating the boundary conditions:

    \[\begin{split} @X_1=0&:u_1=0\\ @X_1=L&:u_1=0 \end{split} \]

leads to the trivial solution u_1 = 0, which will automatically be rejected.

Polynomial of the Second Degree:
Let u_1 be a polynomial of the second degree:

    \[ u_1 = a_0+a_1 X_1+a_2X_1^2 \]

The essential boundary condition can be used to find values for the constants a_0 and a_2 as follows:

    \[\begin{split} @X_1=0&:u_1=0\Rightarrow a_0=0\\ @X_1=L&:u_1=0\Rightarrow a_2=-\frac{a_1}{L} \end{split} \]

Therefore, the displacement function has only one parameter a_1 that can be controlled:

    \[ u_1=a_1X_1\left(1-\frac{X_1}{L}\right) \]

The strain component \varepsilon_{11} can be computed as follows:

    \[ \varepsilon_{11}=\frac{\mathrm{d}u_1}{\mathrm{d}X_1}=a_1\left(1-2\frac{X_1}{L}\right) \]

Substituting into the equation for the potential energy of the system:

    \[\begin{split} PE&=\int_0^L\!\frac{EA}{2}\left(a_1\left(1-2\frac{X_1}{L}\right)\right)^2 \,\mathrm{d}X_1 -\int_0^L\!5X_1^2a_1X_1\left(1-\frac{X_1}{L}\right)\,\mathrm{d}X_1\\ & = \frac{6250}{3}a_1^2-4a_1 \end{split} \]

a_1 can be found by minimizing the potential energy of the system. Therefore:

    \[ \frac{\mathrm{d}PE}{\mathrm{d}a_1}=0\Rightarrow \frac{12500}{3}a_1-4=0\Rightarrow a_1=\frac{3}{3125} \]

Therefore, the best second degree polynomial according to the Rayleigh Ritz method is:

    \[ u_1=\frac{3X_1}{3125}\left(1-\frac{X_1}{2}\right) \]

A parabolic displacement would produce a linear stress:

    \[ \sigma_{11}=96(1-X_1) \]

Polynomials of Higher Degrees:
The results from the following Mathematica code show that the third-degree approximation is very close to the actual solution for both displacements and stresses. The polynomial of the
fourth-degree approximation is in fact exact since the exact solution is a polynomial of the fourth degree. Study the Mathematica code carefully to understand how the second boundary conditions of displacement are incorporated. The graph comparing the results are shown below.

View Mathematica Code
Clear[u, x, a1, a2, a3, a4, a5, a0, Ee, A, L, p, u3, u4]  
p = 5 x^2;
s = DSolve[{Ee*u''[x] == -p/A, u[L] == 0, u[0] == 0}, u, x];
u = u[x] /. s[[1]];
stress = Ee*D[u, x];
(*Second Degree*)
u2 = a2*x^2 + a1*x;
sol = Solve[(u2 /. x -> L) == 0, a2];
u2 = u2 /. sol[[1]];
PE = 1/2*Integrate[Ee*A*D[u2, x]^2, {x, 0, L}] - Integrate[p*u2, {x, 0, L}];
Eq1 = D[PE, a1];
sol1 = Solve[{Eq1 == 0}, {a1}];
u2 = FullSimplify[u2 /. sol1[[1]]];
stress2 = FullSimplify[Ee*D[u2, x]];
(*Third Degree*)
u3 = a3*x^3 + a2*x^2 + a1*x;
sol = Solve[(u3 /. x -> L) == 0, a3];
u3 = u3 /. sol[[1]];
PE = 1/2*Integrate[Ee*A*D[u3, x]^2, {x, 0, L}] -  Integrate[p*u3, {x, 0, L}];
Eq1 = D[PE, a1];
Eq2 = D[PE, a2];
sol1 = Solve[{Eq1 == 0, Eq2 == 0}, {a2, a1}];
u3 = FullSimplify[u3 /. sol1[[1]]];
stress3 = FullSimplify[Ee*D[u3, x]];
(*Fourth Degree*)
u4 = a4*x^4 + a3*x^3 + a2*x^2 + a1*x;
sol = Solve[(u4 /. x -> L) == 0, a4];
u4 = u4 /. sol[[1]];
PE = 1/2*Integrate[Ee*A*D[u4, x]^2, {x, 0, L}] -   Integrate[p*u4, {x, 0, L}];
Eq1 = D[PE, a1];
Eq2 = D[PE, a2];
Eq3 = D[PE, a3];
sol1 = Solve[{Eq1 == 0, Eq2 == 0, Eq3 == 0}, {a1, a2, a3}];
u4 = FullSimplify[u4 /. sol1[[1]]];
stress4 = FullSimplify[Ee*D[u4, x]];
A = (25/100)^2;
Ee = 100000;
L = 2;
(*comparison*)
Plot[{u, u2, u3, u4}, {x, 0, L},  PlotLegends -> {Style["Exact", Bold, 12], Style["u2", Bold, 12],    Style["u3", Bold, 12], Style["u4", Bold, 12]},  AxesLabel -> {"X1 m", "u1(m)"}, BaseStyle -> Directive[Bold, 12]]  
Plot[{stress, stress2, stress3, stress4}, {x, 0, L},  PlotLegends -> {Style["Exact", Bold, 12], Style["u2", Bold, 12],    Style["u3", Bold, 12], Style["u4", Bold, 12]},  AxesLabel -> {"X1 m", "Stress(Pa)"},  BaseStyle -> Directive[Bold, 12]]
RR3Sol

 

Euler Bernoulli Beams under Lateral Loading

For plane Euler Bernoulli beams under lateral loading, the unknown displacement function is y. In this case, the potential energy of the system has the following form (See Euler Bernoulli Beam and energy expressions):

    \[ PE=\mbox{Total Strain Energy} - W=\int\!\frac{EI}{2}\left(\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}\right)^2 \,\mathrm{d}X_1 - W \]

where E is Young’s modulus, I is the moment of inertia, and W is the work done by all the external forces acting on the bar (concentrated forces and distributed loads). A few examples to solve for the displacements, rotations, shear forces, and moments for an Euler Bernoulli beam under lateral loading are presented in this section. Similar to the previous section, the Rayleigh Ritz method starts by choosing a form for the displacement function. This form has to satisfy the “essential” boundary conditions, which for this particular type of problems are those boundary conditions imposed on the displacements and the rotations. The following examples show the method for different external loads, and beam configurations.

Example 1: Cantilever Beam

Let a beam with a length L and a constant cross sectional parameters area A and moment of inertia I be aligned with the coordinate axis X_1. Assume that the beam is subjected to a constant force of value P applied at the end X_1=L in the direction shown in the figure. Assume that the bar is fixed (rotation and displacement are equal to zero) at the end X_1=0. Find the displacement, rotation, bending moment, and shearing forces on the beam by directly solving the differential equation of equilibrium and then using the Rayleigh Ritz method assuming polynomials of the degrees (0,1, 2, 3, and 4). Assume the beam to be linear elastic with Young’s modulus E and that the Euler Bernoulli Beam is an appropriate beam approximation.
RR4

Solution

Exact Solution:
First, the exact solution that would satisfy the equilibrium equations can be obtained. The equilibrium equation as shown in the Euler Bernoulli beam section when E and I are constant is:

    \[ EI\frac{\mathrm{d}^4y}{\mathrm{d}X_1^4}=q \]

For the shown cantilever beam, q=0 and therefore the exact solution for the displacement y has the form:

    \[ y=\frac{C_1X_1^3}{6}+\frac{C_2X_1^2}{2}+C_3X_1+C_4 \]

where the constants can be obtained from the boundary conditions:

    \[ \begin{split} & @X_1=0:y=0\Rightarrow C_4=0\\ & @X_1=0:\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0\Rightarrow C_3=0\\ & @X_1=L:M=0\Rightarrow EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=0\Rightarrow C_2=-C_1L\\ & @X_1=L:V=P\Rightarrow EI\frac{\mathrm{d}^3y}{\mathrm{d}X_1^3}=P\Rightarrow C_1=\frac{P}{EI} \end{split} \]

Therefore, the exact solution for the displacement y, the rotation \theta, the moment M, and the shear V are:

    \[\begin{split} y & =\frac{PX_1^3}{6EI}-\frac{PLX_1^2}{2EI}\\ \theta & = \frac{\mathrm{d}y}{\mathrm{d}X_1}=\frac{PX_1^2}{2EI}-\frac{PLX_1}{EI}\\ M & = EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=PX_1-PL\\ V & = EI\frac{\mathrm{d}^3y}{\mathrm{d}X_1^3}=P \end{split} \]

The Rayleigh Ritz Method:
The first step in the Rayleigh Ritz method finds the minimizer of the potential energy of the system which can be written as:

    \[ PE=\int\!\overline{U}\,dX - W=\int_0^L\!\frac{EI}{2}\left(\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}\right)^2 \,\mathrm{d}X_1 - (-P)y|_{X_1=L} \]

Notice that the potential energy lost by the action of the end force P is equal to the product of -P (P is acting downwards and y is assumed upwards) and the displacement y evaluated at X_1= L. To find an approximate solution, an assumption for the shape or the form of y has to be introduced.
Polynomials of the Zero and First Degrees:
If we assume that y is a constant (polynomial of the zero degree) or a linear (polynomial of the first degree) function, and when attempting to satisfy the essential boundary conditions:

    \[\begin{split} @X_1=0 & :y=0\\ @X_1=0 & :\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0 \end{split} \]

this leads to the trivial solution y= 0, which will automatically be rejected.

Polynomial of the Second Degree:
The approximate solution forms for the displacement y, the rotation \theta, the moment M, and the shear force V functions if y is assumed to be a polynomial of the second degree are:

    \[ \begin{split} y & =a_2X_1^2+a_1X_1+a_0\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =2a_2X_1+a_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =2EIa_2\\ V &=0 \end{split} \]

To satisfy the essential boundary conditions, we have:

    \[ \begin{split} @X_1=0 & :y=0\Rightarrow a_0=0\\ @X_1=0 & :\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0\Rightarrow a_1=0 \end{split} \]

Therefore, the solutions forms have only one parameter a_2 that can be controlled:

    \[ \begin{split} y & =a_2X_1^2\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =2a_2X_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =2EIa_2\\ V &=0 \end{split} \]

Substituting into the equation of the potential energy of an Euler Bernoulli beam:

    \[ PE=\frac{EI}{2}\int_0^L\!\left(2a_2\right)^2\,\mathrm{d}X_1+P\left(a_2L^2\right) \]

The unknown coefficient a_2 can be obtained by minimizing the potential energy:

    \[ \frac{\mathrm{d} PE}{\mathrm{d} a_2} =4EILa_2+PL^2=0\Rightarrow a_2=-\frac{PL}{4EI} \]

Therefore, the “best” parabolic approximation for y along with the corresponding solutions for \theta, M, and V are given by:

    \[ \begin{split} y & =-\frac{PLX_1^2}{4EI}\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =-\frac{PLX_1}{2EI}\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =-\frac{PL}{2}\\ V &=EI\frac{\mathrm{d}^3y}{\mathrm{d}X_1^3} =0 \end{split} \]

Note that because of the simplicity of the approximation, the shear force (which is the third derivative) is always equal to zero.

Polynomial of the Third Degree:
The approximate solution forms for the displacement y, the rotation \theta, the moment M, and the shear force V functions if y is assumed to be a polynomial of the third degree are:

    \[ \begin{split} y & =a_3X_1^3+a_2X_1^2+a_1X_1+a_0\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =3a_3X_1^2+2a_2X_1+a_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =6EIa_3X_1+2EIa_2\\ V &=6EIa_3 \end{split} \]

To satisfy the essential boundary conditions, we have:

    \[ \begin{split} @X_1=0 & :y=0\Rightarrow a_0=0\\ @X_1=0 & :\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0\Rightarrow a_1=0 \end{split} \]

Therefore, the solutions forms have only two parameters a_2 and a_3 that can be controlled:

    \[ \begin{split} y & =a_3X_1^3+a_2X_1^2\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =3a_3X_1^2+2a_2X_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =6EIa_3X_1+2EIa_2\\ V &=6EIa_3 \end{split} \]

Substituting into the equation of the potential energy of an Euler Bernoulli beam:

    \[ PE=\frac{EI}{2}\int_0^L\!\left(6a_3X_1+2a_2\right)^2\,\mathrm{d}X_1+P\left(a_3L^3+a_2L^2\right) \]

The unknown coefficients a_2 and a_3 can be obtained by minimizing the potential energy:

    \[\begin{split} \frac{\partial PE}{\partial a_2} & =\frac{EI}{2}(8a_2L+12a_3L^2)+PL^2=0\\ \frac{\partial PE}{\partial a_3} & =\frac{EI}{2}(12a_2L^2+24a_3L^3)+PL^3=0 \end{split} \]

Solving the above two equations yields:

    \[ a_2=-\frac{PL}{2EI} \qquad a_3=\frac{P}{6EI} \]

Therefore, the “best” cubic approximation for y along with the corresponding solutions for \theta, M, and V are given by:

    \[ \begin{split} y & =\frac{PX_1^3}{6EI}-\frac{PLX_1^2}{2EI}\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =\frac{PX_1^2}{2EI}-\frac{PLX_1}{EI}\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =PX_1-PL\\ V &=EI\frac{\mathrm{d}^3y}{\mathrm{d}X_1^3} =P \end{split} \]

The “approximate” solution obtained here is identical to the “exact” solution obtained above. This is because the trial function for y was assumed to be a cubic polynomial and the exact solution is in fact a cubic polynomial. One should note that when the “form” or “shape” of the approximate solution contains the “form” or “shape” of the exact solution, then the Rayleigh Ritz method renders the exact solution! To see this, we are going to try a finer approximation (A polynomial of the fourth degree).

Polynomial of the Fourth Degree:
The approximate solution forms for the displacement y, the rotation \theta, the moment M, and the shear force V functions if y is assumed to be a polynomial of the fourth degree are:

    \[ \begin{split} y & =a_4X_1^4+a_3X_1^3+a_2X_1^2+a_1X_1+a_0\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =4a_4X_1^3+3a_3X_1^2+2a_2X_1+a_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =12EIa_4X_1^2+6EIa_3X_1+2EIa_2\\ V &=24EIa_4X_1+6EIa_3 \end{split} \]

To satisfy the essential boundary conditions, we have:

    \[ \begin{split} @X_1=0 & :y=0\Rightarrow a_0=0\\ @X_1=0 & :\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0\Rightarrow a_1=0 \end{split} \]

Therefore, the solutions forms have the parameters a_2, a_3, and a_4 that can be controlled:

    \[ \begin{split} y & =a_4X_1^4+a_3X_1^3+a_2X_1^2\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =4a_4X_1^3+3a_3X_1^2+2a_2X_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =12EIa_4X_1^2+6EIa_3X_1+2EIa_2\\ V &=24EIa_4X_1+6EIa_3 \end{split} \]

Substituting into the equation of the potential energy of an Euler Bernoulli beam:

    \[ PE=\frac{EI}{2}\int_0^L\!\left(12a_4X_1^2+6a_3X_1+2a_2\right)^2\,\mathrm{d}X_1+P\left(a_4L^4+a_3L^3+a_2L^2\right) \]

The unknown coefficients a_2, a_3, and a_4 can be obtained by minimizing the potential energy:

    \[ \frac{\partial PE}{\partial a_2} =0\qquad \frac{\partial PE}{\partial a_3} =0 \qquad \frac{\partial PE}{\partial a_4} =0 \]

Solving the above three equations yields:

    \[ a_2=-\frac{PL}{2EI} \qquad a_3=\frac{P}{6EI} \qquad a_4=0 \]

Therefore, assuming a fourth order polynomial gives the same solution as a third order polynomial!

View Mathematica Code
(*Exact Solution*)  
Clear[a,x,V,M,P,y,L,EI]  
M=EI*D[y[x],{x,2}];
V=EI*D[y[x],{x,3}];
th=D[y[x],x];
q=0;
M1=M/.x->0;
M2=M/.x->L;
V1=V/.x->0;
V2=V/.x->L;
y1=y[x]/.x->0;
y2=y[x]/.x->L;
th1=th/.x->0;
th2=th/.x->L;
s=DSolve[{EI*y''''[x]==q,y1==0,th1==0,V2==P,M2==0},y,x];
y=y[x]/.s[[1]]  
th=D[y,x]  
M=EI*D[y,{x,2}]  
V=EI*D[y,{x,3}]  
(*Second Degree*)
Clear[y,x,M,V,th,EI,PE]  
y=a2*x^2  
th=D[y,x];
M=EI*D[y,{x,2}];
V=EI*D[y,{x,3}];
PE=EI/2*Integrate[D[y,{x,2}]^2,{x,0,L}]+P*(y/.x->L);
Eq2=D[PE,a2];
Sol=Solve[{Eq2==0},{a2}]  
y=y/.Sol[[1]]  
th/.Sol[[1]]  
FullSimplify[M/.Sol[[1]]]  
FullSimplify[V/.Sol[[1]]]  
(*Third Degree*)
Clear[y,x,M,V,th,EI,PE]  
y=a3*x^3+a2*x^2  
th=D[y,x];
M=EI*D[y,{x,2}];
V=EI*D[y,{x,3}];
PE=EI/2*Integrate[D[y,{x,2}]^2,{x,0,L}]+P*(y/.x->L);
Eq1=D[PE,a3];
Eq2=D[PE,a2];
Sol=Solve[{Eq1==0,Eq2==0},{a3,a2}]  
y=y/.Sol[[1]]  
th/.Sol[[1]]  
FullSimplify[M/.Sol[[1]]]  
FullSimplify[V/.Sol[[1]]]  
(*Fourth Degree*)
Clear[y,x,M,V,th,EI,PE]  
y=a4*x^4+a3*x^3+a2*x^2  
th=D[y,x];
M=EI*D[y,{x,2}];
V=EI*D[y,{x,3}];
PE=EI/2*Integrate[D[y,{x,2}]^2,{x,0,L}]+P*(y/.x->L);
Eq1=D[PE,a3];
Eq2=D[PE,a2];
Eq3=D[PE,a4];
Sol=Solve[{Eq1==0,Eq2==0,Eq3==0},{a3,a2,a4}]  
y=y/.Sol[[1]]  
th/.Sol[[1]]  
FullSimplify[M/.Sol[[1]]]  
FullSimplify[V/.Sol[[1]]]  

 

Example 2: Simple Beam

Let a beam with a length L=10m and a constant moment of inertia I=4\times 10^8 mm^4 be aligned with the coordinate axis X_1. Assume that the beam is subjected to a distributed load q=-25kN/m.(Notice that the negative sign was added since q is positive when it is in the direction of positive X_2). Assume that the beam is hinged at the ends X_1= 0 and X_1=L. Find the displacement, rotation, bending moment, and shearing forces on the beam by directly solving the differential equation of equilibrium and then using the Rayleigh Ritz method assuming polynomials of the degrees (0, 1, 2, 3 and 4). Assume the beam to be linear elastic with Young’s modulus E= 200 GPa and that the Euler Bernoulli Beam is an appropriate beam approximation. Compare the Rayleigh Ritz method with the exact solution.
RR5

Solution

Exact Solution:
First, the exact solution that would satisfy the equilibrium equations can be obtained. The equilibrium equation as shown in the Euler Bernoulli beam section when E and I are constant is:

    \[ EI\frac{\mathrm{d}^4y}{\mathrm{d}X_1^4}=q \]

For the shown simple beam, q=-25 kN/m and therefore the exact solution for the displacement y has the form:

    \[ y=q\frac{X_1^4}{24}+\frac{C_1X_1^3}{6}+\frac{C_2X_1^2}{2}+C_3X_1+C_4 \]

where the constants can be obtained from the boundary conditions:

    \[ \begin{split} & @X_1=0:y=0\Rightarrow C_4=0\\ & @X_1=0:M=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=0\Rightarrow C_2=0\\ & @X_1=L:y=0\Rightarrow q\frac{L^4}{24}+\frac{C_1L^3}{6}+C_3L=0\\ & @X_1=L:M=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=0\Rightarrow q\frac{L^2}{2}+C_1L=0\\ \end{split} \]

therefore,

    \[ C_1=-\frac{qL}{2}\qquad C_3=\frac{qL^3}{24} \]

Therefore, the exact solution for the displacement y, the rotation \theta, the moment M, and the shear V are:

    \[\begin{split} y & =\frac{qX_1^4}{24EI}-\frac{qLX_1^3}{12EI}+\frac{qL^3X_1}{24EI}\\ \theta & = \frac{\mathrm{d}y}{\mathrm{d}X_1}=\frac{qX_1^3}{6EI}-\frac{qLX_1^2}{4EI}+\frac{qL^3}{24EI}\\ M & = EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=\frac{qX_1^2}{2}-\frac{qLX_1}{2}\\ V & = EI\frac{\mathrm{d}^3y}{\mathrm{d}X_1^3}=qX_1-\frac{qL}{2} \end{split} \]

The Rayleigh Ritz Method:
The first step in the Rayleigh Ritz method finds the minimizer of the potential energy of the system which can be written as:

    \[ PE=\int\!\overline{U}\,dX - W=\int_0^L\!\frac{EI}{2}\left(\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}\right)^2 \,\mathrm{d}X_1 - \int_0^L \! qy\,\mathrm{d}X_1 \]

To find an approximate solution, an assumption for the shape or the form of y has to be introduced.
Polynomials of the Zero and First Degrees:
If we assume that y is a constant (polynomial of the zero degree) or a linear (polynomial of the first degree) function, and when attempting to satisfy the essential boundary conditions:

    \[\begin{split} @X_1=0 & :y=0\\ @X_1=L & :y=0 \end{split} \]

this leads to the trivial solution y= 0, which will automatically be rejected.

Polynomial of Higher Degrees:
Similar to the previous examples, the coefficients are obtained by first satisfying the essential boundary conditions. Then, the remaining coefficients are obtained by minimizing the potential energy of the system. As shown in the plots below, the second and third degree trial functions give the same result, while the fourth and fifth degree trial functions give the same result. The solution is shown to give the accurate result when a polynomial of the fourth or higher degree is used.

RR5example

View Mathematica Code
(*Exact Solution*)  
Clear[q,EI,y,x,c1,c2,c3,c4,L]  
y=(q*x^4/24+c1*x^3/6+c2*x^2/2+c3*x+c4)/EI;
Eq1=y/.x-> 0;
Eq2=y/.x-> L;
Eq3=EI*D[y,{x,2}]/.x-> 0;
Eq4=EI*D[y,{x,2}]/.x-> L;
a=Solve[{Eq1==0,Eq2==0,Eq3==0,Eq4==0},{c1,c2,c3,c4}];
y=y/.a[[1]];
theta=D[y,x]/.a[[1]];
M=FullSimplify[EI*D[y,{x,2}]/.a[[1]]];
V=FullSimplify[EI*D[y,{x,3}]/.a[[1]]];
(*Second Degree*)
y2=a2*x^2+a1*x+a0;
s=Solve[{(y2/.x-> 0)==0,(y2/.x->L)==0},{a0,a1}];
y2=y2/.s[[1]];
theta2=D[y2,x];
M2=EI*D[y2,{x,2}];
V2=EI*D[y2,{x,3}];
PE=EI/2*Integrate[(D[y2,{x,2}])^2,{x,0,L}]-Integrate[q*y2,{x,0,L}];
Eq1=D[PE,a2];
a=Solve[{Eq1==0},{a2}];
y2=y2/.a[[1]];
theta2=D[y2,x]/.a[[1]];
M2=FullSimplify[EI*D[y2,{x,2}]/.a[[1]]];
V2=FullSimplify[EI*D[y2,{x,3}]/.a[[1]]];
(*Third Degree*)
y3=a3*x^3+a2*x^2+a1*x+a0;
s=Solve[{(y3/.x-> 0)==0,(y3/.x->L)==0},{a0,a1}];
y3=y3/.s[[1]];
theta3=D[y3,x];
M3=EI*D[y3,{x,2}];
V3=EI*D[y3,{x,3}];
PE=EI/2*Integrate[(D[y3,{x,2}])^2,{x,0,L}]-Integrate[q*y3,{x,0,L}];
Eq1=D[PE,a2];
Eq2=D[PE,a3];
a=Solve[{Eq1==0,Eq2==0},{a2,a3}];
y3=y3/.a[[1]];
theta3=D[y3,x]/.a[[1]];
M3=FullSimplify[EI*D[y3,{x,2}]/.a[[1]]];
V3=FullSimplify[EI*D[y3,{x,3}]/.a[[1]]];
(*fourth Degree*)
y4=a4*x^4+a3*x^3+a2*x^2+a1*x+a0;
s=Solve[{(y4/.x-> 0)==0,(y4/.x->L)==0},{a0,a1}];
y4=y4/.s[[1]];
theta4=D[y4,x];
M4=EI*D[y4,{x,2}];
V4=EI*D[y4,{x,3}];
PE=EI/2*Integrate[(D[y4,{x,2}])^2,{x,0,L}]-Integrate[q*y4,{x,0,L}];
Eq1=D[PE,a2];
Eq2=D[PE,a3];
Eq3=D[PE,a4];
a=Solve[{Eq1==0,Eq2==0,Eq3==0},{a2,a3,a4}];
y4=y4/.a[[1]];
theta4=D[y4,x]/.a[[1]];
M4=FullSimplify[EI*D[y4,{x,2}]/.a[[1]]];
V4=FullSimplify[EI*D[y4,{x,3}]/.a[[1]]];
(*fifth Degree*)
y5=a5*x^5+a4*x^4+a3*x^3+a2*x^2+a1*x+a0;
s=Solve[{(y5/.x-> 0)==0,(y5/.x->L)==0},{a0,a1}];
y5=y5/.s[[1]];
theta4=D[y5,x];
M5=EI*D[y5,{x,2}];
V5=EI*D[y5,{x,3}];
PE=EI/2*Integrate[(D[y5,{x,2}])^2,{x,0,L}]-Integrate[q*y5,{x,0,L}];
Eq1=D[PE,a2];
Eq2=D[PE,a3];
Eq3=D[PE,a4];
Eq4=D[PE,a5];
a=Solve[{Eq1==0,Eq2==0,Eq3==0,Eq4==0},{a2,a3,a4,a5}];
y5=y5/.a[[1]];
theta5=D[y5,x]/.a[[1]];
M5=FullSimplify[EI*D[y5,{x,2}]/.a[[1]]];
V5=FullSimplify[EI*D[y5,{x,3}]/.a[[1]]];
EI = 200000*1000000*4*100000000/1000000000000;
q = -25*1000;
L = 10;
Plot[{y,y2,y4},{x,0,L},AxesLabel-> {"x (m)","y (m)"},PlotLegends->{Style["Exact y",Bold,12],Style["y2",Bold,12],Style["y4",Bold,12]}]  
Plot[{theta,theta2,theta4},{x,0,L},AxesLabel-> {"x (m)","theta"},PlotLegends-> {Style["Exact theta",Bold,12],Style["theta 2",Bold,12],Style["theta 4",Bold,12]}]  
Plot[{M,M2,M4},{x,0,L},AxesLabel-> {"x (m)","M (N.m.)"},PlotLegends-> {Style["Exact M",Bold,12],Style["M2",Bold,12],Style["M4",Bold,12]}]  
Plot[{V,V2,V4},{x,0,L},AxesLabel-> {"x (m)","V (N)"},PlotLegends->{Style["Exact V",Bold,12],Style["V 2",Bold,12],Style["V4",Bold,12]}]  

Example 3: Cantilever Beam with a Varying Cross Sectional Area

Let a beam with a length L=8m and with a constant width b=0.25m and a linearly varying height (@X_1=0, h=0.5m and @X_1=8, h=0.25m) be aligned with the coordinate axis X_1. Assume that the beam has a fixed left end and had a shearing load acting downwards of 10kN on its right end. Find the displacement, bending moment, and shearing force on the beam by directly solving the differential equation of equilibrium and then using the Rayleigh Ritz method assuming polynomials of the degrees (0, 1, 2, 3, 4, 5, 6, 7). Assume the beam to be linear elastic with Young’s modulus E= 20 GPa and that the Euler Bernoulli Beam is an appropriate beam approximation. Compare the Rayleigh Ritz method with the exact solution.
RR new

Solution

Exact Solution:
First, the exact solution that would satisfy the equilibrium equations can be obtained. The equilibrium equation as shown in the Euler Bernoulli beam section apply only when E and I are constant. Since I is variable, the differential equation has the form:

    \[ \frac{\mathrm{d}^2\left(EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}\right)}{\mathrm{d}X_1^2}=q \]

For the cantilever beam, q=0. For consistency in the units, in all the following calculations E=20,000,000,000 Pa and P=10,000 N. The moment of inertia is a function of X_1 as follows:

    \[ I=\frac{bt^3}{12}=\frac{0.25\left(0.5-0.25\frac{X_1}{8}\right)^3}{12}=\frac{(16-X_1)^3}{1572864} \]

The exact solution can be obtained using Mathematica. The following are the four boundary conditions required:

    \[ \begin{split} & @X_1=0:y=0\\ & @X_1=0:\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0\\ & @X_1=L:V=\frac{\mathrm{d}M}{\mathrm{d}X_1}=\frac{\mathrm{d}\left(EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}\right)}{\mathrm{d}X_1}=10000\\ & @X_1=L:M=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}=0\\ \end{split} \]

Using Mathematica, the above differential equation can be solved to yield the following exact solution for y:

    \[ y=\frac{34.8872+(0.036864X_1-2.96688)X_1+(0.786432X_1-12.5829)\ln[16-X_1]}{X_1-16} \]

The associated moment in N.m. and shearing force in N. are given by:

    \[ M=-80000+10000X_1\qquad V=10000 \]

Obviously, as the beam is statically determinate, the equations for the moment and shear are very simple. However, the variable cross section leads to a complicated form for the displacement y. Note that traditionally, in structural engineering applications, a constant cross section with the intermediate height of two thirds of the maximum height can be assumed to simplify the equations.

The Rayleigh Ritz Method:
The first step in the Rayleigh Ritz method finds the minimizer of the potential energy of the system which can be written as:

    \[ PE=\int\!\overline{U}\,dX - W=\int_0^L\!\frac{EI}{2}\left(\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2}\right)^2 \,\mathrm{d}X_1 - (-P)y|_{X_1=L} \]

Notice that the potential energy lost by the action of the end force P is equal to the product of -P (P is acting downwards and y is assumed upwards) and the displacement y evaluated at X_1= L. To find an approximate solution, an assumption for the shape or the form of y has to be introduced.
Polynomials of the Zero and First Degrees:
If we assume that y is a constant (polynomial of the zero degree) or a linear (polynomial of the first degree) function, and when attempting to satisfy the essential boundary conditions:

    \[\begin{split} @X_1=0 & :y=0\\ @X_1=0 & :\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0 \end{split} \]

this leads to the trivial solution y= 0, which will automatically be rejected.

Polynomial of the Second Degree:
The approximate solution forms for the displacement y, the rotation \theta, the moment M, and the shear force V functions if y is assumed to be a polynomial of the second degree are:

    \[ \begin{split} y & =a_2X_1^2+a_1X_1+a_0\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =2a_2X_1+a_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =2EIa_2=-\frac{9765625}{384}a_2(-16+X_1)^3\\ V &=\frac{\mathrm{d}M}{\mathrm{d}X_1}=-\frac{9765625}{128}a_2(-16+X_1)^2 \end{split} \]

To satisfy the essential boundary conditions, we have:

    \[ \begin{split} @X_1=0 & :y=0\Rightarrow a_0=0\\ @X_1=0 & :\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0\Rightarrow a_1=0 \end{split} \]

Therefore, the solutions forms have only one parameter a_2 that can be controlled:

    \[ \begin{split} y & =a_2X_1^2\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =2a_2X_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =2EIa_2=-\frac{9765625}{384}a_2(-16+X_1)^3\\ V &=\frac{\mathrm{d}M}{\mathrm{d}X_1}=-\frac{9765625}{128}a_2(-16+X_1)^2 \end{split} \]

Substituting into the equation of the potential energy of an Euler Bernoulli beam:

    \[\begin{split} PE&=\int_0^L\!\frac{EI}{2} \left(2a_2\right)^2\,\mathrm{d}X_1+P\left(a_2L^2\right)=\int_0^L\!\frac{9765625a_2^2(16-X_1)^3}{384}\,\mathrm{d}X_1+640000a_2\\ &=640000a_2+390625000a_2^2 \end{split} \]

Note that I is a function of X_1 and so it stays inside the integral. The unknown coefficient a_2 can be obtained by minimizing the potential energy:

    \[ \frac{\mathrm{d} PE}{\mathrm{d} a_2} =0\Rightarrow a_2=\frac{-64}{78125} \]

Therefore, the “best” parabolic approximation for y along with the corresponding solutions for M and V are given by:

    \[ \begin{split} y & =\frac{-64}{78125}X_1^2=-0.0008192X_1^2\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =\frac{125}{6}(X_1-16)^3\\ &=-85333.3 + 16000 X_1 - 1000 X_1^2 + 20.8333 X_1^3\\ V &=\frac{\mathrm{d}M}{\mathrm{d}X_1} =\frac{125}{2}(X_1-16)^2\\ &=16000 - 2000 X_1 + 62.5 X_1^2 \end{split} \]

Polynomial of the Third Degree:
The approximate solution forms for the displacement y, the rotation \theta, the moment M, and the shear force V functions if y is assumed to be a polynomial of the third degree are:

    \[ \begin{split} y & =a_3X_1^3+a_2X_1^2+a_1X_1+a_0\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =3a_3X_1^2+2a_2X_1+a_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =6EIa_3X_1+2EIa_2\\ V &=\frac{\mathrm{d}M}{\mathrm{d}X_1} \end{split} \]

To satisfy the essential boundary conditions, we have:

    \[ \begin{split} @X_1=0 & :y=0\Rightarrow a_0=0\\ @X_1=0 & :\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0\Rightarrow a_1=0 \end{split} \]

Therefore, the solutions forms have only two parameters a_2 and a_3 that can be controlled:

    \[ \begin{split} y & =a_3X_1^3+a_2X_1^2\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =3a_3X_1^2+2a_2X_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =6EIa_3X_1+2EIa_2=-\frac{9765625}{384}(-16+X_1)^3(a_2+3a_3X_1)\\ V &=-\frac{9765625}{128}(-16+X_1)^2(a_2+4a_3(X_1-4)) \end{split} \]

Substituting into the equation of the potential energy of an Euler Bernoulli beam:

    \[ PE=\int_0^L\!\frac{EI}{2}\left(6a_3X_1+2a_2\right)^2\,\mathrm{d}X_1+P\left(a_3L^3+a_2L^2\right) \]

The unknown coefficients a_2 and a_3 can be obtained by minimizing the potential energy:

    \[\begin{split} \frac{\partial PE}{\partial a_2} & =640000+781250000a_2+6500000000a_3=0\\ \frac{\partial PE}{\partial a_3} & =5120000+6500000000a_2+84000000000a_3=0 \end{split} \]

Solving the above two equations yields:

    \[ a_2=-\frac{512}{584375} \qquad a_3=\frac{4}{584375} \]

Therefore, the “best” cubic approximation for y along with the corresponding solutions for \theta, M, and V are given by:

    \[ \begin{split} y & =\frac{4X_1^3}{584375}-\frac{512X_1^2}{584375}\\ &=-0.00087615X_1^2+6.84492(10)^{-6}X_1^3\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =-\frac{3125(-16+X_1)^3(-128+3X_1)}{17952}\\ &=-91265.6 + 19251.3 X_1 - 1470.59 X_1^2 + 47.3485 X_1^3 - 0.522226 X_1^4\\ V &=\frac{\mathrm{d}M}{\mathrm{d}X_1} =-\frac{3125(-36+X_1)(-16+X_1)^2}{1496}\\ &=19251.3 - 2941.18 X_1 + 142.045 X_1^2 - 2.0889 X_1^3 \end{split} \]

Polynomials of the Fourth Degree:
The approximate solution forms for the displacement y, the rotation \theta, the moment M, and the shear force V functions if y is assumed to be a polynomial of the fourth degree are:

    \[ \begin{split} y & =a_4X_1^4+a_3X_1^3+a_2X_1^2+a_1X_1+a_0\\ \theta &=\frac{\mathrm{d}y}{\mathrm{d}X_1} =4a_4X_1^3+3a_3X_1^2+2a_2X_1+a_1\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =12EIa_4X_1^2+6EIa_3X_1+2EIa_2\\ V &=\frac{\mathrm{d}M}{\mathrm{d}X_1} \end{split} \]

To satisfy the essential boundary conditions, we have:

    \[ \begin{split} @X_1=0 & :y=0\Rightarrow a_0=0\\ @X_1=0 & :\theta=\frac{\mathrm{d}y}{\mathrm{d}X_1}=0\Rightarrow a_1=0 \end{split} \]

Therefore, the solutions forms have the parameters a_2, a_3, and a_4 that can be controlled:

    \[ \begin{split} y & =a_4X_1^4+a_3X_1^3+a_2X_1^2\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =12EIa_4X_1^2+6EIa_3X_1+2EIa_2\\ V &=\frac{\mathrm{d}M}{\mathrm{d}X_1} \end{split} \]

Substituting into the equation of the potential energy of an Euler Bernoulli beam:

    \[ PE=\int_0^L\!\frac{EI}{2}\left(12a_4X_1^2+6a_3X_1+2a_2\right)^2\,\mathrm{d}X_1+P\left(a_4L^4+a_3L^3+a_2L^2\right) \]

The unknown coefficients a_2, a_3, and a_4 can be obtained by minimizing the potential energy:

    \[ \frac{\partial PE}{\partial a_2} =0\qquad \frac{\partial PE}{\partial a_3} =0 \qquad \frac{\partial PE}{\partial a_4} =0 \]

Solving the above three equations yields:

    \[ a_2=-0.000704051 \qquad a_3=-0.0000484584 \qquad a_4=-4.01821(10)^{-6} \]

Therefore, assuming a fourth order polynomial gives the following set of solutions:

    \[ \begin{split} y & =-0.000704051X_1^2-0.0000484584X_1^3+4.01821(10)^{-6}X_1^4\\ M &=EI\frac{\mathrm{d}^2y}{\mathrm{d}X_1^2} =-73338.7 - 1392.24 X_1 + 4491.3 X_1^2 - 630.438 X_1^3 + 33.1273 X_1^4 - 0.61313 X_1^5\\ V &=\frac{\mathrm{d}M}{\mathrm{d}X_1}=-1392.24 + 8982.6 X_1 - 1891.32 X_1^2 + 132.509 X_1^3 - 3.06565 X_1^4 \end{split} \]

Approximate solutions using higher order polynomials can be obtained similarly. The following plots show the convergence of the polynomial functions to the exact solution up to a polynomial of the sixth degree (Mathematica code is given below). In the first plot, it is clear that a polynomial of the second degree is capable of capturing the displacement quite accurately. However, the corresponding approximate solutions for the moment and shear are not as accurate. The accuracy of the approximate solution for the bending moment is higher than that for the shearing force and that is because the moment is related to the second derivative of y while the shear is related to the third derivative of y. It takes up to a polynomial of the seventh degree for y for the approximate shearing force to have a flat distribution similar to the exact solution.

RR new 2
View Mathematica Code

Clear[y, b, L, t, Ii, Ee, P]
b = 1/4;
L = 8;
t = (1/2 - 1/4*X1/L);
Ii = b*t^3/12;
Ee = 20*10^9;
P = 10000;
a = DSolve[{D[Ee*Ii*y''[X1], {X1, 2}] == 0, y[0] == 0, y'[0] == 0, y''[L] == 0, (D[Ee*Ii, X1] /. X1 -> L)*y''[L] + (Ee*Ii /. X1 -> L)*y'''[L] == P}, y[X1], X1, Assumptions -> 0 <= X1 <= L]
y = FullSimplify[Chop[N[y[X1] /. a[[1]]]]]
th = D[y, X1];
M = FullSimplify[Ee*Ii*D[y, {X1, 2}]]
V = FullSimplify[D[M, X1]]

Plot[y, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact"}]
Plot[M, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact"}]
Plot[V, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact"}]

Print["Second Order Polynomial"]
y2 = a2*X1^2;
d2y = D[y2, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y2 /. X1 -> L;
sol = Solve[D[PE, a2] == 0, a2]
y2 = y2 /. sol[[1]]
M2 = FullSimplify[Ee*Ii*D[y2, {X1, 2}]]
V2 = FullSimplify[D[M2, X1]]
Plot[{y, y2}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y2"}]
Plot[{M, M2}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M2"}]
Plot[{V, V2}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V2"}]

Print["Third Order Polynomial"]
y3 = a2*X1^2 + a3*X1^3;
d2y = D[y3, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y3 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0}, {a2, a3}]
y3 = y3 /. sol[[1]]
M3 = FullSimplify[Ee*Ii*D[y3, {X1, 2}]]
V3 = FullSimplify[D[M3, X1]]
Plot[{y, y3}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y3"}]
Plot[{M, M3}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M3"}]
Plot[{V, V3}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"},  PlotLegends -> {"V exact", "V3"}]

Print["Fourth Order Polynomial"]
y4 = a2*X1^2 + a3*X1^3 + a4*X1^4;
d2y = D[y4, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y4 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0, D[PE, a4] == 0}, {a2, a3, a4}]
y4 = y4 /. sol[[1]]
M4 = FullSimplify[Ee*Ii*D[y4, {X1, 2}]]
V4 = FullSimplify[D[M4, X1]]
Plot[{y, y4}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y4"}]
Plot[{M, M4}, {X1, 0, L}, PlotLabel -> "Moment M", AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M4"}]
Plot[{V, V4}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V4"}]

Print["Fifth Order Polynomial"]
y5 = a2*X1^2 + a3*X1^3 + a4*X1^4 + a5*X1^5;
d2y = D[y5, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y5 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0, D[PE, a4] == 0, D[PE, a5] == 0}, {a2, a3, a4, a5}]
y5 = y5 /. sol[[1]]
M5 = FullSimplify[Ee*Ii*D[y5, {X1, 2}]]
V5 = FullSimplify[D[M5, X1]]
Plot[{y, y5}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y5"}]
Plot[{M, M5}, {X1, 0, L}, PlotLabel -> "Moment M",  AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M5"}]
Plot[{V, V5}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}}, PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V5"}]

Print["Sixth Order Polynomial"]
y6 = a2*X1^2 + a3*X1^3 + a4*X1^4 + a5*X1^5 + a6*X1^6;
d2y = D[y6, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y6 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0, D[PE, a4] == 0, D[PE, a5] == 0, D[PE, a6] == 0}, {a2, a3, a4, a5, a6}]
y6 = y6 /. sol[[1]]
M6 = FullSimplify[Ee*Ii*D[y6, {X1, 2}]]
V6 = FullSimplify[D[M6, X1]]
Plot[{y, y6}, {X1, 0, L}, PlotLabel -> "Displacement y",  AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y6"}]
Plot[{M, M6}, {X1, 0, L}, PlotLabel -> "Moment M",  AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M6"}]
Plot[{V, V6}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}},  PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"},  PlotLegends -> {"V exact", "V6"}]

Print["Seventh Order Polynomial"]
y7 = a2*X1^2 + a3*X1^3 + a4*X1^4 + a5*X1^5 + a6*X1^6 + a7*X1^7;
d2y = D[y7, {X1, 2}];
PE = Integrate[Ee*Ii/2*d2y^2, {X1, 0, L}] - (-P)*y7 /. X1 -> L;
sol = Solve[{D[PE, a2] == 0, D[PE, a3] == 0, D[PE, a4] == 0, D[PE, a5] == 0, D[PE, a6] == 0, D[PE, a7] == 0}, {a2, a3, a4, a5,a6, a7}]
y7 = y7 /. sol[[1]]
M7 = FullSimplify[Ee*Ii*D[y7, {X1, 2}]]
V7 = FullSimplify[D[M7, X1]]
Plot[{y, y7}, {X1, 0, L}, PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"}, PlotLegends -> {"y exact", "y7"}]
Plot[{M, M7}, {X1, 0, L}, PlotLabel -> "Moment M",  AxesLabel -> {"X1 (m)", "M (N.m)"}, PlotLegends -> {"M exact", "M7"}]
Plot[{V, V7}, {X1, 0, L}, PlotRange -> {{0, L}, {0, 2 P}},  PlotLabel -> "Shear V", AxesLabel -> {"X1 (m)", "V (N)"}, PlotLegends -> {"V exact", "V7"}]


Plot[{y, y2, y3, y4, y5, y6, y7}, {X1, 0, L},  PlotLabel -> "Displacement y", AxesLabel -> {"X1 (m)", "y (m)"},  PlotLegends -> {"y exact", "y2", "y3", "y4", "y5", "y6", "y7"}]
Plot[{M, M2, M3, M4, M5, M6, M7}, {X1, 0, L}, PlotLabel -> "Moment M",  AxesLabel -> {"X1 (m)", "M (N.m)"},  PlotLegends -> {"M exact", "M2", "M3", "M4", "M5", "M6", "M7"}]
Plot[{V,V2,V3,V4,V5,V6,V7}, {X1, 0, L}, PlotRange ->{{0, L}, {0, 2 P}},PlotLabel->"Shear V",AxesLabel -> {"X1 (m)", "V (N)"},PlotLegends -> {"V exact", "V2", "V3", "V4", "V5", "V6", "V7"}]

 

Problems

  1. The exact displacement in meters of the shown Euler Bernoulli beam follows the function:

        \[ y=\begin{cases} 0.0003375X_1^2(3X_1-20) & 0\leq X_1\leq 5\\ -0.0003125(X_1-8)^2(7X_1-20) & 5< X_1\leq 8 \end{cases} \]

    The beam’s Young’s modulus and moment of inertia are E=20,000MPa and I=0.00260417m^4.
    RRp1

    1. Find the strain energy stored in the beam (Answer: 21093.8 N.m.).
    2. Use the Rayleigh Ritz method to find approximate solutions for the displacement across the whole beam. Assume the approximations to be polynomials of the zero, first, second, third, fourth, fifth, and sixth degrees.
    3. Plot the exact solution versus the approximate solutions for y, \theta, M, and V. Comment on the results.
  2. For the following two linear elastic isotropic small deformations beam structures:
    RRp2

    1. Find the exact solution of the displacement.
    2. Find an approximate solution using the Rayleigh Ritz method with polynomial trial functions of zero, first, second, third, and fourth degrees satisfying the essential boundary conditions. Comment on the difference between the approximate and the exact solutions.
  3. An Euler Bernoulli beam of length L with both ends fixed is loaded by a concentrated force P at the centre of the beam. Assume that the beam is made out of a linear elastic isotropic material and that the deformations are small. If E is Young’s modulus and I is the moment of inertia, then find:
    1. the strain energy in terms of the displacement function y.
    2. The expression for the potential energy of the system.
    3. Show that the following series satisfies the essential boundary conditions

          \[ w=\sum_{n=2}^\infty A_n\left(\cos{\frac{n\pi x}{L}}-1\right) \]

      where n=2, 4, 6, \cdots and A_n are constants

    4. Take one term and two terms from the series above and determine the constants using the Rayleigh Ritz method.
    5. Find the exact solution and compare with the approximate solutions obtained using the Rayleigh Ritz method.
  4. For the shown two beam structures
    RRp3

    1. Find the exact solution for the displacement.
    2. Use the Rayleigh Ritz method with polynomials of the zero, one, two, three, four, five, and six degrees to find an approximate solution for the displacement. For the beam under axial loading compare the stresses with the exact solution. For the beam under lateral loading, compare the rotation, moment, and shear with the corresponding exact solutions.
  5. The shown beam under axial loading is subjected to the varying body force given per unit volume. Assume that the area varies linearly across the length of the structure.
    RRp4

    1. Find the exact displacement by solving the differential equation of equilibrium for bars under axial loading.
    2. Use the Rayleigh Ritz method to find a polynomial approximation of the zero, first, second, and third degrees for the displacement function.
  6. The shown plate has a thickness t = 10mm. The plate is in a plane stress condition. The left boundary where X_1= 0 has zero vertical and horizontal displacements. The material is linear elastic with small deformations and Young’s modulus E and Poisson’s ratio \nu are equal to 200 GPa and 0.3, respectively. The right boundary is subjected to a constant normal tensile stress of 20MPa while the bottom boundary is subjected to a normal compressive stress that varies linearly with a maximum value of 20MPa as shown in the figure below. Using the Rayleigh Ritz method and assuming that the displacement of the plate can be approximated using the function:

        \[ u=\left(\begin{array}{c}a_0+a_1X_1+a_2X_2+a_3X_1X_2\\b_0+b_1X_1+b_2X_2+b_3X_1X_2\end{array}\right) \]

    RRp5

    1. Utilize the essential boundary conditions to reduce the number of unknown constants in the displacement function.
    2. Find the infinitesimal strain tensor in terms of the remaining unknown constants.
    3. Write the expression of the potential energy of the system in terms of the remaining unknown constants.
    4. Find the unknown constants using the Rayleigh Ritz method.
    5. Comment on how accurate the above displacement function is in approximating the exact displacement of the plate.
  7. Consider a uniaxial bar of length L, Young’s modulus E, and variable cross sectional area A=A_0\left(1+\frac{x}{L}\right). The bar is fixed at one end and free at the other end. The free end is subjected to a constant force F. Using the Rayleigh Ritz method, and assuming a trial parabolic displacement function u=a_0+a_1x+a_2x^2, determine the approximate displacement function of the bar. Compare the results with the exact solution given by:

        \[ u=\frac{FL}{EA_0}\ln\left(1+\frac{x}{L}\right) \]

  8. Two bars joined at x=L are subjected to a force P at the joint. The cross sectional area of the larger bar is twice that of the smaller bar. Use the Rayleigh Ritz method with a trial displacement function of u=a_0+a_1x+a_2x^2 to approximate the exact displacement of the bar. Assume the area of the smaller bar to be equal to 10mm^2, L=1m, E=200GPa, and P=10kN. Compare with the exact displacement field.
  9. A beam of length L Young’s modulus E, and moment of inertia I is fixed at both ends and subjected to a uniformly distributed load q and a concentrated force P at midspan. Use the Rayligh Ritz method with a trial function u=A\sin\left(\frac{\Pi x}{L}\right) where A is a real number to approximate the midspan deflection of the beam.
  10. Consider a uniaxial bar of length L, Young’s modulus E, and variable cross sectional area A=A_0\left(2-\frac{x}{L}\right). The bar is fixed at one end and free at the other end. The free end is subjected to a constant force P and a uniformly distributed force q. Using the Rayleigh Ritz method, and assuming a quadratic trial displacement function, determine the approximate displacement of the bar.

Leave a Reply

Your email address will not be published.