1D Linear Elasticity Equivalent Stiffness Solution not being solved (6 equations and 6 vars)

Hey guys, I am unable to make my program get a result of w as a function of EI and L

syms w x EI A B C D L

x0=0;

eqns = [EI*diff(w,3)==heaviside(x-x0)+A, EI*diff(w,2)==heaviside(x-x0)*(x-x0)+A*x+B, EI*diff(w)==heaviside(x-x0)*(x-x0)^2/2+A*x^2/2+B*x+C, EI*w==heaviside(x-x0)*(x-x0)^3/6+A*x^3/6+B*x^2/2+C*x+D, w(x==L/2)==0, diff(w(x==-L/2))==0];
ctes = [A B C D];

[A, B, C, D]=solve(eqns, ctes);

w(x==0)

This is what I've coded. As you can see it is 4 equations and 2 boundary conditions (The 2 last ones) I want it to solve the 4 constants A, B, C and D, so I get the full result of the function w(x) with those Constants solved and then finally evaluate w for x=0 (in this case). But it is giving me this result:

ans =
 
[ empty sym ]

I am pretty sure it is something about syntax in the equations definition, but I can't find it. Could you give it a look guys? I would really appreciate. This example should give the following solution w(x=0) = L^3/(192*EI)

Thank you very much!

Answer 1

