Finite-difference method for two layered model

15 次查看（过去 30 天）
Komal Goyal
Komal Goyal 2023-12-5
评论： Komal Goyal 2023-12-7
Ran in:
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.
  7 个评论
显示 5更早的评论
Komal Goyal
Komal Goyal 2023-12-7
@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
Komal Goyal 2023-12-7
@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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Jason Shrand
Jason Shrand 2023-12-6
编辑：Jason Shrand 2023-12-6
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
  2 个评论
Jason Shrand
Jason Shrand 2023-12-6
Ran in:
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
Komal Goyal 2023-12-6
@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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Timing and presenting 2D and 3D stimuli 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by