Using ODE45 to Solve a System of Equations Outputting Linear non-Ideal Results

조회 수: 2 (최근 30일)
John
John 2024년 2월 17일
편집: Torsten 2024년 2월 20일
I am currently modeling a thruster in 1D. And using a set of governing equations found in a paper by Tommaso Misuri, the equations are 3.1 to 3.6. The form I have rearranged the equations in to is:
M(Y) * dY/dz = F(Y)
Where M(Y) is the mass matrix, dY/dz are the collection of diff equations, and F(Y) is the vector of RHS.
The following is the function for solving this system.
function dY_dz = odefcn(z, Y, sigma, rho, R, alpha, muo, j)
v = Y(1);
Tk = Y(2);
B_self = Y(3);
c_v = (3/2)*R; % Specific heat at specific volume (ideal monatomic)
c_p = c_v + R; % Specific heat ratio at specific pressure (ideal monatomic)
% F(Y) Right-Hand Side
rhs_eq = zeros(size(Y));
rhs_eq(1) = (1/rho)*(j*B_self); % (1) Momentum Equation
rhs_eq(2) = (j*v*B_self) + (((j^2)/sigma)*(1 - alpha)); % (2) Energy Equation
rhs_eq(3) = -muo*j; % (3) Maxwell's equation for self induced mag field
% Mass Matrix
M = zeros(length(Y), length(Y));
M(1, 1) = v; % (1) Momentum Equation
M(2, 1) = rho*(v^2); M(2, 2) = rho*v*c_p; % (2) Energy Equation
M(3, 3) = 1; % (3) B_self
dY_dz = M \ rhs_eq;
end
And this is the setup for ODE45:
v_o = velocity;
T_o = Tk;
B_o = B_self;
Y0 = [v_o; T_o; B_o]; % Initial conditions
[z, Y] = ode45(@(z, Y) odefcn(z, Y, sigma, rho, R, alpha, muo, j), ...
z, ... % tspan
Y0); % Initial conditions
And these are the results:
As you can see these are less than ideal, for an electric propulsion thruster velocity increase should be dramatically higher than this (and non-linear). Ideal results should look similar to the paper linked previously.
My initial thoughts as to why they are incorrect is because the function I am using 'reuses' the v, T, and B_self variables, so basically the function does compute a new value but then as soon as it is looped back into the function it just reassigns it back to the initial value. This would explain why the results do not change greatly.
Would this be the reason why my code doesn't work? Or is it a different problem? My knowledge on ODE45 is limited. Any help will be appreciated.
  댓글 수: 10
John
John 2024년 2월 20일
As the paper I'm referencing from states that the continuity equation is: rho*v*A = constant. I just assumed constant is zero. Is this incorrect?
Regarding ode15s and ode23t neither of them worked unfortunately.
Torsten
Torsten 2024년 2월 20일
편집: Torsten 2024년 2월 20일
I just assumed constant is zero. Is this incorrect?
rho*v has unit kg/(m^2*s) and is the mass flow per unit area into the domain. Setting it to 0 is a very uninteresting case because your initial value either for rho or v at z=0 had to be 0.
rho(z)*v(z)*A = rho0*v0*A
where rho0, v0 are density and velocity at z = 0.
Thus your first equation must be
rho(z)*v(z) - rho0*v0 = 0
assuming A does not change with z.

댓글을 달려면 로그인하십시오.

답변 (1개)

Steven Lord
Steven Lord 2024년 2월 17일
The form I have rearranged the equations in to is:
M(Y) + dY/dz = F(Y)
Where M(Y) is the mass matrix, dY/dz are the collection of diff equations, and F(Y) is the vector of RHS.
Can you confirm this is the correct form of the equations? When a mass matrix is involved usually it multiplies the derivative, M(Y)*dy/dz = F(Y), rather than adding to it as you've written it. If you are solving the case with a mass matrix multiplying the derivative, rather than trying to handle it yourself in your ODE function I'd use odeset to set the Mass option and pass that options structure into ode45.
If the addition form is correct, your code is incorrect. The right-hand side should be F(Y)-M(Y).
This documentation page works through an example with a mass matrix; since you say your knowledge of ode45 is limited perhaps you could use this as a guide or model for how to solve an ODE with a mass matrix.
  댓글 수: 2
John
John 2024년 2월 19일
Hello Steven,
Thankyou for pointing this out, that was my error, it should've been multiply. I have changed this.
Torsten
Torsten 2024년 2월 19일
In your code, you multiplied from the very beginning ...

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

태그

제품


릴리스

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by