MATLAB Answers

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

12 views (last 30 days)
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!

  0 Comments

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 27 Jun 2020
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.

  7 Comments

Show 4 older comments
César Solla González
César Solla González on 27 Jun 2020
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
John D'Errico on 27 Jun 2020
While I know the wiki page claims the two ways of formulating the problem are equivalent,
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:
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. :)

Sign in to comment.


Translated by