Simulation of 1D Super-sonic flow through a nozzle using MacCormack method

1D SUPER SONIC NOZZLE FLOW SIMULATION: USING MACCORMACK METHOD

The objective of this project is to simulate a 1D quasi flow through a convergent-divergent nozzle, using the MacCormack method. For this, we need to define and derive:

1. The governing equations for 1D quasi flow through a nozzle

2. MacCormack method: A summary

3. Coding and analysis of results

For a better perspective, we are going to refer heavily from the John D.Anderson book, \"Computational Fluid Dynamics: The basics with applications\" for this problem.

So, let\'s start with deriving the necessary equations.

DERIVATION OF EQUATIONS

Before we proceed we have to first know the two types of equation forms that represent all physical aspects in a governing equation of a particular state or system.

• Non-conservative form
• Conservative form

Basically, there is not much of a difference between these two types of equations in the physical aspect, i.e, taking into account the physics of the solution. Both equations are the same and when a few steps of derivation are done, the conservative form of an equation can be turned into a non-conservative form and vice-versa. But numerically speaking, there are some differences and each method is used in one or the other situation with regards to the applications of CFD in a physical system.

The non-conservative form of an equation, with respect to computational techniques, handles variables of the flow parameters as direct effect in the formulation and discretization of the equation computationally.

Whereas, the conservative form of an equation, with respect to computational techniques, handles variables of the flow parameters as indirect effect in the formulation and discretization of the equation computationally, i.e, they are formulated more generally into something called \'solution vectors\' and we have to discretize terms using them and then derive the variables or \'primitive\' variables from those solution vectors.

This kind of methodology makes conservative equations more difficult to code or simulate computationally, but as we will learn along, it does have it\'s own advantages with respect to CFD applications. Simply speaking, conservative forms have dependent variables to be derived after, to get primitive flux variables and represents a more \'conservative\' approach in solving the problem, thereby conserving most of the variables\' physical properties throughout the problem. Non-conservative forms can be directly discretized in it\'s own variable format and are mostly conserved along the system time life, but are not preferred in complex problems where, for example, we need to account for shocks or chemical reactions, in which case conservative form is used.

Non-conservative form of equations:

We also know that, all physical problems can be represented within three equations or physical statements such as,

1. The continuity equation, where mass is conserved.
2. The momentum equation, where newton\'s second law is used.
3. The energy equation, where energy is conserved.

Now, for our problem, we need to acquire both the non-conservative and conservative forms of the governing equations that describe the physics of a 1D quasi nozzle flow. As derived from the JDA book, we can now write the non-conservative form of the governing equations such as,

The continuity equation:

`(deltarho\')/(deltat\') = -(rho\'deltaV\')/(deltax\') - (rho\'V\'(delta(lnA\'))/(deltax\')) - (V\'deltarho\')/(deltax\')`

The momentum equation:

`(deltaV\')/(deltat\') = -(V\'deltaV\')/(deltax\') - (1/gamma)[(deltaT\')/(deltax\') +(T\')/(rho\')(deltarho\')/(deltax\')]`

The energy equation:

`(deltaT\')/(deltat\') = -(V\'deltaT\')/(deltax\') - (gamma-1)T\'[(deltaV\')/(deltax\') + V\'(delta(lnA\'))/(deltax\')]`

where, all the physical terms are non-dimensional quantities.

`(rho^ast)/(rho_0)` = `rho\'` = density

`(V^ast)/(V_0)` = `V\'` = velocity

`(A^ast)/(A_0)` = `A\'` = Area

`(T^ast)/(T_0)` = `T\'` = Temperature

`(x^ast)/(x_0)` = `x\'` = length of nozzle

`t^ast/t_0 = t\'` = Time taken

and, `rho^ast, T^ast` etc., are values at the throat section and `rho_0, T_0` etc., are initial values.

These equations can be directly represented in code using MacCormack\'s method. We use finite volume method to derive these equations from their integral form and then use finite difference techniques in discretizing them with the MacCormack\'s method.

Initial conditions are specifically given as a function of time and nozzle length, which complement the physics of the given problem. Therefore, initial conditions must be chosen in a fashion that qualitatively behaves in the same way as flow parameters change in the physical system. In this case, we know that density and temperature decreases and velocity increases inside the nozzle. Hence, let us assume linear variations of the flow parameters such that,

