periodic boundary conditions for advection for a<0

9 views (last 30 days)
lulu
lulu on 20 Dec 2020
Commented: lulu on 21 Dec 2020
The following shows the codes for advection of population move to right and left . I got the correct BCs for right moving that I am able to see the solution profile leaving the domain at one end, and re-entering at the other end. Whereas, I am struggled to get a correct BCs for left moving. Could you please help.
clear all;
% numerical simulation for ur_t+a*ur_x=0 and ul_t-a*ul_x=0 where ur: population moving to right, ul: population moving to left
xl=0; xr=1; %domain[xl,xr]
J=1000; % J: number of dividion for x
dx=(xr-xl)/J; % dx: mesh size
tf=50; % final simulation time
Nt=10000; %Nt:number of time steps
dt=tf/Nt;
c=50; % paremeter for the solution
lambdaR=0.005; %turning rate
lambdaL=0.003; %turning rate
a=0.1; % speed
mu=a*dt/dx;
if mu>1.0 %make sure dt satisfy stability condition
error('mu should<1.0!')
end
%Evaluate the initial conditions
x=xl:dx:xr; %generate the grid point
f=exp(-c*(x-0.5).^2); %dimension f(1:J+1)
%store the solution at all grid points for all time steps
ur=zeros(J+1,Nt);
ul=zeros(J+1,Nt);
%find the approximate solution at each time step
for n=1:Nt
t=n*dt; % current time
if n==1 % first time step
for j=2:J % interior nodes
% u(j,n)=f(j)- mu *(f(j) - f(j-1)); %upwind scheme
ur(j,n)=f(j)-0.5*mu*(f(j+1)-f(j-1))+0.5*mu^2*(f(j+1)-2*f(j)+f(j-1))+lambdaR*f(j)-lambdaL*f(j); %lax wendroff
ul(j,n)=f(j)+0.5*mu*(f(j+1)-f(j-1))+0.5*mu^2*(f(j+1)-2*f(j)+f(j-1))-lambdaR*f(j)+lambdaL*f(j);
end
ur(J+1,1)=ur(J,1); %peridic BC for right moving
ur(1,1)=ur(J+1,1); %peridic BC for right moving
ul(J+1,1)=ul(1,1); %peridic BC for left moving
ul(J,1)=ul(J-1,1); %peridic BC for left moving
else
for j=2:J % interior nodes
% u(j,n)=u(j,n-1)- mu *(u(j,n) - u(j-1,n)); %upwind scheme
ur(j,n)=ur(j,n-1)-0.5*mu*(ur(j+1,n-1)-ur(j-1,n-1))+0.5*mu^2*(ur(j+1,n-1)-2*ur(j,n-1)+ur(j-1,n-1))+lambdaR*ul(j,n-1)-lambdaL*ur(j,n-1);
ul(j,n)=ul(j,n-1)+0.5*mu*(ul(j+1,n-1)-ul(j-1,n-1))+0.5*mu^2*(ul(j+1,n-1)-2*ul(j,n-1)+ul(j-1,n-1))-lambdaR*ul(j,n-1)+lambdaL*ur(j,n-1);
end
ur(J+1,n)=ur(J,n); %peridic BC for right moving
ur(1,n)=ur(J+1,n); %peridic BC for right moving
ul(J+1,n)=ul(1,n); %peridic BC for left moving
ul(J,n)=ul(J-1,n); %peridic BC for left moving
end
%calculate the analytical solution
for j=1:J+1
xj=xl+(j-1)*dx;
u_ex(j,n)=exp(-c*(xj-a*t-0.5)^2);
end
end
%plot the analytical and numerical solution at different times
figure;
hold on;
n=1;
plot(x,ur(:,n),'r-');
plot(x,ul(:,n),'r--');
%plot(x,ur(:,n),'r-',x,ul(:,n),'k-',x ,u_ex(:,n),'r.'); %red
n=1000;
plot(x,ur(:,n),'g-');
plot(x,ul(:,n),'g--');
%plot(x,ur(:,n),'g-',x,ul(:,n),'b-',x, u_ex(:,n),'g.'); %green
n=1600;
plot(x,ur(:,n),'b-'); %blue
plot(x,ul(:,n),'b--');
n=2000;
plot(x,ur(:,n),'m-');
plot(x,ul(:,n),'m--');
n=3000;
plot(x,ur(:,n),'k-'); %black
plot(x,ul(:,n),'k--'); %black
% %%% time=n*\deta t=
legend('Numericalrght-moving t=0','Numericalleft-moving t=0','Numericalrght-moving t=5','Numericalleft-moving t=5','Numericalrght-moving t=8','Numericalleft-moving t=8','Numericalrght-moving t=10','Numericalleft-moving t=10','Numericalrght-moving t=15','Numericalleft-moving t=15');% title('Numerical and Analytical Solutions at t=0.1,0.3,0.5');

Answers (1)

Alan Stevens
Alan Stevens on 20 Dec 2020
For left moving instead of
ul(J,n)=ul(J-1,n); %peridic BC for left moving
shouldn't you have
ul(J,n)=ul(J+1,n); %peridic BC for left moving
assuming J+1 is to the right of J.
  5 Comments
lulu
lulu on 21 Dec 2020
Hello, turning rate which allow to turn right and left. Where the flow of right moving> the flow of left moving.

Sign in to comment.

Categories

Find more on Graphics Object Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by