# Euler-Bernoulli cantilever beam with variable flexural rigidty

I am trying to find the deflection of cantilever beam with variable flexural rigidity (i.e. EI(x) instead of constant EI) starting from this equation: $$EI(x)\frac{\partial^2y}{\partial x^2} = M(x)$$, where $$M(x)$$ itself is found by integrating forces (distributed or concentrated) twice. I consider a rectangular distributed force in the form $$q(x)=\frac{w}{l} x-w$$ which can also be any other possible distributions. The variable flexural rigidity in my case is in the following form:

EI[x_]:= 0.26 + 9.12*10^6 (4.21875*10^-11 + 3.01401*10^-7 (2.68801 Sqrt[0.107 + x] - 2.06557 (0.107 + x) - 18.8981 (0.107 + x)^2 + 50.1011 (0.107 + x)^3 - 58.6458 (0.107 + x)^4)^2);


I write the main equation in Mathematica as follows:

qq=(w/L)*x - w;
eq1 = EI[x]*y''[x] == Integrate[Integrate[qq , x], x] + c1*x + c2;


where c1 and c2 are the expected integration constants. For constant flexural rigidity, I easily get an answer using this code, but it fails with no solution for nearly all variable flexural rigidties:

sol1 = First@DSolve[{eq1, y[x], x];
y1 = y[x] /. sol1;
y1 = y1 /. {C[1] -> c3, C[2] -> c4};(*Renaming the new constants*)
First@Solve[{(y1 /. x -> 0) == 0, (D[y1, x] /. x -> 0) == 0, (D[D[y1, x], x] /. x -> L) == 0, (D[D[D[y1, x], x], x] /. x -> L) == 0}, {c1, c2, c3, c4}]; (*Four boundary conditions are defined, i.e. no deflection and slope at x = 0 and no shear force and moment at x = L*)
y1 = y1 /. %


I am wondering if there exists another method to solve such problems effectively or rewriting my equations.

Edit

The provided solutions by Hugh and Rudy did actually yield a response, although the plot of deflection (either based on absolute value or real parts) becomes discontinuous around 0.2 - 0.4 when using Mathematica 11.0. Using Mathematica 12.2.0, the discontinuities are gone.

ClearAll[EI, y, x, qq, mm];
EI[x_] := 0.26 + 9.12*10^6 (4.21875*10^-11 + 3.01401*10^-7 (2.68801 Sqrt[0.107 + x] - 2.06557 (0.107 + x) - 18.8981 (0.107 + x)^2 + 50.1011 (0.107 + x)^3 - 58.6458 (0.107 + x)^4)^2);
Plot[EI[x], {x, 0, 1}, Frame -> True,
FrameLabel -> {"Position", "Bending Stiffness"}]
qq = (w/L)*x - w;
vv = Integrate[qq, x] + c1;
vv1 = vv /. First@Solve[vv == 0 /. x -> L, c1];
mm = Integrate[vv1, x] + c2;
mm1 = mm /. First@Solve[{mm == 0 /. x -> L}, {c2}];
sol = y[x] /.
First@DSolve[{EI[x] y''[x] == mm1, y[0] == 0, y'[0] == 0}, y[x], x];
Plot[Evaluate[Abs[sol //. {L -> 1, c1 -> (L w)/2, w -> 1}]], {x, 0,
0.8}, Frame -> True, FrameLabel -> {"Position", "Deflection"}]


## 2 Answers

Why have you integrated your shear force twice? I think the bending moment is equal to one integral of the shear force (if I remember my beam theory correctly).

Here is my calculation:

 ClearAll[EI, y, x, qq, mm];
EI[x_] :=
0.26 + 9.12*10^6 (4.21875*10^-11 +
3.01401*10^-7 (2.68801 Sqrt[0.107 + x] - 2.06557 (0.107 + x) -
18.8981 (0.107 + x)^2 + 50.1011 (0.107 + x)^3 -
58.6458 (0.107 + x)^4)^2);
Plot[EI[x], {x, 0, 1}, Frame -> True,
FrameLabel -> {"Position", "Bending Stiffness"}]


So the bending stiffness increases dramatically as we get to the end of the beam. I am assuming a 1 m long beam.

Now integrate the shear force to get the bending moment

 qq = (w/L)*x - w;
mm = Integrate[qq, x] + c1

(*   c1 - w x + (w x^2)/(2 L)  *)


Let's assume this is a cantilever and at x = L the bending moment is zero

Solve[mm == 0 /. x -> L, c1]

(*  {{c1 -> (L w)/2}}  *)


Now solve the differential equation assuming clamped at x = 0:

sol = y[x] /.
First@DSolve[{EI[x] y''[x] == mm, y[0] == 0, y'[0] == 0}, y[x], x];


I am not printing out the solution as it is very complicated. Let's plot it for some guessed values. There are some small imaginary parts (due to numerical error) which I remove.

Plot[Evaluate[Re[sol //. {L -> 1, c1 -> (L w)/2, w -> 1}]], {x, 0, 1},
Frame -> True, FrameLabel -> {"Position", "Deflection"}]


This looks about right. Does that help?

• Thanks Hugh, this indeed helped me to get the answer easily. Actually, integrating from distributed forces (say, pressure p(x)) generates shear forces, and integrating them yields moments. By doing so, I obtained the new moment and then I followed your solution. However, although the length of my beam is less than 0.2, when I follow your exact code, I obtain discontinuities in the last plot, whereas yours is continuous. Am I missing something? Commented Jan 6, 2021 at 18:57
• @Masoud Ah I understand. The qq is pressure not shear force. Integrating twice is then correct. I suggest you edit your question and add your current problem. Go to the end of the question type edit in bold and show us what your problem with discontinuities is.
– Hugh
Commented Jan 6, 2021 at 19:02
• I was actually referring to your final plot. I simply type your codes, but the plot I get has some discontinuities between 0.2 and 0.4, which is not reflected in your final plot. The same thing occurs with my new moments as I get some "Infinite expression 1/0 encountered". Commented Jan 6, 2021 at 19:28
• Don't know why that happens. The fact that there are discontinuities in the moments needs investigating first. You have integrated twice in your new version? I still can't do much unless I see your code.
– Hugh
Commented Jan 6, 2021 at 19:32
• I tried Mathematica 12.2 instead of my 11.0 version, and everything seems fine now. Thanks again. Commented Jan 6, 2021 at 20:58

I think I see three things going on.

1. A mismatched { in the definition of sol1
2. For algebraic solutions, I think it might work better to use = instead of := in the definition of EI[x_]
3. The substitution is being applied before the derivative in your second to last line. Keeping y1 as function and using prime notation cleans it up a bit.

Putting that all together I get:

ClearAll[y1]
ClearAll[EI]
EI[x_] = 0.26 + 9.12*10^6 (4.21875*10^-11 + 3.01401*10^-7 (2.68801 Sqrt[0.107 + x] - 2.06557 (0.107 + x) - 18.8981 (0.107 + x)^2 + 50.1011 (0.107 + x)^3 - 58.6458 (0.107 + x)^4)^2);

qq = (w/L)*x - w;
eq1 = EI[x]*y''[x] == Integrate[Integrate[qq, x], x] + c1*x + c2

sol1 = First@DSolve[eq1, y[x], x]
y1[x_] = Evaluate[y[x] /. sol1]
y1[x_] = y1[x] /. {C[1] -> c3, C[2] -> c4}
Solve[{y1[0] == 0, y1'[0] == 0, y1''[L] == 0, y1'''[L] == 0}, {c1, c2, c3, c4}]


I'm not sure the y1[x_] = y1[x] /. {C[1] -> c3, C[2] -> c4} line is necessary.

• Thanks Rudy. The mentioned tips were enough to easily get an answer in less than a minute. And the final plot by your method was the same as that generated using Hugh's approach. Commented Jan 6, 2021 at 19:11