Linear Convection for varying grid size and understanding significance of grid size on solution stability

Objective : To calculate Ufinal">Ufinal of Linear Convection equation for set of nodes & plot Ufinal">Ufinal vs X">X

Theory :

Linear Convection part of the Navier stokes equation is given by

∂u∂t+c⋅∂u∂x=0">∂u∂t+c⋅∂u∂x=0

Indexing for 'u' to characterize space & time marching elements :

Uxn">Unx

where

n is the Time index

x is the space index

Let us consider a square wave form with nodes & specific data to apply time & space marching approach

Applying forward differencing of 2nd order accuracy to time derivative term & backward differencing of 2nd order accuracy to spatial derivative term

[Uin+1-UinΔT]=-c⋅[Uin-Ui-1nΔx]">[Un+1i−UniΔT]=−c⋅[Uni−Uni−1Δx]

∴Uin+1=Uin-(c⋅ΔTΔx)⋅(Uin-Ui-1n)">∴Un+1i=Uni−(c⋅ΔTΔx)⋅(Uni−Uni−1)

Now c⋅ΔTΔx">c⋅ΔTΔx

= k

∴Uin+1-Uin=-(k)⋅(Uin-Ui-1n)→">∴Un+1i−Uni=−(k)⋅(Uni−Uni−1)→Explicit Form

∴Uin+1-Uin=-(k)⋅(Uin+1-Ui-1n+1)→">∴Un+1i−Uni=−(k)⋅(Un+1i−Uni−1+1)→Implicit Form

Given Data & Problem to be solved:

Numerical solution to the One-dimensional Linear Convection equation.

1. Assume that the domain length is L = 1m

2. The initial velocity profile is a step function. It is equal to 2m/s between x= 0.1 and 0.3 and 1m/s everywhere else

3. Use first-order forward differencing for the time derivative

4. Use first-order backward/rearward differencing for the space term.

5. Constant convection coefficient, c = 1 m/s

6. Time step, dt = 0.01

7. Simulation run time , t = 0.4 seconds.

Perform the simulation for Number of nodes, n=20,40,80 and 160.

Matlab code to calculate U final for array of Nodes :

%% Approximating U_final with variable nodes

% Input Data
L = 1;              % Length of domain
n = [20 40 80 180]; %No of grid points
c = 1;              % Constant convection cooefficient
dt = 0.01;          % Time step 
t = 0.4;            % Integration time
n_start = 0.1;      % Starting point of change in Velocity
n_stop = 0.3;       % starting & ending points of change in Velocity
nt = t/dt;          % No of time steps

% Initial setup for X & U 

x_result = 0;
u_result = 0;

for i = 1:length(n)
    [x_result,u_result] = linear_convection_var_function(L,n(i),n_start,n_stop,nt,c,dt,i);
end

Function code to calculate U final & X :

%% Function to calculate x,u_final,n_start & n_stop for variation in nodes & other parameters

function [x,u] = linear_convection_var_function(L,n,n_start,n_stop,nt,c,dt,i)
                 
           % Calculation of Parameters
                      x  = linspace(0,L,n); % Meshing for length L
                      dx = x(2)-x(1);      % Grid Size
                 
           % Calculating starting & ending points of change in Velocity
           % Assuming values of x are in A.P with difference dx
           % n_start = 0 + (f_start - 1)*dx (nth term of an AP)
           % n_stop  = 0 + (f_stop - 1)*dx 
                 f_start = round((n_start/dx)+1);
                 f_stop  = round(n_stop/dx)+1;
                 
           % Initializing U array & setting BCs
                       u = 1*ones(1,n);
       u(f_start:f_stop) = 2;
       
           % Boundary Condition 
                    u(1) = 1;
                    uold = u;
                       u_initial  = u;
           % Time Loop
            for j  = 1:nt
                
           % Space Loop    
                for k = 2:n
                    u(k) = uold(k) - (c*(dt/dx)*(uold(k) - uold(k-1)));
                end
                
           % Update old Velocities
           
                    uold = u;
                                   
           % Plotting subplots U_final vs U-initial
           figure(i)
           subplot(2,1,1)
           plot(x,u_initial,'b');
           axis([0 1 1 2])
           xlabel('X')
           ylabel('Velocity')
           legend('U initial profile')
           grid on 
           hold on
           
           subplot(2,1,2)
           plot(x,u,'k');
           axis([0 1 1 2])
           xlabel('X')
           ylabel('Velocity')
           legend('U final profile')
           grid on 
           hold off
           title('Final Velocity profile for varying number of nodes') 
           
            end
