FDM using succsesive overrelaxation method,
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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')
colorbar;
title('Velocity distribution');
xlabel('X-direction');
ylabel('Y-direction');
axis equal;
4 commentaires
Réponses (1)
Alan Stevens
le 17 Déc 2023
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 commentaire
Voir également
Catégories
En savoir plus sur Spline Postprocessing dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!