Solve 1D Burgers Equation via Semi-Implicit Adams-Bashforth/Crank-Nicolson Scheme
7 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hello world,
I'm trying to solve the 1D Burgers equation (nonlinear convection-diffusion equation) by applying the explicit Adams-Bashforth scheme to the nonlinear convective term and the implicit Crank-Nicolson scheme to the diffusive term. The equation in dimensionless form is
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804344/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/804349/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804354/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804359/image.png)
The exact solution is
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804364/image.png)
so that the initial condition and steady state solutions respectively are
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804369/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804374/image.png)
Note that there is a discontinuity in the exact solution at
. To avoid this, we can start the solution at
where
is a small parameter on the order of
or
. Thus, the initial condition becomes
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804379/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804384/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804389/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804394/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804399/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804404/image.png)
I wish to construct a solution using 100 spatial steps in
so that
, and 278 time steps in
so that
. This makes the Courant number (convective stability)
where
since
is the maximum velocity at
found by a perturbation of
. Note that there is no restriction on viscous stability here since the plan is to use an implicit method on the diffusion term.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804349/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804414/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804419/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804424/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804429/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804434/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804439/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804444/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804384/image.png)
I have solved the 1D diffusion equation using the implicit Crank-Nicolson method, I have solved the 1D Burgers equation using the explicit Euler method, and I have solved several ODEs using the explicit two-step Adams-Bashforth method. I think I understand the theory of this problem, and on paper I have written out the time discretized equations, but my skill in MATLAB is novice and I have no idea how to write code that combines these methods to solve a single PDE. Here, since I'm trying to use the Crank-Nicolson method on the diffusive term in the Burgers equation, I suspect that I can just apply my code from the diffusion equation (in conjunction with some unwritten code for the Adams-Bashforth method) to these initial and boundary conditions.
The MATLAB code from my Crank-Nicolson solution to the 1D diffusion equation applied directly to the initial and boundary conditions of the 1D Burgers equation is:
clear all; close all; clc;
% Construct spatial mesh
Nx = 100; % Number of spatial steps
xl = -6; % Left x boundary
xr = 6; % Right x boundary
dx = (xr-xl)/Nx; % Spatial step
x = xl:dx:xr; % x
% Construct temporal mesh
tf = 2; % Final time
dt = 0.0072; % Timestep
Nt = round(tf/dt); % Number of timesteps
t = 0:dt:tf; % t
% Parameters
umax = 15; % Found by a perturbation of t=10^-2
C = umax*dt/dx; % Convective Stability / Courant Number
disp("Courant Number: "+C)
V = dt/(dx*dx); % Viscous Stability / Diffusion Number
disp("Diffusion Number: "+V)
% Define Boundary & Initial Conditions
% Boundary conditions (Dirichlet)
for j = 1:Nt+1 % Loop through all t
u(1,j) = 2; % u(-6,t) = 2
u(Nx+1,j) = -2; % u(6,t) = -2
end
% Initial condition (t = t0, NOT t = 0)
for i = 1:Nx+1 % Loop through all x
u(i,1) = (-2*sinh(x(i)))/(cosh(x(i))-exp(0.0001)); % u(x,t0)
end
% Define the right (MMr) and left (MMl) matrices in the Crank-Nicolson method
aal(1:Nx-2) = -V; % Below the main diagonal in MMl
bbl(1:Nx-1) = 2+2*V; % The main diagonal in MMl
ccl(1:Nx-2) = -V; % Above the main diagonal in MMl
MMl = diag(bbl,0)+diag(aal,-1)+diag(ccl,1); % Building the MMl matrix
aar(1:Nx-2) = V; % Below the main diagonal in MMr
bbr(1:Nx-1) = 2-2*V; % The main diagonal in MMr
ccr(1:Nx-2) = V; % Above the main diagonal in MMr
MMr = diag(bbr,0)+diag(aar,-1)+diag(ccr,1); % Building the MMr matrix
% Implementation of the Crank-Nicholson method
for j = 1:Nt
u(2:Nx,j+1) = inv(MMl)*MMr*u(2:Nx,j);
end
Plotting this solution:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/804454/image.png)
Clearly this is problematic, but hopefully a good start. In addition to incorporating the Adams-Bashforth scheme, I think I also need to incorporate the Neumann boundary conditions. If anyone out there is familiar with this semi-implicit AB/CN method, I would be very appreciative of some guidance. I know this is an elaborate questions, so please don't hesitate to write comments and communicate with me about anything.
Thank you!
0 commentaires
Réponses (1)
Mat Hunt
le 14 Juil 2023
Make sure your vaiue of D*dt/dx^2 isn't too large, ideally, it should be around 1.
0 commentaires
Voir également
Catégories
En savoir plus sur Boundary Conditions dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!