end
    

Main Code Explaination :

• Defining input parameters for L, n, c, dt & t
• Starting & ending points of change in velocity are also mentioned by n_start & n_stop
• No of time steps nt is calculated by dividing Integration time by Time step
• Initial Setup values for X & U is set to 0
• A 'for loop' is introduced to calculate 'X & U' for set of nodes in the array

Function Code Explaination :

• A function is created with an input of L, n, c, n_start, n_stop, dt & nt and with X & U being output values
• Defining meshing for length 'L' of the domain with 'n' grid points
• Defining grid size by the difference between two consecutive mesh points
• Calculating starting & ending points of change in Velocity by using Arithmetic Progression method

Nth Term of AP :

n_start = 0 + (f_start - 1)*dx →">→

Initial Value is set to 0
n_stop = 0 + (f_stop - 1)*dx

• To avoid floats, rounding off is done by using 'round'
• Initialising U array by creating the array of Ones of 1 row & n colums
• Initialising the value of u at n_start & n_stop = 2 m/s
• Setting up the Boundary condition u(1) = 1 m/s
• Assigining U as 'uold' value & creating a duplicate called 'u_initial' for further plot reference which is also the same as U
• To perform 'nt' time steps of operation, a for loop is introduced with 'j' as a counter running through the length of 'nt'
• Calculating U final by implementing Forward differencing scheme for Time marching term & Backward diff scheme for Space Marching term

[Uin+1-UinΔT]=-c⋅[Uin-Ui-1nΔx]">[Un+1i−UniΔT]=−c⋅[Uni−Uni−1Δx]

• Updating old velocities to final velocity value
• Plotting U_initial vs U_final under subplots

Observation :

CASE I : Grid Points = 20

Interpretation : U initial value is constant at 2 m/s from 0.1 to 0.3 X . There is a clear visual of Drop down of U final value from 2 to nearly 1.7 m/s due to numerical diffusion effect which is caused due to ignorance of second order derivatives in the approximation of Taylor series.

This can be eliminated by decreasing the mesh size or by increasing nodes in the same length of domain (L).

CASE II : Grid Points = 40

Interpretation : U initial value is constant at 2 m/s from 0.1 to 0.3 X . There is a clear visual of Drop down of U final value from 2 to nearly 1.9 m/s which is comparitively less drop down than for 20 nodes, this is due to numerical diffusion effect which is caused due to ignorance of second order derivatives in the approximation of Taylor series.

This can be eliminated by decreasing the mesh size or by increasing nodes in the same length of domain (L).

CASE III : Grid Points = 80

Interpretation : U initial value is constant at 2 m/s from 0.1 to 0.3 X . There is no visual of Drop down of U final value .This is due to numerical diffusion effect is reduced completely by increasing the no of nodes in away decreasing the mesh size so that the information travels instantaneously giving accurate information.

CASE IV : Grid Points = 160

Interpretation : Distortion in the result plot is due information reaching the next neighbour mesh point (Grid point) in less than step time. This means that we reducing the length of grid size so much that information is travelling to fast that computer cannot compute it as it gets no time to calculate because the information has already passed that mesh point (Grid Point)

Understanding Stability of the solution through Courant–Friedrichs–Lewy Condition :

The Courant–Friedrichs–Lewy or CFL condition is a condition for the stability of unstable numerical methods that model convection or wave phenomena

The CFL condition expresses that the distance that any information travels during the timestep length within the mesh must be lower than the distance between mesh elements. In other words, information from a given cell or mesh element must propagate only to its immediate neighbors.

In a CFD simulation, the Courant number broadly indicates how much the information travels across a computational grid cell in a unit of time . If the Courant number is greater than one, it means that the information propagates through more than one grid cell at each time step, making the solution inaccurate and potentially leading to nonphysical results or divergence of the solution in certain integration schemes.

• For increasing no of nodes on the same length of domain, the artificial/numerical diffussion effect decreases due to decrease in grid size.

• For increase in no of nodes more than 120, it provides solution blowup which is due to CFL condition being violated. Hence appropriate mesh size or no of nodes has to limited for physically significant results