Problems Implementing ddesd with a system of equations.

6 次查看（过去 30 天）

Megan Frisbey 2023-6-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1988933-problems-implementing-ddesd-with-a-system-of-equations

编辑： Torsten 2023-6-29

I am attempting to model a system of equations to capture a single rigid body dynamics problem. I want to use functions for the forces and moments acting on the body thus creating a state dependent function. First question - this does make it a DDE problem correct? Next I have implemented a simple rudemetary model using the examples given online but I am missing/not understanding something when it comes to the history function. I took this as a way to define the initial conditions, but after seeing the results I believe this is not correct or I have not defined it properly? To back out what was going on I dusted off my ODE solver using ode45 with the same equations and such and determined that the ddesd solver is using the history function and forcing those values for my delay values -- Why?!? It doesn't matter how I change my delay functions it has no effect on the solution and appears to be only using the value for the first position of the vector I define in the histroy function. Obvously I'm missing something, but what?

I appreciate the help!

below is the dde system I've tried to implement.

t0 = 0                                                                                                                                                                                
tf = 10
t = linspace(t0,tf,200)
  
sol = ddesd(@ddefun,@delay,@history,t)
  
t = sol.x
y = sol.y
%I've ommited all the plotting but can add this back in if needed
function d = delay(t,y)
   d = [delayForce(t,y) delayMoment(t,y)]
end
 
function dF = delayForce(t,y)
   g = 32.174        %ft/s^2
   mB = 50;          %lbs
   rho = 0.002377    %SL english
 % FLaunch = getLaunchForce(t);
% FAero = getAeroForce(y,rho)
% FWeight = mB  %mB*g
% FThrust = -FAero(1)*20
% Fbx = FLaunch + FAero(1) + FThrust;
% Fby = FAero(2);
% Fbz = FAero(3) + -FWeight;
  
   Fbx = 0;  %getLaunchForce(t)
   Fby = 0;
   Fbz = 0;
  
dF = [Fbx Fby Fbz];
  
end
  
function dM = delayMoment(t,y)
 
   Lbx = 0
   Mbx = 0
   Nbx = 0
 
   dM = [Lbx Mbx Nbx]
 
end 
function v = history(t)
% U = y(1); V = y(2); W = y(3); x = y(7); y = y(8); z = y(9);
% P = y(4); Q = y(5); R = y(6); phi = y(10); theta = y(11); psi = y(12);
   g = 32.174        %ft/s^2
   v = transpose([10 0 0 0 0 0 0 0 0 0 0 0]);
end
function dydt = ddefun(t,y,Z)
   dydt = zeros(12,1);
 
   mB = 50/32.174;          %lbs
   mPlate = 6;
   rPlate = 6;
   hPlate = 1;
  
   IBxx = 1/12*mB*(3*rPlate*rPlate + hPlate*hPlate);
   IByy = 1/12*mPlate*(3*rPlate*rPlate + hPlate*hPlate);
   IBzz = 1/2*mPlate*rPlate*rPlate;
   
 % u1 = y(1); u2 = y(2); u3 = y(3); q1 = y(7); q2 = y(8); q3 = y(9);
 % u4 = y(4); u5 = y(5); u6 = y(6); q4 = y(10); q5 = y(11); q6 = y(12); 
 
   dydt(1) = Z(1,1)*(1/mB) + y(2)*y(6) - y(3)*y(5);  
   dydt(2) = Z(1,2)*(1/mB) + y(3)*y(4) - y(1)*y(6);
   dydt(3) = Z(1,3)*(1/mB) + y(1)*y(5) - y(2)*y(4);
   dydt(4) = (Z(1,4) + (IByy - IBzz)*y(5)*y(6))*(1/IBxx);  
   dydt(5) = (Z(1,5) + (IBzz - IBxx)*y(4)*y(6))*(1/IByy);
   dydt(6) = (Z(1,6) + (IBxx - IByy)*y(4)*y(5))*(1/IBzz);
 
   dydt(7) =  y(1) + y(6)*y(8) - y(5)*y(9);
   dydt(8) =  y(2) + y(4)*y(9) - y(6)*y(7);
   dydt(9) =  y(3) + y(5)*y(7) - y(4)*y(8);
   dydt(10) = y(4) + y(5)*sin(y(10))*tan(y(11)) +y(6)*cos(y(10))*tan(y(11));
   dydt(11) = y(5)*cos(y(10)) -y(6)*sin(y(10));
   dydt(12) = y(6)*(cos(y(10))/cos(y(11)))+y(5)*(sin(y(10))/cos(y(11)));
  
 
end

Below here is the code where I use ode45 and having tuned the force functions so that it matches the incorrect output from ddesd.

function [T,Y] = call_osc_1BodyNE()                                                                                                                                                   
tspan = [0 10];
initVec = [10 0 0 0 0 0 0 0 0 0 0 0]
y1_0 = 2;
y2_0 = 0;
  
[T,Y] = ode45(@osc_1BodyNE_alt,tspan,initVec);
%I've ommited all the plotting but can add this back in if needed

function dydt = osc_1BodyNE_alt(t,y)                                                                                                                                                  
 
dydt = zeros(12,1);
 
