Finite-difference method for two layered model

I am a beginner in MATLAB, and I am attempting to code a two-layered system using the finite-difference method. In my code, region I is defined for 0 <= x1 < 1/2, and region II for 1/2 < x2 <= 1, where 0 <= y <= 1. My goal is to create a contour plot for the temperature distribution in a single figure window covering the entire domain (0 <= x <= 1 and 0 <= y <= 1). I want region I to display the temperature T1 and region II to display the temperature T2. However, my code contains errors that I'm struggling to identify and correct. I would greatly appreciate your assistance in resolving these issues. Please feel free to seek clarification if needed. Thank you in advance for your help.

clc;
clear all;


%% Geometric parameters for the domain
L = 1;                                                                                %Length of the duct in y-direc(in m)
H = 1;                                                                                %Width of the duct in x-direc(in m)
Nx1 = 101;                                                                            %Number of grid point in region-I in x-direc
Nx2 = 101;                                                                            %Number of grid point in region-II in x-direc
Ny = 101;                                                                             %Number of grid point in y-direc
dx1 = H/(Nx1-1);                                                                      %Length of the element in x-direc in region-I(in m)
dx2 = H/(Nx2-1);                                                                      %Lenght of the element in x-direc in region-II(in m)
dy = L/(Ny-1);                                                                        %Lenght of the element in y-direc(in m)

%% 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)+((phi1*ro1*beta1)/(rof1*betaf1)));
B4 = (ka1+2*kaf1-(2*phi1*(kaf1-ka1)))/(ka1+2*kaf1+(phi1*(kaf1-ka1)));

%% Boundary and initial conditions 
T1 = ones(Nx1,Ny);                                                                    %Creating an array for temperature in region-I of order Nx1*Ny
T2 = ones(Nx2,Ny);                                                                    %Creating an array for temperature in region-I of order Nx1*Ny 
w1 = ones(Nx1,Ny);                                                                    %Creating an array for velocity in region-I of order Nx1*Ny
w2 = ones(Nx2,Ny);                                                                    %Creating an array for velocity in region-II of order Nx2*Ny
w1(:,1) = -w1(:,1);                                                                   %Velocity conditon for left wall region-I(y=0, 0<=x1<A1)
T1(:,1) = -1-T1(:,1);                                                                 %Temperature conditon for left wall region-I(y=0, 0<=x1<A1)
w1(:,Ny) = -w1(:,Ny);                                                               %Velocity conditon for right wall region-I(y=1, 0<=x1<A1)
T1(:,Ny) = 1-T1(:,Ny);                                                              %Temperature conditon for right wall region-I(y=1, 0<=x1<A1)
w1(1,:) = -w1(1,:);                                                                   %Velocity condition at lower adiabtic wall(x=0,0<=y<=1)
T1(1,:) = T1(1,:);                                                                    %Temperature condition at lower adiabtic wall(x=0,0<=y<=1) 
w2(Nx2,:) = w1(Nx1,:)+w1(Nx1,:)-w2(Nx2,:);                                        %Continuity of velocity at interface(x=A1, 0<=y<=1)
w1(Nx1,:) = (((la*B1*dx1)/(dx2*A1))*(w2(Nx2,:)-w2(Nx2,:)))+w1(Nx1,:);             %Continuity of tangential velocity/shear stress at the interface(x=A1, 0<=y<=1)
T2(Nx2,:) = T1(Nx1,:)+T1(Nx1,:)-T2(Nx2,:);                                        %Continuity of Temperature at interface(x=A1, 0<=y<=1)
T1(Nx1,:) = (((ka*B4*dx1)/(dx2*A4))*(T2(Nx2,:)-T2(Nx2,:)))+T1(Nx1,:);             %Continuity of heat flux at the interface(x=A1, 0<=y<=1)
w2(:,1) = -w2(:,1);                                                                   %Velocity conditon for left wall region-II(y=0, A1<x2<=A2)
T2(:,1) = -1-T2(:,1);                                                                 %Temperature conditon for left wall region-II(y=0, A1<x2<=A2)
w2(:,Ny) = -w2(:,Ny);                                                               %Velocity conditon for right wall region-II(y=1, A1<x2<=A2)
T2(:,Ny) = 1-T2(:,Ny);                                                              %Temperature conditon for right wall region-II(y=1, A1<x2<=A2)
w2(1,:) = -w2(1,:);                                                                   %Velocity condition at upper adiabtic wall(x=A1+A2,0<=y<=1)
T2(1,:) = T2(1,:);                                                                    %Temperature condition at upper adiabtic wall(x=A1+A2,0<=y<=1) 

