FDM using succsesive overrelaxation method,

3 次查看（过去 30 天）
Komal Goyal
Komal Goyal 2023-12-15
评论： Komal Goyal 2023-12-18
Ran in:
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;
  4 个评论
显示 2更早的评论
Torsten
Torsten 2023-12-18
编辑：Torsten 2023-12-18
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
Komal Goyal 2023-12-18
@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;)

请先登录，再进行评论。

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

回答（1 个）

Alan Stevens
Alan Stevens 2023-12-17
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?
  1 个评论
Komal Goyal
Komal Goyal 2023-12-17
编辑：Komal Goyal 2023-12-17
@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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fluid Mechanics 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by