Solve 1D Burgers Equation via Semi-Implicit Adams-Bashforth/Crank-Nicolson Scheme
7 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(1 个)
Mat Hunt
2023-7-14
Make sure your vaiue of D*dt/dx^2 isn't too large, ideally, it should be around 1.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Boundary Conditions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!