%% Convergence
Epsilon = 1e-8;
Error = 5;                                                                            %Error is any value greater than epsilon
  
%% Plotting the results
% Plotting the results for Region-1 (w1, T1)
x1 = linspace(0,0.5,Nx1);
x2 = linspace(0.5,1,Nx2);
y = 0:dy:L;
figure;
contourf(x1, y, T1, 'LineStyle', 'none');
hold on
contourf(x2, y, T2, 'LineStyle', 'none');
hold off
colorbar;
title('Temperature distribution in Region-1');
xlabel('X-direction');
ylabel('Y-direction');
axis equal;
%% Computation
Iter = 0; 

while(Error>Epsilon)
      
    Iter = Iter + 1;
    disp(Iter);

    for j = 2:Ny-1 
    for i = 2:Nx1-1
    % Equations for region-I    
eq1 = (w1(i+1,j)-2*w1(i,j)+w1(i-1,j)/(dx1^2))+(w1(i,j+1)-2*w1(i,j)+w1(i,j-1)/(dy^2))+((Gr*A3*T1(i,j))/A1)-(P/A1);
eq2 = (T1(i+1,j)-2*T1(i,j)+T1(i-1,j)/(dx1^2))+(T1(i,j+1)-2*T1(i,j)+T1(i,j-1)/(dy^2))+((Br*A1/A4)*(((w1(i+1,j)-w1(i-1,j))/(2*dx1))^2+((w1(i,j+1)-w1(i,j-1))/(2*dy))^2));
     end
    for j = 2:Ny-1
     for i = 2:Nx2-1
         % Equations for region-II
       eq3= (w2(i+1,j)-2*w2(i,j)+w2(i-1,j)/(dx2^2))+(w2(i,j+1)-2*w2(i,j)+w2(i,j-1)/(dy^2))+((Gr*B3*beta*n*T2(i,j))/B1)-(P/(B1*la));
eq4 = (T2(i+1,j)-2*T2(i,j)+T2(i-1,j)/(dx2^2))+(T2(i,j+1)-2*T2(i,j)+T2(i,j-1)/(dy^2))+((Br*la*B1/ka*B4)*(((w2(i+1,j)-w2(i-1,j))/(2*dx2))^2+((w2(i,j+1)-w2(i,j-1))/(2*dy))^2));
    end
    end
    end

    %Error = sqrt(sumsqr(T1 - T1_old) + sumsqr(T2 - T2_old)); % Use T1_old and T2_old to store the previous values

    % Update T1_old and T2_old for the next iteration
    %T1_old = T1;
    %T2_old = T2;
    %disp(Error);
end

Index in position 2 is invalid. Array indices must be positive integers or logical values.

[Torsten]

What is the equation you try to solve ?

What are the boundary conditions at x = 0, x = 1, y = 0 and y = 1 ?

What is the transmission condition at x = 1/2 ?

[Komal Goyal]

@Torsten I have shared one image file with all the governing equations and boundary conditions. Kindly look into that.

[Torsten]

Do you have the problem description in a non-discretized form ?

[Komal Goyal]

Yes, @Torsten, sir. That also, I have attached. Kindly see that once also.

[Torsten]

At a first glance, your equations and boundary conditions look reasonable.

But from your code, I get the impression that your experience with MATLAB and numerical methods will not suffice to solve this 2d-problem if you try to program the numerical method on your own. You will have to implement a nonlinear solver because the discrete expressions for (dw1/dx1)^2, (dw1/dy)^2, (dw2/dx2)^2 and (dw2/dy)^2 will make your discretized system of equations nonlinear - and this for a 2d-problem with a layer interface.

Maybe the PDE toolbox can help, but I'm not sure whether it supports your kind of equations and the general transmission conditions you want to impose.

[Komal Goyal]

@Torsten, thankyou for answering. Actually, you are correct sir, I am the beginner in the field. I am trying to learn and program the numerical method at my own. Can you please help me in that. I think there is some problem in the indexing, and some more minor errors if we will try to rectify that I think the problem will be solved

[Komal Goyal]

@Torsten I have tried as per you said,sir. Now, the code is working fine. But, I am thinking some error are there specially at interface. Can you please see that once.

