10 views (last 30 days)

I try to solve the ODE shown on the figure. EI and GK are constants related to the beam cross section and properties, and 'm' is a value of an applied distirbuted external torque applied in the middle of a beam that has fixed-fixed conditions (i.e. at z=L/2). I want to obtain the distribution of Circulatory Torque and Warping Torque, as shown in the figure.

The equation to solve is:

Importantly, it has been derived from:

I believe I need to solve the first equation, with the boundary conditions, in order to plot Tc(z) and Tw(z), which is also given here:

The beam starts at z=0 and ends at z=252.

L=252 inches

EI = 1.1688e+10

GK=1.3907e+06 kip-inch^2

GK, EI and L have all been verified.

I haven't really been able to get very far, as I'm a bit of newbie when it comes to MATLAB. Any help would gratefully be appreciated.

- The Warping Moment is 0 at z=L/4
- The Warping Moment is 0 at z=3L/4
- Theta θ is 0 at z=0
- Theta θ is 0 at z=L

I have also included a figure related to the boundary conditions, if it could help...

syms t(z)

GK=1.3907e+06;

EI=1.1688e+10;

L=252;

m=1;

c1=t(0)==0;

c2=t(L)==0;

t2=diff(s,z,2);

c3=t2(L/4)*EI==0;

c4=t2(3*L/4)*EI==0;

t(z)=dsolve(diff(t,z,4)==m/(EI)+(GK/EI)*diff(t,z,2),c1,c2,c3,c4) %theta as a function of z

dsdz(z)=diff(s(z),z) %ds/dz

d2sdz2(z)=diff(dsdz(z),z) %d^2s/dz^2

d3sdz3(z)=diff(d2sdz2(z),z) %d^3s/dz^3

counter=1;

for z=0:1:252/2

T_circulatory(counter,1)=(GK)*dsdz(z)-EI*d3sdz3(z);

T_warp(counter,1)=-EI*d3sdz3(z);

counter=counter+1;

end

z=0:1:252/2;

plot(z,T_circulatory)

Bjorn Gustavsson
on 1 Nov 2019

Edited: Bjorn Gustavsson
on 1 Nov 2019

In general you convert a high-order ODE to a set of first order ODEs like in this example below (I will use a simple equation of motion ):

function dydt = ODE1orderex(t,y,CKG)

dydt(1) = y(2);

dydt(2) = -CKG(1)*y(2)^2 - CKG(2)*y(1) - CKG(3)

dydt = dydt(:);

end

You will have to do the same for your fourth-order ODE, the procedure is simple once you get to it.

As for solving your problem numerically you will have to look at the help and documentation for the boundary-value problem solvers bvp4c and bpv5c if you cannot get the symbolic solver to do this for you. Otherwise these beam-problems typically have "reasonably neat" solutions...

HTH

Bjorn Gustavsson
on 3 Nov 2019

What you can do is to calculate all the differences with matlabs gradient function and plug in those you need into the fourth-order ODE-equation. If that is reasonably satisfied you're good to go! (You might have to compare the deviation from equality with the amplitude of the different terms if they are very big, etc...)

HTH

Bjorn Gustavsson
on 22 Nov 2019

Sorry for not getting back to you. In addition to checking the gradients in the beam, you might

have to check the boundary conditions (all of them), and possibly also the values of your constants (I regularly forget what scaling I've decided to use for a particular experiment). Just to proceed with something. Then if you've made this much you surely have some-one in-house (at your department?) that could take a more competent look at this than I can...

HTH

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.