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

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.

[Torsten]

Nothing moves without a driving force - in this case a pressure gradient (that you did not implement).

[Sam Chak]

Hi @John

While I don't consider myself an expert in magnetohydrodynamics, I would like to assist you by suggesting that you insert the correct values for the parameters (Sets #1 and #2) in your code. By doing so, we can verify if any mistakes are present and ensure the accuracy of the results in Figure 4.

%% Parameters (set #1)
velocity = 1;
Tk       = 1;
B_self   = 1;

v_o     = velocity;
T_o     = Tk;
B_o     = B_self;
zspan   = [0, 0.12];
Y0      = [v_o; T_o; B_o];  % Initial conditions
[z, Y]  = ode45(@(z, Y) odefcn(z, Y), zspan, Y0);
plot(z, Y), grid on

function dY_dz = odefcn(z, Y)
    %% Parameters (set #2)
    sigma    = 1;
    rho      = 1;
    R        = 1;
    alpha    = 1;
    muo      = 1;
    j        = 1;

    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

[John]

Hello,

Good idea Sam, I've been working through this for the past few hours - I haven't got the results I wanted but some progress has been made. I have a feeling there are some missing elements from the paper I linked. As Torsten suggested I didn't include the pressure gradient term. I didn't include this dp/dz term as the equation of state was not differentiated, and I set dp/dz in momentum equation to 1. So I will follow Torsten's advice and differentiate the equation of state eq and see what effect that has, notable on the temperature and velocity.

Outside of this, I am getting some results for velocity and tempature terms, although the magnetic field term B is stagnating. This can be seen below. I can assume the B-field is a straight line because its just doing 1 \ -(muo*j).

Below is the slighly modified function. I should note that the figure above was obtained when E = 50 (I think). The E-field value is what defines the model from what I've found. If E is set to one then the entire plot is equal to 1 across z.

function dY_dz = odefcn(z, Y)
    v         = Y(1);
    Tk        = Y(2);
    B_self    = Y(3);

    % Parameters (set #2)
    A = 0.07*0.02; % [m]
    E = 1; % [Volts/m]
    sigma    = ((1.5085e-2)*Tk^(3/2))/(10); % ((1.5085e-2)*Tk^(3/2))/(10)
    rho      = v*A;
    R        = 8.314;
    alpha    = 0.5; % between 0.2 to 0.5
    muo      = 1.25663706e-6;
    j        = sigma*(E-(v*B_self)); % sigma*(E-(v*B_self))

    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 -((1/rho)*(rho*R*Tk))+
    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

Will try to play around with the equations and see if my results change at all.

[Torsten]

As Torsten suggested I didn't include the pressure gradient term. I didn't include this dp/dz term as the equation of state was not differentiated, and I set dp/dz in momentum equation to 1. So I will follow Torsten's advice and differentiate the equation of state eq and see what effect that has, notable on the temperature and velocity.

You can leave the equation of state in its algebraic form. Most of the integrators accept differential-algebraic systems of equations, thus a mixture of differential and algebraic equations.

[John]

As there is a dp/dz in the momentum equation, wouldn't this make the mass matrix 4x4, and then the RHS (F(Y)) 3x1 vector? May seem like a stupid question, but honestly I lack experience in modelling so this is difficult to get my head round.

[Torsten]

The row of the mass matrix for the equation of state is a 1x4 vector of zeros. The right-hand side component is p-rho*R*T. Multiplying out gives 0 = p-rho*R*T.

[John]

Right I'm following. The following code is my updated function, including each one of the governing equations provided in the linked reference. Is this the correct approach?

function dY_dz = odefcn(z, Y)
    rho     = Y(1); %     rho         = v*A;
    v       = Y(2);
    T       = Y(3);
    p       = Y(4);
    j       = Y(5);
    B_self  = Y(6);

    % Parameters (set #2)
    A           = 0.07*0.02; % [m]
    E           = 250; % [Volts/m]
    R           = 8.314;
    alpha       = 0.2; % between 0.2 to 0.5
    muo         = 1.25663706e-6;
    sigma       = ((1.5085e-2)*T^(3/2))/(10);

    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) = rho*v*A;                                        % (1) Continuity Equation
    rhs_eq(2) = (1/rho)*(j*B_self);                             % (2) Momentum Equation
    rhs_eq(3) = (j*v*B_self)+(((j^2)/sigma)*(1-alpha));         % (3) Energy Equation
    rhs_eq(4) = (rho*R*T)-p;                                    % (4) Equation of State
    rhs_eq(5) = (sigma*(E-(v*B_self)))-j;                       % (5) Ohm's law
    rhs_eq(6) = -muo*j;                                         % (6) Maxwell equation
    
    % Mass Matrix
    M         = zeros(length(Y), length(Y));
                                                            % (1) Continuity Equation
    M(2, 2)   = v;              M(2, 4) = 1/rho;            % (2) Momentum Equation
    M(3, 2)   = rho*(v^2);      M(3, 3)   = rho*v*c_p;      % (3) Energy Equation
                                                            % (4) Equation of State
                                                            % (5) Ohm's law
    M(6, 6)   = 1;                                          % (6) Maxwell equation


    dY_dz     = M\rhs_eq;
end

The problem now is that I'm not getting any results at all. The Y matrix comes back as NaN.

[Torsten]

You want to have rho*v*A = 0 ? That's what the first row of your DAE system says.

If the code does not work after the necessary change for the first equation, please include your call to "ode15s" or "ode23t" so that your code becomes executable. Note that "ode45" cannot be used since your matrix M is singular.

[John]

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]

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.

Answer 1

[Steven Lord]

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.

[John]

Hello Steven,

Thankyou for pointing this out, that was my error, it should've been multiply. I have changed this.

[Torsten]

In your code, you multiplied from the very beginning ...