3. (20 pts) For the 2D continuous pollutant model assume the decay rate depends

Question

3. (20 pts) For the 2D continuous pollutant model assume the decay rate depends on the spatial coordinates (x,y), ie., = (z,y). (a) Derive the corresponding finite difference model to account for this. (b) Modify the Matlab code flow2d2.m to account for the case when = 0.1 +0.05e- cos(). and use it to study this case for L = 2, W-3, m-40. n-45,T-20. mazk 400. v1- 0.15, v2 0.2 and with the initial condition u(x, y, 0) = 0.1 and the upwind boundary condi- tions u(0, y, t) = 0.1, u(z, 0, t) = 0. sinArx)| if (2, t) E [0,× [0, ) and u(z, 0, t-0.1 otherwise. Use the Matlab subplot utility to plot your results in the (r.y. u) space at the times t = 0, 5, 10, 15, 20 in a single figure. (c) If we change v from 0.15 to -0.15 in (b) while keep all other parameters the same, derive the corresponding FD model and modify the upwind boundary conditions to account for this change. Modify the Matlab code flow2d2.m to reflect the changes of the FD model and upwind boundary conditions and use the modified Matlab code to approximate the pollutant concentration u(x,y,t) and plot your results in the (r,y, u) space at the times t = 0,5, 10, 15, 20 in a single figure using subplot utility, and explain any differences you have observed from the two plots in (b) and (c).

Explanation / Answer

See modified flow2d2 function below:

-------------------------------------

function [u,x,y] = flow2d2(L,W,m,n,T,maxk,v1,v2)
%%%%%%%%%%%%%%%%%%
% compute decay lambda
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lambda = 0.1 + 0.05*exp((-10*x)/L)*cos((2*pi*y)/W);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute derived parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt = T/maxk;
dx = L/m;
dy = W/n;
%stability condition
a = 1-v1*(dt/dx) - v2*(dt/dy) - lambda*dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check stability condition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (a<=0)
fprintf(' (a,v1,v2,lambda) = (%e,%e,%e,%e),',a,v1,v2,lambda)
fprintf(' so the stability condition is NOT satisfied. ')
else
fprintf(' (a,v1,v2,lambda) = (%e,%e,%e,%e),',a,v1,v2,lambda)
fprintf(' so the stability condition is satisfied. ')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set initial conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:m+1
x(i) = (i-1)*dx;
for j = 1:n+1
y(j) = (j-1)*dy;
u(i,j,1) = 0.1;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set upwind boundry conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:maxk+1
t(k) = (k-1)*dt;
for j = 1:n+1
y(j) = (j-1)*dy;
u(0,j,k) = 0.1;
end
for i = 1:m+1
if ((x(i) > 0 & x(i) < (L/4)) & (t(k) > 0 & t(k) < (T/4))
u(i,0,k) = 0.1 + abs(sin(4*pi*x(i));
else
u(i,0,k) = 0.1;
end  
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Execute the explicit finite difference method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:maxk
for i = 2:m+1
for j=2:n+1
u(i,j,k+1) = v1*(dt/dx)u(i-1,j,k)+v2(dt/dy)u(i,j-1,k)+(1-v1(dt/dx)-v2*(dt/dy)-dt*lambda)*u(i,j,k);
end
end
end

end

--------------------------------------------

Related Questions

Navigate

• Browse (All)
• Browse 3
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.