1 view (last 30 days)

Show older comments

Nx = 100;

Ny = 20;

Nt = 400;

x = 0:10/Nx:10;

y = -1:2/Ny:1;

t = 0:1/Nt:1;

Dx = 10/Nx; % delta x

Dy = 2/Ny; % delta y

Dt = 1/Nt; % delta t

%%Numerical Solution using alpha (a) = 1

% Constants

a = 1.0; % alpha

Xd = a*Dt/(Dx^2); % Diffusion Length Scale in x

Yd = a*Dt/(Dy^2); % Diffusion Length Scale in y

c = Dt/(2*Dx); % Convection Length Scale

% Boundary Conditions

T = zeros(Ny+1,Nx+1,Nt+1); % Scalar field

T(1,1:end,:) = 0; % Bottom Wall

for m = 1:Ny+1

T(m,1,:) = (1 - (y(m).*y(m))).^2; % Inlet

end

% Forward Euler in Time, Central Differencing in Space

maxiter = 500;

for k = 1:maxiter

Tlast = T; % saving the last guess

T(:,:,1) = Tlast(:,:,end); % initialize the scalar field at t = 0 to the last guess

for i = 2:Nt+1

T(2:end-1,2:end-1,i) = (1-(2*Xd)-(2*Yd))*T(2:end-1,2:end-1,i-1) + (Xd-((3/2)*c*(1-(y(2:end-1).^2)))).*T(2:end-1,3:end,i-1) + (Xd+((3/2)*c*(1-(y(2:end-1).^2)))).*T(2:end-1,1:end-2,i-1) + Yd*T(3:end,2:end-1,i-1) + Yd*T(1:end-2,2:end-1,i-1);

T(end,1:end,i) = (4*T(end-1,1:end,i) - T(end-2,1:end,i))/3; % Top Wall (dT/dy = 0)

T(1:end,end,i) = (4*T(1:end,end-1,i) - T(1:end,end-2,i))/3; % Outlet (dT/dx = 0)

end

err(k) = max(abs(T(:)-Tlast(:))); % finding the residual value between two successive iterations

if err(k) < 1E-04

break; % stopping the solution is the residual value is small

end

end

Walter Roberson
on 3 Dec 2012

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

Start Hunting!