Answer 1

[Jason Shrand]

The first issue I see is on line 17, in the line of code here:

w1(:,0) = -w1(:,1);   

. In MATLAB, array indices start at 1, not 0. But here, and in several other locations, you're trying to set index 0 of an array

[Jason Shrand]

I found a few issues and addressed them in the updated code below. You'll still need to compute the "Error" variable in your while loop, right now it uses an undefined variable T. I assume you want either T1 or T2?

clc;
clear all;
%% Geometric parameters for the domain
L = 1;                                                                                %Length of the duct in y-direc(in m)
H = 1;                                                                                %Width of the duct in x-direc(in m)
Nx1 = 101;                                                                            %Number of grid point in region-I in x-direc
Nx2 = 101;                                                                            %Number of grid point in region-II in x-direc
Ny = 101;                                                                             %Number of grid point in y-direc
dx1 = H/(Nx1-1);                                                                      %Length of the element in x-direc in region-I(in m)
dx2 = H/(Nx2-1);                                                                      %Lenght of the element in x-direc in region-II(in m)
dy = L/(Ny-1);                                                                        %Lenght of the element in y-direc(in m)
%% 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)+((phi1*ro1*beta1)/(rof1*betaf1)));
B4 = (ka1+2*kaf1-(2*phi1*(kaf1-ka1)))/(ka1+2*kaf1+(phi1*(kaf1-ka1)));

%% Boundary and initial conditions 
T1 = ones(Nx1,Ny);                                                                    %Creating an array for temperature in region-I of order Nx1*Ny
T2 = ones(Nx2,Ny);                                                                    %Creating an array for temperature in region-I of order Nx1*Ny 
w1 = ones(Nx1,Ny);                                                                    %Creating an array for velocity in region-I of order Nx1*Ny
w2 = ones(Nx2,Ny);                                                                    %Creating an array for velocity in region-II of order Nx2*Ny
w1(:,1) = -w1(:,1);                                                                   %Velocity conditon for left wall region-I(y=0, 0<=x1<A1)
T1(:,1) = -1-T1(:,1);                                                                 %Temperature conditon for left wall region-I(y=0, 0<=x1<A1)
w1(:,Ny) = -w1(:,Ny);                                                               %Velocity conditon for right wall region-I(y=1, 0<=x1<A1)
T1(:,Ny) = 1-T1(:,Ny);                                                              %Temperature conditon for right wall region-I(y=1, 0<=x1<A1)
w1(1,:) = -w1(1,:);                                                                   %Velocity condition at lower adiabtic wall(x=0,0<=y<=1)
T1(1,:) = T1(1,:);                                                                    %Temperature condition at lower adiabtic wall(x=0,0<=y<=1) 
w2(Nx2,:) = w1(Nx1,:)+w1(Nx1,:)-w2(Nx2,:);                                        %Continuity of velocity at interface(x=A1, 0<=y<=1)
w1(Nx1,:) = (((la*B1*dx1)/(dx2*A1))*(w2(Nx2,:)-w2(Nx2,:)))+w1(Nx1,:);             %Continuity of tangential velocity/shear stress at the interface(x=A1, 0<=y<=1)
T2(Nx2,:) = T1(Nx1,:)+T1(Nx1,:)-T2(Nx2,:);                                        %Continuity of Temperature at interface(x=A1, 0<=y<=1)
T1(Nx1,:) = (((ka*B4*dx1)/(dx2*A4))*(T2(Nx2,:)-T2(Nx2,:)))+T1(Nx1,:);             %Continuity of heat flux at the interface(x=A1, 0<=y<=1)
w2(:,1) = -w2(:,1);                                                                   %Velocity conditon for left wall region-II(y=0, A1<x2<=A2)
T2(:,1) = -1-T2(:,1);                                                                 %Temperature conditon for left wall region-II(y=0, A1<x2<=A2)
w2(:,Ny) = -w2(:,Ny);                                                               %Velocity conditon for right wall region-II(y=1, A1<x2<=A2)
T2(:,Ny) = 1-T2(:,Ny);                                                              %Temperature conditon for right wall region-II(y=1, A1<x2<=A2)
w2(1,:) = -w2(1,:);                                                                   %Velocity condition at upper adiabtic wall(x=A1+A2,0<=y<=1)
T2(1,:) = T2(1,:);                                                                    %Temperature condition at upper adiabtic wall(x=A1+A2,0<=y<=1) 
%% Convergence
Epsilon = 1e-8;
Error = 5;                                                                            %Error is any value greater than epsilon

