Explicit Euler Method Solution to 1D Nonlinear Convection-Diffusion Equation
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Adam Harris
am 16 Nov. 2021
Kommentiert: Wahyu Widhi
am 2 Jan. 2022
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801669/image.png)
where the second term is nonlinear. Dirichlet boundary conditions on
are
, while those Neumann are
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801674/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801679/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801684/image.png)
The exact solution is:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801689/image.png)
The initial and steady state solutions are:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801694/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801699/image.png)
The fully explicit scheme must satisfy convective stability
and viscous stability
, where
is the maximum velocity at
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801704/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801709/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801714/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801719/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/801724/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/802754/image.png)
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!
0 Kommentare
Akzeptierte Antwort
Adam Harris
am 17 Nov. 2021
Bearbeitet: Adam Harris
am 17 Nov. 2021
5 Kommentare
Wahyu Widhi
am 2 Jan. 2022
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!
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Boundary Conditions finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!