MATLAB Answers

Solving 4th order ode problem

10 views (last 30 days)
Henry Chinaski
Henry Chinaski on 1 Nov 2019
Commented: Bjorn Gustavsson on 22 Nov 2019
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.
  1. The Warping Moment is 0 at z=L/4
  2. The Warping Moment is 0 at z=3L/4
  3. Theta θ is 0 at z=0
  4. 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)

  0 Comments

Sign in to comment.

Answers (1)

Bjorn Gustavsson
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

  4 Comments

Show 1 older comment
Bjorn Gustavsson
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
Henry Chinaski
Henry Chinaski on 11 Nov 2019
Thanks for your suggestion Bjorn. I've done exactly that to check the "solution" - the gradient function is indicating that the different rows of theta are equal to the different derivatives (see all of the plots).
Unfortunately, I haven't been able to get a solution that is similar to that on the worksheet :-( I'm guessing the initial assumptions and boundary conditions are responsible, but I can't see what I have done wrong.
I spent some time putting together the moment function ("m") using the heaviside function, which is based on the distribution given above.. but still no hope! I would be very grateful for any more suggestions...
function [theta,z,T,aaa,ccc] = mat4bvp3
close all
clear all
L=252;
GK=1.3907e+06; %0.4*Ec*K
EI=1.1688e+10; %Youngs Modulus multiplied by Sectorial
options = bvpset('RelTol', 1e-3, 'AbsTol', 1e-3, 'NMax', 250);
solinit = bvpinit(linspace(0,L,20),[0 0 0 0]);
sol = bvp4c(@mat4ode,@mat4bc,solinit,options);
z = linspace(0,L);
theta = deval(sol,z);
figure
plot(z,theta(1,:)); hold on
figure
plot(z,theta(2,:)); hold on
aaa=gradient(theta(1,:),z);
plot(z,aaa(1,:)); hold on
figure
plot(z,theta(3,:)); hold on
bbb=gradient(aaa,z);
plot(z,bbb(1,:)); hold on
figure
plot(z,theta(4,:)); hold on
ccc=gradient(bbb,z);
plot(z,ccc(1,:)); hold on
T=(GK)*theta(2,:)-(EI)*theta(4,:);
figure
plot(z,T);
xlim([0 L/2])
function dxdz = mat4ode(~,z)
GK=1.3907e+06; %0.4*Ec*K
EI=1.1688e+10; %Youngs Modulus multiplied by Sectorial
T=1; %assume a torque of 1 kip-inch has been applied to the beam
L=252; %length of the beam, inches
%m=((T/2)*z(1)-(T*L/8));
m=((heaviside((L/2)-z(1))*2-1)*((T/2)*z(1)))-(T*L/8)+(heaviside(z(1)-(L/2)))*((T*L)/2);
dxdz=[z(2) %dtheta/dz
z(3) %d^2theta/dz^2
z(4)
m/EI+(GK/EI)*z(3)]; %4th order equation
function res = mat4bc(za,zb)
res=[za(1) % theta=0 at z=0
za(2) % dtheta/dz=0 at z=0
zb(1) % theta=0 at z=L
zb(2)]; % dtheta/dz=0 at z=L
Bjorn Gustavsson
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.