Solving 4th order ode problem

4 vues (au cours des 30 derniers jours)
Henry Chinaski
Henry Chinaski le 1 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)

Réponses (1)

Bjorn Gustavsson
Bjorn Gustavsson le 1 Nov 2019
Modifié(e) : Bjorn Gustavsson le 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 commentaires
Henry Chinaski
Henry Chinaski le 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 le 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

Connectez-vous pour commenter.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by