1
$\begingroup$

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"}]

Discontinuities of plot

enter image description here

$\endgroup$

2 Answers 2

2
$\begingroup$

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"}]

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"}]

Deflection

This looks about right. Does that help?

$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – Masoud
    Commented Jan 6, 2021 at 18:57
  • $\begingroup$ @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. $\endgroup$
    – Hugh
    Commented Jan 6, 2021 at 19:02
  • $\begingroup$ 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". $\endgroup$
    – Masoud
    Commented Jan 6, 2021 at 19:28
  • $\begingroup$ 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. $\endgroup$
    – Hugh
    Commented Jan 6, 2021 at 19:32
  • $\begingroup$ I tried Mathematica 12.2 instead of my 11.0 version, and everything seems fine now. Thanks again. $\endgroup$
    – Masoud
    Commented Jan 6, 2021 at 20:58
1
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$
    – Masoud
    Commented Jan 6, 2021 at 19:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.