[John D'Errico]

These are differential equations, a boundary value problem.

We use dsolve to solve it. For example, a simple example is:

syms x(t)
Eq = diff(x,t,2) == 1 + t;

dsolve(Eq,x(0) == 0,x(1) == 2)
ans =
(t*(t^2 + 3*t + 8))/6

Now, what happens when you just throw it at solve?

solve(Eq,x(0) == 0,x(1) == 2)
Warning: Unable to find explicit solution. For options, see help. 
> In sym/solve (line 317) 
ans =
Empty sym: 0-by-1

Finally, even though you used solve to try to find a solution, that last line, w(x==0) does not do what you think it does. In fact, in your code, even though you tried to return A,B,C,D, those values are not magically inserted into the relationship for w.

Oh, this is not how you create a boundary condition:

w(x==L/2)==0

Again, that does NOT set x to a value of L/2, and then substitute it into w. Instead, you would formulate that boundary condition for dsolve as

w(L/2) == 0

So, first, solve the ODE as a differential equation. Use dsolve.

[El Vapor]

Hey friend. Thank you very much, I didn't know about that command! I am trying to apply what you've told me. However, I am now getting a message saying that matlab is unable to find a solution :/ Do you know what could be happening? This is the brand new code:

syms w(x) x EI A B C D L

x0=L;

eqn = EI*diff(w,x,4)== dirac(x-x0);
Dw = diff(w,x);
DDw = diff(w,x,2);
DDDw = diff(w,x,3);
cond = [w(0)==0, Dw(0)==0, DDDw(L)==0, DDw(L)==0];

S = dsolve(eqn, cond);

w(0)

And the message I am receiving when trying to solve:

>> Elastica
Warning: Unable to find symbolic solution. 
> In dsolve (line 216)
  In Elastica (line 11) 
 
ans =
 
w(0)

[John D'Errico]

Ok. This is at least making slightly more sense. Start with the beam equation.

What is the domain? Is it the interval [0,L]? That appear to be the case.

We now see 4 boundary conditions, necessary for a 4th order ODE.

w(0) = 0 % pinned at 0, so no displacement at x==0.

w'(0) = 0 % so no slope at x==0.

At the right end point, thus X==L, the second and third derivatives are set to zero. So a "natural", or "free" boundary of sorts?

It appears as if you are trying to create a cantilevered beam, pinned at x ==0, but free at the right end?

Lets see, if we consider a picture, we can find an example of exactly such a cantilevered beam

https://en.wikipedia.org/wiki/Euler–Bernoulli_beam_theory

in the section on boundary conditions.

What loads are you wanting to place on the beam? At x == L, you want to place an infinite point load? And then you want to solve for the deflection of the beam at x==0? Since you have already established w(0) to be zero, if you could solve this problem, you would see w(0) must be zero.

So still a highly confusing question. Please explain the problem in terms of the engineering. What beam do you want to model? What are the loads? What do you know at the ends?

[El Vapor]

Sorry, sure I haven't explained much about the problem itself:

What I am trying to do is modelizing the equivalent stiffness of a beam into a spring. That stiffness depends on the different boundary conditions that the beam is supporting. As well as it depends on the distance between those conditions and the point of interest where I want to modelize the spring. I am taking E*I constant with x. I've made you a drawing so you can understand it better.

Finally, thank you for the time you are putting in helping me. I am constantly having to analize this kind of systems (vibrations) and would be very useful to automatize this proccess since I need to do many more different calculations later

PD: The final "w(0)", yeah I know that's equal to 0, it was just a test to check if the code could calculate it right. My very real concern is the calculation of w(x0), being x0 the distance to the point of interest where I want to modelize the spring

[John D'Errico]

Ok. Getting there. You want to treat a cantilevered beam (with initially constant EI over the length of the beam) as a spring? So a unit shear load at X == L? The deflection of the beam as a function of x will be

syms w(x) W EI L F
ODE = EI*diff(w,x,4) == 0;
dw = diff(w,1);
ddw = diff(w,2);
dddw = diff(w,3);

conds = [w(0) == 0, dw(0) == 0, ddw(L) == 0, -EI*dddw(L) == F];

wsol = dsolve(ODE,conds)
wsol =
-(F*x^3 - 3*F*L*x^2)/(6*EI)

So as a function of the load at L, we see a deflection of:

subs(wsol,x,L)
ans =
(F*L^3)/(3*EI)

Thus proportional to the shear load F.

As for the verificatinon at x == 0, we see:

subs(wsol,x,0)
ans =
0

So it is clamped. We could compute the first derivative too.

subs(diff(wsol),x,0)
ans =
0

And now you want to think of that as a spring. Are we getting closer?

[El Vapor]

Man, thank you very much, for real. Thanks to your example I've managed to do it with some different examples:

syms w(x) x EI L keq

x0=0;
assume(L,'positive')
EDO = EI*diff(w,x,4) == dirac(x-x0);
Dw = diff(w,x);
DDw = diff(w,x,2);
DDDw = diff(w,x,3);
CC = [w(-L/2) == 0, Dw(-L/2) == 0, w(L/2) == 0, Dw(L/2) == 0];

W = dsolve(EDO, CC);
keq = 1/subs(W,x,x0)

>> Elastica
 
keq =
 
(192*EI)/L^3

However, for some reason, I still cannot make it work with second and third derivatives of w(x) when using dirac's function. (I am using F=1, that's why F doesn't appear in my equation), I am also following the indications in the page you sent me:

I am using dirac's since it is much easier for me having all the boundary conditions == 0, Here is the code and the error I am getting:

syms w(x) x EI L keq

x0=L;
assume(L,'positive')
EDO = EI*diff(w,x,4) == dirac(x-x0);
Dw = diff(w,x);
DDw = diff(w,x,2);
DDDw = diff(w,x,3);
CC = [w(0) == 0, Dw(0) == 0, DDw(L) == 0, DDDw(L) == 0];

W = dsolve(EDO, CC)

>> Elastica
Warning: Unable to find symbolic solution. 
> In dsolve (line 216)
  In Elastica (line 11) 
 
W =
 
[ empty sym ]

[John D'Errico]

While I know the wiki page claims the two ways of formulating the problem are equivalent,

https://en.wikipedia.org/wiki/Euler–Bernoulli_beam_theory

it seems possible that dsolve is not happy with using a dirac delta at an end point to formulate the problem. For example, reading the answer to this question:

https://www.mathworks.com/matlabcentral/answers/18350-impulse-response-using-dirac-t

What you tried

syms w(x) x EI L keq
x0=L;
assume(L,'positive')
EDO = EI*diff(w,x,4) == dirac(x-x0);
Dw = diff(w,x);
DDw = diff(w,x,2);
DDDw = diff(w,x,3);
CC = [w(0) == 0, Dw(0) == 0, DDw(L) == 0, DDDw(L) == 0];

must fail though. As even the euler-bernouli page suggests it will. See this statement on the wiki page:

"Note that shear force boundary condition (third derivative) is removed, otherwise there would be a contradiction."

I'll admit that I tried a few things with no success to get dsolve to happily use dirac to implement this as an impulse response at the right hand end point.

So here, with 3 boundary conditions. Strangly, I could not get dsolve to solve this over the interval [0,L]. So here I solved it over a unit interval, and then transformed the problem. (I was also lazy and implicitly left EI == 1. That should not be very significant, except to make things a bit simpler to read.)

W1 = dsolve(diff(w,4) == dirac(x - 1),w(0) == 0,dw(0) == 0,ddw(1) == 0)
W1 =
x*(sign(x - 1)/4 + 1/4) - sign(x - 1)/12 - x^2*(sign(x - 1)/4 - C1 + 1/4) + x^3*(sign(x - 1)/12 - C1/3 + 1/12) - 1/12

As you should see, there is a constant of integration remaining in there: C1.

Now, transform to the interval [0,L].

syms u
Wu = subs(W1,x,u/L)
Wu =
(u*(sign(u/L - 1)/4 + 1/4))/L - sign(u/L - 1)/12 - (u^2*(sign(u/L - 1)/4 - C1 + 1/4))/L^2 + (u^3*(sign(u/L - 1)/12 - C1/3 + 1/12))/L^3 - 1/12

The sign function in there can be viewed as a heaviside function.

Does it satisfy the boundary conditions?

>> subs(Wu,u,0)
ans =
0

>> subs(diff(Wu,1),u,0)
ans =
0

>> subs(diff(Wu,2),u,L)
ans =
NaN

>> limit(diff(Wu,2),u,L,'left')
ans =
0

As you can see, on the right end I needed to take the limit as u --> L from the left.

If I now take the limit of the third derivative, as u approaches L from the left side...

>> limit(diff(Wu,3),u,L,'left')
ans =
-(2*C1)/L^3

I get something that is consistent with the use of an impulse response on that left end.

But honestly, I felt I had to work much harder to get to this point than I did by just imposing a shear load directly in the boundary conditions. I'm not sure it is worth it. :)

[El Vapor]

You have helped me a lot, I really appreciate it man. Thank you very much!

Answer 2

[Pablo García Fogeda Núñez]

Esto no te va a ayudar a aprobar Vibraciones.