Week 4 - Linear Convection

Effect of time step size on the solution:

Objective:

• To represent the wave equcation for convection @ t = .4 Second for time step .1. .01. .001, .0001 in single plot.
• Ploting of simulation time for different time step

The body of the program (time_step_dt) call for the function time_step_loop for different input values of the time step dt. the function gives output for velocity (u) @ different points (x).

To calculate the simulation time we can use in build option tic , toc in Matlab.

The wave equcation can be be obtain for all time step but stability will be changes

Stability for the simple linear convection problem by explicit method can be find by CFL (Courant Freidrich Lewy) condition

`mu = a*((Deltax)/(Delta2))(1-(a*(Deltat)/(Deltax)))`

`mu =`numerical diffusion

Courant Number, `CFL = (C*(Deltat)/(Deltax))`

For a stable solution CFL <= 1

In this linear convection problem we are taking following parameters

L = 1
n = 80
x = linspace(0,L,n)
C = 1
dx = x(2) - x(1)
nt = .4

@ dt = .1 , CFL = 7.8740 which is to much greater then 1. hence solution become unstable and wigles are form as shown in plot 1 by blue line, this is called as gibbs phenomeon. Solution never converges and wigles will form continuesly leading to unstable solution.

@dt = .01, CFL = .7874 approx to 1. However, solution is stable but accuracy is very low. This leads to trapezium wave form as shown by orange line.

@dt = .001, CFL = .0787 is much less then 1. solution become stable but approx equcal to square wave form as shown by yellow line.

@dt = .0001, CFL = .0079 approx to 1. solution is stable and accuracy is also high. This leads to almost square wave form as shown by purple line.

Plot 2 shows the simulation time for different time step (dt)

Main body of program

clear all close allclc tic[x0,u0] = time_step_loop(0.1)simulation_time_t1 = toctic[x1,u1] = time_step_loop(0.01)simulation_time_t2 = toctic[x2,u2] = time_step_loop(0.001)simulation_time_t3 = toctic[x3,u3] = time_step_loop(0.0001)simulation_time_t4 = tocfigure(1)plot(x0,u0,'linewidth',1.5)hold on;plot(x1,u1,'linewidth',1.5)hold on;plot(x2,u2,'linewidth',1.5)hold on;plot(x3,u3,'linewidth',1.5)hold on;axis([0 1 0 2.5])legend('velocity t=.4 @ dt=.1','velocity t=.4 @ dt=.01','velocity t=.4 @ dt=.001','velocity t=.4 @ dt=.0001')grid on;title('velocity profile @ t=.4 for dt')ylabel('velocity')xlabel('length of domain')t = [simulation_time_t1,simulation_time_t2,simulation_time_t3,simulation_time_t4]figure(2)bar3(t,1,'b')zlabel('simulation time')ylabel('time steps dt')title('Simulation time VS time step')

Sub function: time_step_loop

function [x,u] = time_step_loop(dt)L = 1n = 80x = linspace(0,L,n)C = 1dx = x(2) - x(1)nt = .4%velocity @ start & end% x_start  = .1% x_end = .3%Locattion of start & end velocity profile% i = n_start% j = n_end[m,i] = min(abs(x-.1))[s,j] = min(abs(x-.3))n_last = find(x==1)%Intailiazing array and applying BCu = ones(1,n)u(1:i) = 1 %velocity from node 1 to ith nodeu(i:j) = 2 %velocity from ith node to jth nodeu(j:n_last) = 1 % velocity from jth node to last node% u_intial = u;uold = u;% time step nt%for a = 1:length(dt)for l = 0.1:nt for k = 2:n    u(k) = uold(k)-(C*dt(1)/dx)*(uold(k)-uold(k-1)); end    uold = u;endend

Plot 1

Plot 2: simulation time for different dt