`rho = 1 - 0.3146x`

`T = 1 - 0.2314x`

`V = (0.1+1.09x)T^(1/2)`

Area is also given as a function of nozzle length, depicting a parabolic area distribution,

`A = 1+2.2(x-1.5)^2`, for `0<= x<= 3`

Boundary conditions are given along with the initial conditions as,

Inlet:

`V(1) = 2*V(2) - V(3)`

Outlet:

`V(n) = 2*V(n-1)-V(n-2)`

`rho(n) = 2*rho(n-1)-rho(n-2)`

`T(n) = 2*T(n-1)-T(n-2)`

Conservative form of equations:

Similarily, the conservative form of governing equations are given as,

The continuity equation:

`(delta(rho\'A\'))/(deltat\') = -(delta(rho\'A\'V\'))/(deltax\')`

The momentum equation:

`(delta(rho\'A\'V\'))/(deltat\') = -(delta[(rho\'A\'V\'^2)+(1/gamma)rho\'A\'])/(deltax\') + (1/gamma) p\'(deltaA\')/(deltax\')`

The energy equation:

`(delta[rho\'((T\')/(gamma-1) + (gamma/2)V\'^2)A\'])/(deltat\') = -(delta[rho\'((T\')/(gamma-1) + (gamma/2)V\'^2)A\'+rho\'A\'V\'])/(deltax\')`

As we can see, we cannot directly use the primitive variables in these equations as code to solve the problem, and hence we use the solution vectors, that can be written as functions of those dependant variables itself.

Now, we can redefine our equations in the form of these solution vectors as,

The continuity equation:

`(deltaU_1)/(deltat\') = -(deltaF_1)/(deltax\')`

The momentum equation:

`(deltaU_2)/(deltat\') = -(deltaF_2)/(deltax\') + J_2`

The energy equation:

`(deltaU_3)/(deltat\') = -(deltaF_3)/(deltax\')`

Remember that for the conservation form, we solve with respect to the solution vectors and then after solving, we update the primitive variables, which are directly derived from the solution vectors itself such as,

Similarily, the pure form of the solution vectors and flux terms are also given as,

Initial conditions are also given, keeping in mind the dependant variables with the solution vectors,

The primitive variables are updated again and then the boundary conditions are given as,

Inlet:

`U_1(1) = rho(1)*A(1)`

`U_2(1) = 2*U_2(2) - U_2(3)`

`V(1) = (U_2(1))/(U_1(1))`

`U_3(1) = U_1(1)*[(t(1))/(gamma-1)+[(gamma/2)*(V(1)^2)]]`

Outlet:

`U_1(n) = 2*U_1(n-1)-U_1(n-2)`

`U_2(n) = 2*U_2(n-1)-U_2(n-2)`

`U_3(n) = 2*U_3(n-1)-U_3(n-2)`

MACCORMACK METHOD

MacCormack method is a method which also imploys time marching procedure. It consists of a flow of process which can be defined as:

1. Initial conditions
2. Predictor step
3. Solution update
4. Corrector step
5. Average calculation
6. Time marching update

Let us assume that a certain function, f, is used to represent the flow of computation across the MacCormack method. Function f is a governing equation of a specific physical process.

• The initial conditions are first defined as per the problem
• The predictor step is the step in which the PDE\'s are discretized using the forward differencing method(FWD).