%% Governing equations
%(w1(i+1,j)-2*w1(i,j)+w1(i-1,j)/(dx1^2))+(w1(i,j+1)-2*w1(i,j)+w1(i,j-1)/(dy^2))+((Gr*A3*T1(i,j))/A1)-(P/A1)==0;
%(T1(i+1,j)-2*T1(i,j)+T1(i-1,j)/(dx1^2))+(T1(i,j+1)-2*T1(i,j)+T1(i,j-1)/(dy^2))+((Br*A1/A4)*(((w1(i+1,j)-w1(i-1,j))/(2*dx1))^2+((w1(i,j+1)-w1(i,j-1))/(2*dy))^2))==0;
%(w2(i+1,j)-2*w2(i,j)+w2(i-1,j)/(dx2^2))+(w2(i,j+1)-2*w2(i,j)+w2(i,j-1)/(dy^2))+((Gr*B3*beta*n*T2(i,j))/B1)-(P/(B1*la))==0;
%(T2(i+1,j)-2*T2(i,j)+T2(i-1,j)/(dx2^2))+(T2(i,j+1)-2*T2(i,j)+T2(i,j-1)/(dy^2))+((Br*la*B1/ka*B4)*(((w2(i+1,j)-w2(i-1,j))/(2*dx2))^2+((w2(i,j+1)-w2(i,j-1))/(2*dy))^2))==0;
%% Plotting the results
x1 = linspace(0, 0.5, Nx1); %x1 = 0:dx1:1/2;
x2 = linspace(0.5, 1, Nx2); % x2 = 1/2:dx2:H;
y = 0:dy:L;
colormap(jet);
contourf(x1,y,T1);

contourf(x2,y,T2);
colorbar;
title('Temperature distribution');
xlabel('X-direction');
ylabel('Y-direction');

%% Computation
Iter = 0; 
while(Error>Epsilon)
      
    Iter = Iter + 1;
    disp(Iter);
    % TODO: Iterating from 2 to N-1 is a temporary work-around just to get something up and
    % running. You'll need to update your logic to handle the boundary
    % conditions (i=1, j=1, i=N, j=N) appropriately.
    for j = 2:(Ny-1)
        for i = 2:(Nx1-1)
            eq1 = (w1(i+1,j)-2*w1(i,j)+w1(i-1,j)/(dx1^2))+(w1(i,j+1)-2*w1(i,j)+w1(i,j-1)/(dy^2))+((Gr*A3*T1(i,j))/A1)-(P/A1);
            eq2 = (T1(i+1,j)-2*T1(i,j)+T1(i-1,j)/(dx1^2))+(T1(i,j+1)-2*T1(i,j)+T1(i,j-1)/(dy^2))+((Br*A1/A4)*(((w1(i+1,j)-w1(i-1,j))/(2*dx1))^2+((w1(i,j+1)-w1(i,j-1))/(2*dy))^2));
            end
         for i = 2:(Nx2-1)
           eq3 = (w2(i+1,j)-2*w2(i,j)+w2(i-1,j)/(dx2^2))+(w2(i,j+1)-2*w2(i,j)+w2(i,j-1)/(dy^2))+((Gr*B3*beta*n*T2(i,j))/B1)-(P/(B1*la));
           eq4 = (T2(i+1,j)-2*T2(i,j)+T2(i-1,j)/(dx2^2))+(T2(i,j+1)-2*T2(i,j)+T2(i,j-1)/(dy^2))+((Br*la*B1/ka*B4)*(((w2(i+1,j)-w2(i-1,j))/(2*dx2))^2+((w2(i,j+1)-w2(i,j-1))/(2*dy))^2));
        end
    end
    Error = sqrt(sumsqr(T(i,:) - T(i-1,:)));
    disp(Error);
end

     1

Unrecognized function or variable 'T'.

[Komal Goyal]

@Jason Shrand Thank you for your answer. Can you tell me If we want to plot both temperature T1 and T2 in the same figure T1 is from 0<=x1<1/2 and T2 is from 1/2<x2<=1. So, how can we plot. Moreover, I have attached one image with the governing equation and boundary conditions in reply to @Torsten. Can you please check whether I wrote all the equations correctly or not. Though, I checked but I have some doubt. Thank you for your help once again sir