Explicit Euler Method Solution to 1D Nonlinear Convection-Diffusion Equation

Hello world,

I'm trying to solve the 1D Nonlinear Convection-Diffusion equation (Burgers equation) using the Explicit FTCS "Euler" method. The equation has been nondimensionalized and is written as

where the second term is nonlinear. Dirichlet boundary conditions on are , while those Neumann are .

The exact solution is:

The initial and steady state solutions are:

and

The fully explicit scheme must satisfy convective stability and viscous stability , where is the maximum velocity at .

My MATLAB code so far is as follows:

clear 
close all
clc

%%% Initialize & Define Parameters

% Construct spatial mesh

Nx = 100;               % Number of spatial steps
xl = -6;                % Left boundary
xr = 6;                 % Right boundary
dx = (xr-xl)/Nx;        % Spatial step
x = xl:dx:xr;           % x 

% Construct temporal mesh

tf = 2;                 % Final time
dt = 0.007;             % Time step
Nt = round(tf/dt);      % Number of time steps
t = 0:dt:2;             % t 

% Parameters

umax = 15;              % Maximum velocity at t=0 found by a peturbation of t = 10^-2
C = umax*dt/dx;         % Convective Stability / Courant Number
disp("Convective Stability: "+C)
V = dt/(dx*dx);         % Viscous Stability / Diffusion Number
disp("Viscous Stability: "+V)

%%% Define Boundary & Initial Conditions

% Boundary conditions (Dirichlet)

for j = 1:Nt+1
    u(1,j) = 2;                % u(-6,t) = 2
    u(Nx+1,j) = -2;            % u(6,t) = -2 
end

% Initial condition

for i = 1:Nx+1
    u(i,1) = -2*sinh(x(i))/(cosh(x(i))-1);   % u(x,0) 
end

%%% Implementation of Explicit Euler Method

for j = 1:Nt                    % Time Loop
    for i = 2:Nx                % Space loop excluding boundaries
        u(i,j+1) = u(i,j) - 0.5*(dt/dx)*u(i,j)*(u(i+1,j)-u(i-1,j)) + (dt/(dx*dx))*(u(i+1,j)+u(i-1,j)-2*u(i,j));
    end
end


%%% Define Exact Solution & Select Times of Interest

tau = [0.1 0.2 0.5 1 2];                              % Times of interest

syms xi

for k = 1:length(tau)
    ua(k) = (-2*sinh(xi))/(cosh(xi)-exp(-tau(k)));    % Define analytic solutions at times of interest
    ts(k) = round(tau(k)/dt);                         % Select the indices from t corresponding to to times of interest
end

%%% Plot Solutions

figure(1)
    fplot(ua(1),[-6,6],':k','Linewidth',2);
    hold on
    plot(x,u(:,ts(1)),'-','Linewidth',1.5);  
    plot(x,u(:,ts(2)),'-','Linewidth',1.5); 
    plot(x,u(:,ts(3)),'-','Linewidth',1.5); 
    plot(x,u(:,ts(4)),'-','Linewidth',1.5);
    plot(x,u(:,ts(5)),'-','Linewidth',1.5);
    fplot(ua(2),[-6,6],':k','Linewidth',2); 
    fplot(ua(3),[-6,6],':k','Linewidth',2); 
    fplot(ua(4),[-6,6],':k','Linewidth',2);
    fplot(ua(5),[-6,6],':k','Linewidth',2)
    hold off
    title("Explicit Euler Method Solutions to 1D Convection-Diffusion Equation")
    xlabel('\it{x}')
    ylabel('\it{u(x)}')
    legend('Analytic Solutions','t = 0.1','t = 0.2','t = 0.5','t = 1','t = 2')
    grid on

With this code, the solution appears as follows:

The dotted lines represent the exact solution at times of interest and the colored lines are (should be) the Euler solutions at times of interest. The command print out is of the solution u. It looks like the solution is starting off correctly, then breaks near the origin and quickly as time elapses. I have tinkered with the dt/dx balance a lot, but this is pretty much how the solution always looks. I was very confident in my code but clearly there is a big problem. I'm at an impasse and am hoping someone out there can see the issues that I cannot.

Any advice or suggestions are welcome and all help is appreciated, thanks in advance!

Answer 1

[Adam Harris]

I solved the issue. There is a singularity in the initial condition so I needed to start the solution from a small time , this avoids the singularity at . Thus, the initial condtion is:

for i = 1:Nx+1
    u(i,1) = -2*sinh(x(i))/(cosh(x(i))-exp(-0.0001));   % u(x,t0)
end

And the solution is:

[Bjorn Gustavsson]

Well maybe you might have done. However your equation is not a convection-diffusion equation the way it is written.

The standard form of the CDE for some quantity C that convects with velocity u and diffuses with diffusion-coefficient D is:

In 1-D this simplifies, without sources or losses constant diffusion-coefficient and convection-velocity, to

Here we note that there are no explicit nonlinearities. Your equation as written would be interpreted as a "convection-diffusion-equation for the convection-velocity" which makes no sense, for the fluid velocity you would have to solve the fluid momentum-equation. Somewhere you have mixed together the convection-velocity and the convecting-diffusing quantity. Before your equation makes sense you'll have to untangle that.

HTH

[Adam Harris]

@Bjorn Gustavsson

[Bjorn Gustavsson]

I've never seen an equation of motion with a diffusion-term before , so I learnt something new - again.

[Adam Harris]

@Bjorn Gustavsson This was new to me as well! Perhaps you'd be interested in my next problem: https://www.mathworks.com/matlabcentral/answers/1589549-solve-1d-burgers-equation-via-semi-implicit-adams-bashforth-crank-nicolson-scheme

[Wahyu Widhi]

I have 2 questions, first question is what is your reference on doing this? Second, how do you change the boundary conditions into Neumann boundary condition? Thank you!