FDM using succsesive overrelaxation method,

Here w1 is coming fine but w2 contour graph is not coming. when I checked it in the variables w2 is taking as NaN except all the boundary conditions. Can anyone help me in plotting w2. Same thing is coming for T2. If you have doubt kindly feel free to ask. Your help will be highly appreciable.

Thanking You in Advance

clc;
clear;
close all;
%% Geometric parameters of domain
L = 1;                                                                     %Length along y-direc
H = 1;                                                                     %Length along x-direc
Nx1 = 100;                                                                 %Number of divisions in 0<=x(1)<=1/2
Nx2 = 100;                                                                 %Number of divisions in 1/2<=x(2)<=1 
Ny = 100;                                                                  %Number of divisions in y-direc
dx1 = 1/(2*Nx1);                                                           %Length of element in x-direc in Region-I
dx2 = 1/(2*Nx2);                                                           %Length of element in x-direc in Region-II
dy = 1/Ny;                                                                 %Length of element in y-direc
h = 0.005;
x1 = 0:dx1:1/2;
x2 = 1/2:dx2:1;
y = 0:dy:1;
%% Parameter
Gr = 10; 
Br = 0.1;
P = -0.1;
phi1 = 0.01;
phi2 = 0.01;
ro1 = 8933;                                                                           %Density of nanoparticle in region-I
beta1 = 1.67;                                                                         %Thermal expansion coefficient of nanoparticle in region-I
ka1 = 401;                                                                            %Thermal conductivity of nanoparticle in region-I
ro2 = 10500;                                                                          %Density of nanoparticle in region-II
beta2 = 1.89;                                                                         %Thermal expansion coefficient of nanoparticle in region-II
ka2 = 429;                                                                            %Thermal conductivity of nanoparticle in region-II
rof1 = 884;                                                                           %Density of base fluid in region-I
betaf1 = 70;                                                                          %Thermal expansion coefficient of base fluid in region-I
kaf1 = 0.145;                                                                         %Thermal conductivity of base fluid in region-I
muf1 =0.486;                                                                          %Viscosity of fluid in region-I
rof2 = 920;                                                                           %Density of base fluid in region-II
betaf2 = 64;                                                                          %Thermal expansion coefficient of base fluid in region-II
kaf2 = 0.121;                                                                         %Thermal conductivity of base fluid in region-II
muf2 = 0.0145;                                                                        %Viscosity of fluid in region-II
la = muf2/muf1;                                                                       %Ratio of viscosity
beta = betaf2/betaf1;                                                                 %Ratio of thermal expansion coefficient
n = rof2/rof1;                                                                        %Ratio of density
ka = kaf2/kaf1;                                                                       %Ratio of thermal conductivity
A1 = 1/((1-phi1)^2.5);
A3 = ((1-phi1)+((phi1*ro1*beta1)/(rof1*betaf1)));
A4 = (ka1+2*kaf1-(2*phi1*(kaf1-ka1)))/(ka1+2*kaf1+(phi1*(kaf1-ka1)));
B1 = 1/((1-phi2)^2.5);
B3 = ((1-phi2)+((phi2*ro2*beta2)/(rof2*betaf2)));
B4 = (ka2+2*kaf2-(2*phi2*(kaf2-ka2)))/(ka2+2*kaf2+(phi2*(kaf2-ka2)));
%% Boundary and initial conditions
T1 = zeros(Ny+1,Nx1+1);
w1 = zeros(Ny+1,Nx1+1);
T2 = zeros(Ny+1,Nx2+1);
w2 = zeros(Ny+1,Nx2+1);
w1(:,1) = 0;                                                               %Velocity at lower adiabatic wall(Region-I)
T1(:,1) = T1(:,2);                                                         %Temperature at lower adiabtic wall (Region-I)
w1(1,:) = 0;                                                               %Velocity at left wall in region-I
T1(1,:) = -1;                                                              %Temperature at left wall in region-I
w1(Ny+1,:) = 0;                                                            %Velocity at right wall in region-I
T1(Ny+1,:) = 1;                                                            %Temperature at right wall in region-I
w1(:,Nx1+1) = w2(:,1);                                                     %Continuity of velocity at interface
T1(:,Nx1+1) = T2(:,1);                                                     %Continuity of Temperature at interface
w1(:,Nx1+1) = (((la*B1*dx1)/(dx2*A1))*(w2(:,2)-w2(:,1)))+w1(:,Nx1);        %Continuity of shear stress
T1(:,Nx1+1) = (((ka*B4*dx1)/(dx2*A4))*(T2(:,2)-T2(:,1)))+T1(:,Nx1);        %Continuity of heat flux
w2(1,:) = 0;                                                               %Velocity at left wall in region-II
T2(1,:) = -1;                                                              %Temperature at left wall in region-II
w2(Ny+1,:) = 0;                                                            %Velocity at right wall in region-II
T2(Ny+1,:) = 1;                                                            %Temperature at right wall in region-II
w2(:,Nx2+1) = 0;                                                           %Velocity at upper adiabatic wall(Region-II)
T2(:,Nx2+1) = T2(:,Nx2);                                                   %Temperature at upper adiabtic wall (Region-II)
%% Convergence
Epsilon = 1.e-8;
Error = 2*Epsilon;                                                                            %Error is any value greater than epsilon
omega1 = 1.8;                                                                 %Relaxation parameter
omega2 = 2*0.9;                                                                 %Relaxation parameter of temperature in region-I
omega3 = 1.8;                                                                 %Relaxation parameter of velocity in region-I
omega4 = 2*0.9;                                                                 %Relaxation parameter of temperature in region-II
%% SOR Loop
while (Error > Epsilon)
    w1_old = w1; T1_old = T1; w2_old = w2; T2_old = T2; 
    for j = 2:Ny
        for i = 2:Nx1
            w1(i,j) = (w1(i,j)*(1-omega1))+((2*omega1)/5)*(w1(i+1,j)+w1(i-1,j)+(w1(i,j+1)/4)+(w1(i,j-1)/4)+((Gr*A3*h^2*T1(i,j))/A1)-(P*h^2/A1)); 
            T1(i,j) = (T1(i,j)*(1-omega2))+((2*omega2)/5)*(T1(i+1,j)+T1(i-1,j)+(T1(i,j+1)/4)+(T1(i,j-1)/4)+((Br*A1/A4)*(((w1(i+1,j)-w1(i-1,j))^2)/4)+(((w1(i,j+1)-w1(i,j-1))^2)/16)));
            %f(i,j) =  (((w1(i+1,j)-w1(i-1,j))^2)/4)+(((w1(i,j+1)-w1(i,j-1))^2)/16);
        end
    end
    Error_w1=max(abs(w1(:)-w1_old(:)));
    Error_T1=max(abs(T1(:)-T1_old(:)));
    for j = 2:Ny
        for i = 2:Nx2
            w2(i,j) = (w2(i,j)*(1-omega3))+((2*omega3)/5)*(w2(i+1,j)+w2(i-1,j)+(w2(i,j+1)/4)+(w2(i,j-1)/4)+((Gr*B3*h^2*beta*n*T2(i,j))/la*B1)-(P*h^2/la*B1)); 
            T2(i,j) = (T2(i,j)*(1-omega4))+((2*omega4)/5)*(T2(i+1,j)+T2(i-1,j)+(T2(i,j+1)/4)+(T2(i,j-1)/4)+((Br*la*B1/ka*B4)*(((w2(i+1,j)-w2(i-1,j))^2)/4)+(((w2(i,j+1)-w2(i,j-1))^2)/16)));
            %g(i,j) =  (((w2(i+1,j)-w2(i-1,j))^2)/4)+(((w2(i,j+1)-w2(i,j-1))^2)/16);
        end 
    end 
    Error_w2=max(abs(w2(:)-w2_old(:)));
    Error_T2=max(abs(T2(:)-T2_old(:)));
    Error=max([Error_w1, Error_T1, Error_w2, Error_T2]);