`(deltaf)/(deltat_p) =` Forward differencing of space terms (In our case of 1D quasi nozzle flow, the space term is just the \'x\' co-ordinate).

• The solution update is the step where the variables solved using predictor step are updated.

`f^ast = f_(0) + (deltaf)/(deltat_p)*Deltat`

where,

`f_0` = previous or initial value of solution function.

`Deltat` = Time step defined.

• The corrector step is the step in which the PDE\'s are discretized using the backward differencing method(BWD).

`(deltaf^ast)/(deltat_c) = Backward differencing of space terms (In our case of 1D quasi nozzle flow, the space term is just the \'x\' co-ordinate).

As we can see, the updated values are used in the corrector step.

• The average of both the predictor and corrector step is calculated.
• The time marching update is where the next time step value is calculated by adding the initial value with the average value calculated in the last step, multiplied by the time step value defined.

`f^(n+1)= f_(0) + (deltaf)/(deltat_(avg))*Deltat`

This type of discretization to compute the solution is known as MacCormack\'s method.

Now, we can start coding for the problem.

We have to account for and analyze, by:

• Writing separate functions for conservative and non-conservative forms
• Compare the normalized mass flow rate between the conservative forms and non-conservative forms - make suitable conclusions
• Figure out the minimum number of cycles for which the simulation should be run in order for convergence
• Which method is faster?

CODING: MAIN PROGRAM

%% Simulation of 1D Supersonic nozzle flow using Macormack method
% Implementing both non-conservative and conservative forms of governing
% equations
clear all
close all
clc

% inputs
n = 31; % this is given to obtain specific initial values according to JDA book
x = linspace(0,3,n);
dx = x(2)-x(1);
gamma = 1.4;
th = ((n-1)/2)+1; % Throat value definition along length of nozzle

% specific initial values for non-conservative form
rho = 1-0.3146*x; % density
t = 1-0.2314*x; % here t is temperature
v = (0.1+1.09*x).*t.^0.5; % velocity
a = 1+2.2*(x-1.5).^2; % area
C = 0.5; % arbitrary constant Courant number for time step manipulation
m_exact = linspace(0.579,0.579,n); % Exact mass flow rate is constant across the nozzle length
% all these values are non dimensionalized

% time steps
nt = 1400;
nta = 1:nt; % Time steps in an array

% specific initial values for conservative form
for i=1:n
    if (x(i)>=0 && x(i)<=0.5)
        rho_c(i) = 1;
        t_c(i) = 1;
    end
    
    if (x(i)>=0.5 && x(i)<=1.5)
    rho_c(i) = 1-0.366*(x(i)-0.5); 
    t_c(i) = 1-0.167*(x(i)-0.5); 
    end
    
    if (x(i)>=1.5 && x(i)<=3.5)
        rho_c(i) = 0.634-0.3879*(x(i)-1.5); 
        t_c(i) = 0.833-0.3507*(x(i)-1.5); 
    end
end
v_c = 0.59./(rho_c.*a); %0.59 because it is closer to exact Mass flow rate = 0.579
% all these values are non dimensionalized

%% Calling functions
tic
[density_noncon,temp_noncon,mach,pr,m,mach_nt,pr_nt,t_nt,rho_nt] = noncon(n,nt,C,rho,v,t,a,dx,gamma,th);
st_noncon = toc;

tic
[density_cons,temp_cons,mach_c,pr_c,m_c,machc_nt,prc_nt,tc_nt,rhoc_nt] = cons(n,nt,C,rho_c,v_c,t_c,a,dx,gamma,th);
st_cons = toc;

%% Plot for comparing flow parameters
figure(1)
subplot(4,1,1)
plot(x,mach)
hold on
plot(x,mach_c,\'*\')
hold off
grid on
xlabel(\'Length of nozzle\')
ylabel(\'Mach number\')
title(\'Variation of flow parameters: Mach number, Density, Pressure, Temperature\')
legend(\'Non-Conservative form\',\'Conservative form\')

subplot(4,1,2)
plot(x,density_noncon)
hold on
plot(x,density_cons,\'*\')
hold off
grid on
xlabel(\'Length of nozzle\')
ylabel(\'Density\')
legend(\'Non-Conservative form\',\'Conservative form\')

subplot(4,1,3)
plot(x,pr)
hold on
plot(x,pr_c,\'*\')
hold off
grid on
xlabel(\'Length of nozzle\')
ylabel(\'Pressure\')
legend(\'Non-Conservative form\',\'Conservative form\')

subplot(4,1,4)
plot(x,temp_noncon)
hold on
plot(x,temp_cons,\'*\')
hold off
grid on
xlabel(\'Length of nozzle\')
ylabel(\'Temperature\')
legend(\'Non-Conservative form\',\'Conservative form\')

%% Plot for comparing mass flow rates
figure(2)
plot(x,m)
hold on
plot(x,m_c)
hold on
plot(x,m_exact)
hold off
grid on
xlabel(\'Length of nozzle\')
ylabel(\'Mass flow rate\')
legend(\'Non-Conservative form\',\'Conservative form\',\'Exact analytical value\')
title(\'Normalized mass flow rate comparison\')

%% Plot for comparing simulation time
figure(3)
st = [st_noncon st_cons];
bar(st)
grid on
xlabel(\'Comparison of Simulation time\')
ylabel(\'Simulation time\')
labels = {\'Non-Conservative form\';\'Conservative form\'};
set(gca,\'xticklabel\',labels);
title(\'Time taken for Simulation: Which is faster?\')

%% Plot for finding cycles needed for convergence
figure(4)
subplot(4,1,1)
plot(nta,mach_nt)
hold on
plot(nta,machc_nt)
hold off
grid on
xlabel(\'Time steps taken\')
ylabel(\'Mach number\')
title(\'Convergence of flow parameters along the timesteps: Minimum cycles required\')
legend(\'Non-Conservative form\',\'Conservative form\')

subplot(4,1,2)
plot(nta,rho_nt)
hold on
plot(nta,rhoc_nt)
hold off
grid on
xlabel(\'Time steps taken\')
ylabel(\'Density\')
legend(\'Non-Conservative form\',\'Conservative form\')

subplot(4,1,3)
plot(nta,pr_nt)
hold on
plot(nta,prc_nt)
hold off
grid on
xlabel(\'Time steps taken\')
ylabel(\'Pressure\')
legend(\'Non-Conservative form\',\'Conservative form\')

subplot(4,1,4)
plot(nta,t_nt)
hold on
plot(nta,tc_nt)
hold off
grid on
xlabel(\'Time steps taken\')
ylabel(\'Temperature\')
legend(\'Non-Conservative form\',\'Conservative form\')

NON-CONSERVATIVE FORM FUNCTION

function [density_noncon,temp_noncon,mach,pr,m,mach_nt,pr_nt,t_nt,rho_nt] = noncon(n,nt,C,rho,v,t,a,dx,gamma,th)

% Calculating local time step
    for b = 1:n
    dt(b) = C*dx/(sqrt(t(b))+v(b));
    end
    dt = min(dt);

% Outer time loop

for k=1:nt
    
    rho_old = rho;
    vold = v;
    told = t;
    aold = a;
    
    % Non - Conservative equations
    % Predictor method - use FWD

    % First define and discretize term wise variables
    for j=2:n-1
       dvdx = (v(j+1)-v(j))/dx;
       drhodx = (rho(j+1)-rho(j))/dx;
       dlogadx = (log(a(j+1))-log(a(j)))/dx;
       dtdx = (t(j+1)-t(j))/dx;
       
       % Continuity equation
       drhodt_p(j) = - rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
       
       % Momentum equation
       dvdt_p(j) = - v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
       
       % Energy equation
       dtdt_p(j) = - v(j)*dtdx - (gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
       
       % Solution update
       v(j) = v(j) + dvdt_p(j)*dt;
       rho(j) = rho(j) + drhodt_p(j)*dt;
       t(j) = t(j) + dtdt_p(j)*dt;
    end
    
    % Corrector method - use BWD

    % First define and discretize term wise variables
    for j=2:n-1
       dvdx = (v(j)-v(j-1))/dx;
       drhodx = (rho(j)-rho(j-1))/dx;
       dlogadx = (log(a(j))-log(a(j-1)))/dx;
       dtdx = (t(j)-t(j-1))/dx;
       
       % Continuity equation
       drhodt_c(j) = - rho(j)*dvdx - rho(j)*v(j)*dlogadx - v(j)*drhodx;
       
       % Momentum equation
       dvdt_c(j) = - v(j)*dvdx - (1/gamma)*(dtdx + (t(j)/rho(j))*drhodx);
       
       % Energy equation
       dtdt_c(j) = - v(j)*dtdx - (gamma-1)*t(j)*(dvdx + v(j)*dlogadx);
    end
    
    % Calculating next time step using average of predictor and corrector
    
    drhodt = 0.5*(drhodt_p + drhodt_c);
    dvdt = 0.5*(dvdt_p + dvdt_c);
    dtdt = 0.5*(dtdt_p + dtdt_c);
    
    % Final solution update
    for i=2:n-1
        rho(i) = rho_old(i) + drhodt(i)*dt;
        v(i) = vold(i) + dvdt(i)*dt;
        t(i) = told(i) + dtdt(i)*dt;
    end
    
    % Calculation
    mach = v./(t.^0.5);
    pr = rho.*t;
    m = rho.*a.*v;
    % Calculation of specific values at each time step along the throat
    mach_nt(k) = mach(th);
    pr_nt(k) = pr(th);
    rho_nt(k) = rho(th);
    t_nt(k) = t(th);
    
    % Apply boundary conditions
    % Inlet
    v(1) = 2*v(2) - v(3);
    
    % Outlet
    v(n) = 2*v(n-1)-v(n-2);
    rho(n) = 2*rho(n-1)-rho(n-2);
    t(n) = 2*t(n-1)-t(n-2);
end
velocity_noncon = v;
density_noncon = rho;
temp_noncon = t;
end

CONSERVATIVE FORM FUNCTION

function [density_cons,temp_cons,mach_c,pr_c,m_c,machc_nt,prc_nt,tc_nt,rhoc_nt] = cons(n,nt,C,rho_c,v_c,t_c,a,dx,gamma,th)

% Calculating local time step
    for b = 1:n
    dt(b) = C*dx/(sqrt(t_c(b))+v_c(b));
    end
    dt = min(dt);
    
    u1 = rho_c.*a;
    u2 = rho_c.*a.*v_c;
    u3 = u1.*((t_c/(gamma-1))+((gamma/2).*(v_c.^2)));
    pr_c = rho_c.*t_c;

% Outer time loop

for k=1:nt
    
    u1old = u1;
    u2old = u2;
    u3old = u3;
    
    % Conservative equations
    % Predictor method - use FWD

    % First define and discretize term wise variables
    for i=1:n
    f1(i) = u2(i);
    j2(i) = (((gamma-1)/gamma)*(u3(i)-((gamma/2)*((u2(i)^2)/u1(i)))));
    f2(i) = (u2(i)^2/u1(i)) + j2(i);
    f3(i) = ((gamma*u2(i)*u3(i))/u1(i))-(gamma*(gamma-1)*(u2(i)^3))/(2*(u1(i)^2));
    end
    
    for j=2:n-1
       % Discretize the derivatives 
       df1dx = (f1(j+1)-f1(j))/dx;
       df2dx = (f2(j+1)-f2(j))/dx;
       df3dx = (f3(j+1)-f3(j))/dx;
       dadx = (a(j+1)-a(j))/dx;
       j3(j) = (1/gamma)*(rho_c(j)*t_c(j))*dadx;
       
       % Continuity equation
       du1dt_p(j) = - df1dx;
       
       % Momentum equation
       du2dt_p(j) = - df2dx + j3(j);
       
       % Energy equation
       du3dt_p(j) = - df3dx;
       
       % Solution update
       u1(j) = u1(j) + du1dt_p(j)*dt;
       u2(j) = u2(j) + du2dt_p(j)*dt;
       u3(j) = u3(j) + du3dt_p(j)*dt;
    end
    
    % Primitive variables update
    rho_c = u1./a;
    v_c = u2./u1;
    t_c = (gamma-1)*((u3./u1)-((gamma*(v_c.^2))/2));
    
    % First redefine term wise variables before corrector step
    for i=1:n
    f1(i) = u2(i);
    j2(i) = (((gamma-1)/gamma)*(u3(i)-((gamma/2)*((u2(i)^2)/u1(i)))));
    f2(i) = (u2(i)^2/u1(i)) + j2(i);
    f3(i) = ((gamma*u2(i)*u3(i))/u1(i))-(gamma*(gamma-1)*(u2(i)^3))/(2*(u1(i)^2));
    end
    
    % Corrector method - use BWD

    % First define and discretize term wise variables
    for m=2:n-1
       % Discretize the derivatives 
       df1dx = (f1(m)-f1(m-1))/dx;
       df2dx = (f2(m)-f2(m-1))/dx;
       df3dx = (f3(m)-f3(m-1))/dx;
       dadx = (a(m)-a(m-1))/dx;
       j3(m) = (1/gamma)*(rho_c(m)*t_c(m))*dadx;
       
       % Continuity equation
       du1dt_c(m) = - df1dx;
     
       % Momentum equation
       du2dt_c(m) = - df2dx + j3(m);
       
       % Energy equation
       du3dt_c(m) = - df3dx;
    end
    
    % Calculating next time step using average of predictor and corrector
    
    du1dt = 0.5*(du1dt_p + du1dt_c);
    du2dt = 0.5*(du2dt_p + du2dt_c);
    du3dt = 0.5*(du3dt_p + du3dt_c);

    
    % Final solution update
    for i=2:n-1
        u1(i) = u1old(i) + du1dt(i)*dt;
        u2(i) = u2old(i) + du2dt(i)*dt;
        u3(i) = u3old(i) + du3dt(i)*dt;
    end
    
    % Apply boundary conditions
    % Inlet
    u1(1) = rho_c(1)*a(1);
    u2(1) = 2*u2(2) - u2(3);
    v_c(1) = u2(1)/u1(1);
    u3(1) = u1(1)*((t_c(1)/(gamma-1))+((gamma/2)*(v_c(1)^2))); 
    
    % Outlet
    u1(n) = 2*u1(n-1)-u1(n-2);
    u2(n) = 2*u2(n-1)-u2(n-2);
    u3(n) = 2*u3(n-1)-u3(n-2);
    
    % Final variables update
    rho_c = u1./a;
    v_c = u2./u1;
    t_c = (gamma-1).*((u3./u1)-((gamma.*(v_c.^2))/2)); 
    
    % Calculation of specific terms
    m_c = rho_c.*a.*v_c;
    mach_c = real(v_c./(t_c.^0.5));
    pr_c = rho_c.*t_c;
    % Calculation of specific values at each time step along the throat
    machc_nt(k) = mach_c(th);
    prc_nt(k) = pr_c(th);
    rhoc_nt(k) = rho_c(th);
    tc_nt(k) = t_c(th);
end
velocity_cons = v_c;
density_cons = rho_c;
temp_cons = t_c;
end

RESULTS AND INFERENCE

For accuracy, we have utilized a Courant number of 0.5 and calculated the time step according to stability considerations using the formula,

`Deltat = C * (Deltax)/(sqrt(T) + V)`

We could also use the array values with the given grid specifications, n =31, and take the minimum value of time step for the analysis.

Now, let\'s compare and analyze results.

Variation of flow parameters: Conservative vs Non-Conservative form

We have got the desired results from both forms of equations and as we can clearly see, there isn\'t much of a difference. But the consistency and accuracy of the simulation will be clear with other parameter calculations too, like mass flow rate.

Mass flow rate comparison:

From the above graph, the oscillation of values until attaining convergence is very clear and the conservative form has a more consistent and accurate profile, nearing the exact analytical value of the problem, than the non-conservative form. The latter has lot more oscillations in the start and end of the nozzle and values are not stable for an accurate representation of mass flow rate. Hence, conservation form has a clear advantage in this perspective.

Convergence-Minimum cycles required:

Another proof of stability of the conservation form is given by the fact that iteration of time steps needed for convergence is also relatively lower compared to the non-conservative form, which again has lot more oscillations in flow property values. For example, the flow properties of mach number for non-conservative form, converge almost around the 400-500 time step range after initial oscillations, whereas for the conservative form, it converges around 200-300 time step range itself. Minor oscillations and a consistent profile provides an overall edge to the conservative form of equations, numerically.

Simulation time:

As discussed in various problems, simulation time is a factor of computational load and also the speed of the computer processing the simulation. Hence, as a direct effect of more computational load, the conservative form marginally takes higher simulation time than the non-conservative form, which has relatively lower computational load and has no solution vectors or flux variables to update and account for. But considering the minimum cycles required for convergence and overall consistency and accuracy, the conservative form of equations provide a more cost-effective and efficient way of solving problems than the non-conservative form of equations.

But in a broader sense, both the forms of equations provide accurate results and a clear superiority cannot be established with just one physical problem. Mathematically speaking, both sets of equations are the same and they vary marginally with regards to numerical schemes and methods employed.

GRID DEPENDANCE:

While working with the values, we could also see whether the increase in grid changes the accuracy of the system or not. According to the JDA book, the values are taken non-dimensionally and the results are drawn out after increasing grid size to 61.

As the values indicate, a marginal change in values for the non-conservation part is visible but still it\'s in the order of 1e-2 and hence can be neglected. As far as the conservative part, it almost shows a kind of independance from the grid size and values are unchanged. They are closer to analytical values but still a grid independance could be infered, even after doubling the grid points.