g = 32.174;
weight = 50;
rho = 0.002377;
massB = weight/g;
mB = massB    %Generic
mPlate = 6;
rPlate = 6;
hPlate = 1;
%Calculate Forces
% FLaunch = getLaunchForce(t);
% FAero = getAeroForce(y,rho)
% FWeight = massB*g
% FThrust = FAero(1)
% 
% F_Applied(1) = FLaunch + FAero(1)+FThrust;
% F_Applied(2) = FAero(2);
% F_Applied(3) = FAero(3) + -FWeight;
  
% F_Applied(1) = getLaunchForce(t);
F_Applied(1) = 10;
F_Applied(2) = 10;
F_Applied(3) = 10;
  
% + FAero(3)*1.2;  %
% Lift = FAero(3)
 
%Mass Moments
Ib = zeros(3,3);
Ib(1,1) = 1/12*massB*(3*rPlate*rPlate + hPlate*hPlate);
Ib(2,2) = 1/12*mPlate*(3*rPlate*rPlate + hPlate*hPlate);
Ib(3,3) = 1/2*mPlate*rPlate*rPlate;
IBxx = Ib(1,1);
IByy = Ib(2,2);
IBzz = Ib(3,3);
 
%Calculate Moments
 
M_Applied(1) = 10;
M_Applied(2) = 10;
M_Applied(3) = 10;
  
%body
% U = y(1); V = y(2); W = y(3); x = y(7); y = y(8); z = y(9);
% P = y(4); Q = y(5); R = y(6); phi = y(10); theta = y(11); psi = y(12);
  
dydt(1) = F_Applied(1)*(1/mB) + y(2)*y(6) - y(3)*y(5);
dydt(2) = F_Applied(2)*(1/mB) + y(3)*y(4) - y(1)*y(6);
dydt(3) = F_Applied(3)*(1/mB) + y(1)*y(5) - y(2)*y(4);
 
dydt(4) = (M_Applied(1) + (IByy - IBzz)*y(5)*y(6))*(1/IBxx);
dydt(5) = (M_Applied(2) + (IBzz - IBxx)*y(4)*y(6))*(1/IByy);
dydt(6) = (M_Applied(3) + (IBxx - IByy)*y(4)*y(5))*(1/IBzz);
 
dydt(7) =  y(1) + y(6)*y(8) - y(5)*y(9);
dydt(8) =  y(2) + y(4)*y(9) - y(6)*y(7);
dydt(9) =  y(3) + y(5)*y(7) - y(4)*y(8);
  
dydt(10) = y(4) + y(5)*sin(y(10))*tan(y(11)) +y(6)*cos(y(10))*tan(y(11));
dydt(11) = y(5)*cos(y(10)) -y(6)*sin(y(10));
dydt(12) = y(6)*(cos(y(10))/cos(y(11)))+y(5)*(sin(y(10))/cos(y(11)));
end

27 个评论
显示 25更早的评论隐藏 25更早的评论

Megan Frisbey 2023-6-29

Obtaining accurate solutions for fluid-structure interaction problems is very complex. The equations that rule CFD is the Navier Stokes equations - PDE's, and also for structures more PDE's. I am attempting a much more simplified less refined method of doing this by precomputing the values of coefficient of drag and lift on a body at different flow field states, and by only considering the rigid body dyanimcs (no bending etc). For a very crude model I can reduce this down to some peice wise function. For a better model it will be table look up based off of conditions. I could attempt to map the solutions surface and formulate into some form of an equation but this I believe is more trouble than it is worth. If the time step is small enough for this problem using the previous state values to determine the current force and moment values would be "good enough"

As to the time step not being quantified I'm not sure what you mean. Plenty of numerical schemes don't quanitfy a time step outside of stability bounds and accuracy of solution...What are you meaning by this?

Megan Frisbey 2023-6-29

In my question I specifically referred to using ode45 which is only dependent on the previous time step

"ode45 is based on an explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair. It is a single-step solver – in computing y(tn), it needs only the solution at the immediately preceding time point, y(tn-1) [1], [2]."

Torsten 2023-6-29

编辑：Torsten 2023-6-29

That's true. But it will call your ode function not only for t(n-1) and y(t(n-1)) to compute y(t(n)), but also for some "intermediate" points. You can see this here for the classical Runge-Kutta method - it's similarily done in ode45:

https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Sulaymon Eshkabilov 2023-6-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1988933-problems-implementing-ddesd-with-a-system-of-equations#answer_1264348

在 MATLAB Online 中打开

Why you are not using the time span with a predefined time step, in call_osc_1BodyNE(), that would give you consistent time steps in your found solutions from ode45 e.g.:

function [T,Y] = call_osc_1BodyNE()                                                                                                                                                   
t0 = 0;                                                                                                                                                                              
tf = 10;
tspan = linspace(t0,tf,200);
initVec = [10 0 0 0 0 0 0 0 0 0 0 0];
y1_0 = 2;
y2_0 = 0;
  
[T,Y] = ode45(@osc_1BodyNE_alt,tspan,initVec);
%I've ommited all the plotting but can add this back in if needed

2 个评论
显示无隐藏无

Megan Frisbey 2023-6-29

No good reason I suppose other than it isn't really necessary? Is there a compelling reason to do this using ode45? The solution doesn't change and ode45 seems to be doing fine managing the time step...

Megan Frisbey 2023-6-29