end 
%% Plotting the results
% Plotting the results for Region-1 (w1, T1)
figure;
contourf(x1, y, w1, 'LineStyle', 'none')
colorbar

figure;
contourf(x2, y, w2, 'LineStyle', 'none')

Warning: Contour not rendered for constant ZData

colorbar;
title('Velocity distribution');
xlabel('X-direction');
ylabel('Y-direction');
axis equal;

[Torsten]

Your resulting system of equations in the unknowns w and T is nonlinear. Why do you use a method to iteratively solve linear systems of equations ? Use "fsolve" instead.

Further look at what you do under the part %% Boundary and initial conditions. You set relations at the boundaries and at the interface for matrices that are all zero at the beginning. But the conditions have to be updated permanently in the while loop because your conditions at the boundaries and at the interface change with iteration.

[Komal Goyal]

@Torsten Actually I havr ro solve this using iterative technique, futher in order to compare the accuracy of my results I will use "fsolve". By comparing the outcomes by both the technique.

Yes sir my concern is also the same, since I set all the boundary and interface coniditions zero initially, and the value should be updated in the while loop, in the case of region-I, values are updating but I dont know what mistake I did in region-II, in which except boundary and interface it taking value as NaN

[Torsten]

As said, you don't update boundary and interface conditions in the while loop.

And there is no guarantee that SOR will always converge.

[Komal Goyal]

@Torsten Yes, you are correct, I know i have to update boundary and interface conditions in the while loop but I am not able to do that. Can you please help me by doing some corrections in the code so that so I learn how can we do that.

Thanking you @Torsten for your suggestions always;)

Answer 1

[Alan Stevens]

For w1 and T1 you have

w1(:,Nx1+1) = (((la*B1*dx1)/(dx2*A1))*(w2(:,2)-w2(:,1)))+w1(:,Nx1);        %Continuity of shear stress
T1(:,Nx1+1) = (((ka*B4*dx1)/(dx2*A4))*(T2(:,2)-T2(:,1)))+T1(:,Nx1);        %Continuity of heat flux

There are no equivalent assignments for w2 and T2. Should there be?

[Komal Goyal]

@Alan Stevens This condition is for continuity at the first order derivative of velocity in region-1 equal to (((la*B1*dx1)/(dx2*A1)) time first order derivative of velocity in region 2, i.e continuity of shear stress at the interface of two regions, You can see the in the previous code comment section https://in.mathworks.com/matlabcentral/answers/2056659-finite-difference-method-for-two-layered-model?s_tid=srchtitle sir, that I have 4 equations of second order in two independent variables, and for them I have 16 boundary and interface conditions so it should plot the graph. SInce in the region-I contour plot is coming exactly fine, what is the mistake I am doing in region-II. While having discussion with @Torsten in the previous code I have two attachments. Kindly refer to them for any doubt.

Thanking you for your sugesstions.