Hi. I am attempting to model natural convection in a parallelogram enclosure. I have a coordinate transformation and I have set the code to solve the temperature, stream, velocitu and Vorticity as shown below. I do not understand where the source of the NaN values arise. The teperature, vorticity, velocity and stream equation are functions of each other and an initial temperature(Poisson equation is assumed) to enable fluid computation. The code is attached below.
clear;
clc;
thetad = 90;
thetar = (pi()/180)*thetad;
L = 1;
H = 2;
nx = 50;
ny=nx;
dx = L/(nx+1);
dy = H/(ny+1);
x = 0:dx:L;
y = 0:dy:H;
nE = nx; %Chosen in dimensionless domain
nN = ny;
E = (x-y.*cot(thetar))/L;
N = (y./(H*sin(thetar)));
dE = E(2)-E(1); %Node spacing in the Dimensionless domain
dN = N(2)-N(1);
Gr = 100000;
Pr = 0.70;
%% Initial vectors
vp = zeros(nx+2,ny+2);
sp = zeros(nx+2,ny+2);
tpi = zeros(nx+2,ny+2);
u = zeros(nx+2,ny+2);
v = zeros(nx+2,ny+2);
xir = repmat(E,nx+2);
etar = repmat(N,ny+2);
%% Grouping
AA = 1/((L^2)*(dE^2));
BB = 1/((H^2)*(dN^2));
CC = 1/(2*L*dE);
DD = 1/(2*H*dN);
EE = (1+(cot(thetar))^2)/((dE^2)*(L^2));
FF = 1/(((sin(thetar))^2)*(H^2)*(dN^2));
GG = cos(thetar)/(2*L*H*((sin(thetar))^2)*(dN*dE));
HH = (1/(sin(thetar)))*DD;
II = (cot(thetar))*CC;
%% Relaxation Factors
w1 = 1.85;
w2 = 1.85;
w3 = 1.85;
w4 = 1.85;
%% Initial Explicit Temperature profile
%BC's
%Left and Right
for jj = 2:(nx+1)
tpi(jj,1) = 100;
tpi(jj,ny+2) = 50;
end
for iter1 = 1:50000
for ii = 2:(nx+1)
for jj = 2:(ny+1)
tpi_ii = tpi(ii+1,jj) + tpi(ii-1,jj);
tpi_jj = tpi(ii,jj+1) + tpi(ii,jj-1);
tpi_ij = tpi(ii+1,jj+1) + tpi(ii-1,jj-1) - tpi(ii-1,jj+1) - tpi(ii+1,jj-1);
Num1 = EE*tpi_ii + FF*tpi_jj -GG*tpi_ij;
Den1 = (2*(EE+FF)) - (xir(ii,jj)*L + etar(ii,jj)*H*cos(thetar));
tpi(ii,jj) = tpi(ii,jj) + w1*((Num1/Den1) - tpi(ii,jj));
end
end
end
%% Main Loop
for iter2 = 1:10000
%Stream function
for ii = 2:(nx+1)
for jj = 2:(ny+1)
sp_i = sp(ii+1,jj) - sp(ii-1,jj);
sp_j = sp(ii,jj+1) - sp(ii,jj-1);
sp_ii = sp(ii+1,jj) + sp(ii-1,jj);
sp_jj = sp(ii,jj+1) + sp(ii,jj-1);
sp_ij = sp(ii+1,jj+1) + sp(ii-1,jj-1) - sp(ii-1,jj+1) - sp(ii+1,jj-1);
Num2 = EE*sp_ii + FF*sp_jj - GG*sp_ij + vp(ii,jj);
Den2 = 2*(EE + FF);
sp(ii,jj) =sp(ii,jj) +w2*((Num2/Den2) - sp(ii,jj));
end
end
%Velocity function
for ii = 2:(nx+1)
for jj = 2:(ny+1)
sp_i = sp(ii+1,jj) - sp(ii-1,jj);
sp_j = sp(ii,jj+1) - sp(ii,jj-1);
u(ii,jj) = HH*sp_j - II*sp_i;
v(ii,jj) = -1*(CC)*sp_i;
end
end
%Vorticity boundary conditions
for ii = 1:(nx+2)
vp(1,ii) = -((8*sp(2,ii) - sp(3,ii))/(2*(dy^2))); %top
vp(ny+2,ii) = -((8*sp(ny+1,ii) - sp(ny,ii))/(2*(dy^2))); %Bottom
end
for jj = 2:(ny+1)
vp(jj,1) = -((8*sp(jj,2) - sp(jj,3))/(2*(dy^2))); %left
vp(jj,ny+2) = -((8*sp(jj,ny+1) - sp(jj,ny))/(2*(dy^2))); %right
end
%vorticity profile
for ii = 2:(nx+1)
for jj = 2:(ny+1)
vp_i = vp(ii+1,jj) - vp(ii-1,jj);
vp_j = vp(ii,jj+1) - vp(ii,jj-1);
vp_ii = vp(ii+1,jj) + vp(ii-1,jj);
vp_jj = vp(ii,jj+1) + vp(ii,jj-1);
vp_ij = vp(ii+1,jj+1) + vp(ii-1,jj-1) - vp(ii-1,jj+1) - vp(ii+1,jj-1);
tpi_ii = tpi(ii+1,jj) + tpi(ii-1,jj);
tpi_jj = tpi(ii,jj+1) + tpi(ii,jj-1);
tpi_i = tpi(ii+1,jj) - tpi(ii-1,jj);
tpi_j = tpi(ii,jj+1) - tpi(ii,jj-1);
tpi_ij = tpi(ii+1,jj+1) + tpi(ii-1,jj-1) - tpi(ii-1,jj+1) - tpi(ii+1,jj-1);
sp_i = sp(ii+1,jj) - sp(ii-1,jj);
sp_j = sp(ii,jj+1) - sp(ii,jj-1);
Num3 = (Gr/2)*(HH*tpi_j - II*tpi_i) + (vp_i*CC)*(HH*sp_j - II*sp_i) -(CC*sp_i)*(HH*vp_j - II*vp_i) - EE*vp_ii - FF*vp_jj + GG*vp_ij;
Den3 = -2*(EE+FF);
vp(ii,jj) = vp(ii,jj) + w3*((Num3/Den3) - vp(ii,jj));
end
end
for ii = 2:(nx+1)
for jj = 2:(ny+1)
tpi_ii = tpi(ii+1,jj) + tpi(ii-1,jj);
tpi_jj = tpi(ii,jj+1) + tpi(ii,jj-1);
tpi_i = tpi(ii+1,jj) - tpi(ii-1,jj);
tpi_j = tpi(ii,jj+1) - tpi(ii,jj-1);
tpi_ij = tpi(ii+1,jj+1) + tpi(ii-1,jj-1) - tpi(ii-1,jj+1) - tpi(ii+1,jj-1);
sp_i = sp(ii+1,jj) - sp(ii-1,jj);
sp_j = sp(ii,jj+1) - sp(ii,jj-1);
Num4 = AA*tpi_ii + BB*tpi_jj - GG*tpi_ij -Pr*((CC*tpi_i)*(HH*sp_j - II*sp_i) - (sp_i*CC)*(HH*tpi_j - II*tpi_i));
Den4 = 2*(EE+FF);
tpi(ii,jj) = tpi(ii,jj) + w4*((Num4/Den4) - tpi(ii,jj));
end
end
%Adiabatic conditions
%Top
for jj = 2:(nx+1)
tpi(1,jj) = (-Pr*(u(1,jj)*CC*(tpi(1,jj+1) - tpi(1,jj-1))) + AA*(tpi(1,jj+1) - tpi(1,jj-1)) + 2*BB*tpi(2,jj))/(2*(AA+BB));
end
%Bottom
for jj = 2:(nx+1)
tpi(nx+2,jj) = (-Pr*(u(nx+2,jj)*CC*(tpi(nx+2,jj+1) - tpi(nx+2,jj-1))) + AA*(tpi(nx+2,jj+1) - tpi(nx+2,jj-1)) + 2*BB*tpi(nx+1,jj))/(2*(AA+BB));
end
end
Any help will be greatly appreciated.
