How to implement a mixed boundary conditions into 2D steady state heat conduction equation?
12 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hi
I have 2D steady heat conduction equation on the unit square subject to the following mixed Dirichlet/Neumann boundary conditions.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/364405/image.png)
T(x,1)=sin(
x) T(1,y)=1-y
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/364408/image.png)
Use uniform grid in x and y of 21 points in each direction. Apply an initial condition of T(x,y)=0 . Iterate until the maximum change is less than 0.1. And the tolerance value 1.0e-4
here is my code but I'm having difficulties implement the Neumann boundary condition.
for I=1:2
tic
nx = 21; %step size x -direction
ny = 21; %step size y -direction
Lx = 1;
Ly = 1;
x = linspace(0,Lx,nx);
y = linspace(0,Ly,nx);
dx = x(2)-x(1);
dy = y(2)-y(1);
k = 0;
% Intial condition
T = zeros(nx,ny);
% BC
T(1:nx,nx) = (sin(pi*x(:,1)));
T(1:nx,1) = 5;
T(ny,1:ny) = 0;
T(1,1:ny) = (1-y(1,:));
Told = T;
tol = 1.0e-4;
error = 1;
k = k+1;
% iterative solver
iterative_solver = input('iterative_solver(1. jacobi 2. gauss) = ');
% Jacobi
if iterative_solver ==1
jacobi_iter = 1;
while(error>tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = 0.25*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1));
end
end
error = max(max(abs(Told-T)));
jacobi_iter = jacobi_iter+1;
Told=T;
end
contourf(x,y,T,'ShowText','on')
colorbar;
title_text=sprintf('total jacobi iteration = %d @steady state',jacobi_iter);
title(title_text)
xlabel('x-axis');
ylabel('y-axis');
end
% Gauss
if iterative_solver ==2
gs_iter =1;
while (error>tol)
for i = 2:nx-1
for j = 2:ny-1
T(i,j) = (0.25)*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
end
end
error = max(max(abs(Told-T)));
gs_iter=gs_iter+1;
Told = T;
end
contourf(x,y,T,'ShowText','on')
colorbar;
xlabel('x-axis')
ylabel('y-axis')
title_text=sprintf('total gauss seidel iteration @steady state = %d',gs_iter);
title(title_text)
end
if I==1
ct_j =toc;
fprintf('%d jacobi iterations for computational time of %0.3g secondsn',k,ct_j)
elseif I==2
ct_gs = toc;
fprintf('%d gauss seidel iterations for computational time of %0.3g secondsn',k,ct_gs)
end
end
0 Kommentare
Antworten (1)
Alan Stevens
am 23 Sep. 2020
Like this? I've used: (T(1:nx,2) - T(1:nx,1))/dy = 5 and rearranged this to get T(1:nx,1). In addition I've set the y values to decrease from ny-1, instead of increasing from 2 and shifted the calculation of T(1:nx,1) to the end of the loop. Check carefully!!
for I=1:2
tic
nx = 21; %step size x -direction
ny = 21; %step size y -direction
Lx = 1;
Ly = 1;
x = linspace(0,Lx,nx);
y = linspace(0,Ly,nx);
dx = x(2)-x(1);
dy = y(2)-y(1);
k = 0;
% Intial condition
T = zeros(nx,ny);
% BC
T(1:nx,nx) = (sin(pi*x(:,1)));
%T(1:nx,1) = 5;
T(ny,1:ny) = 0;
T(1,1:ny) = (1-y(1,:));
Told = T;
tol = 1.0e-4;
error = 1;
k = k+1;
% iterative solver
iterative_solver = input('iterative_solver(1. jacobi 2. gauss) = ');
% Jacobi
if iterative_solver ==1
jacobi_iter = 1;
while(error>tol)
for i = 2:nx-1
for j = ny-1:-1:2 %%%%%%%%%%%%%%%%
T(i,j) = 0.25*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1));
end
end
T(1:nx,1) = T(1:nx,2) - 5*dy; %%%%%%%%%%%%%
error = max(max(abs(Told-T)));
jacobi_iter = jacobi_iter+1;
Told=T;
end
contourf(x,y,T,'ShowText','on')
colorbar;
title_text=sprintf('total jacobi iteration = %d @steady state',jacobi_iter);
title(title_text)
xlabel('x-axis');
ylabel('y-axis');
end
% Gauss
if iterative_solver ==2
gs_iter =1;
while (error>tol)
for i = 2:nx-1
for j = ny-1:-1:2 %%%%%%%%%%%%%%%%%%%
T(i,j) = (0.25)*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
end
end
T(1:nx,1) = T(1:nx,2) - 5*dy; %%%%%%%%%%%%%%%%%%%
error = max(max(abs(Told-T)));
gs_iter=gs_iter+1;
Told = T;
end
contourf(x,y,T,'ShowText','on')
colorbar;
xlabel('x-axis')
ylabel('y-axis')
title_text=sprintf('total gauss seidel iteration @steady state = %d',gs_iter);
title(title_text)
end
if I==1
ct_j =toc;
fprintf('%d jacobi iterations for computational time of %0.3g secondsn',k,ct_j)
elseif I==2
ct_gs = toc;
fprintf('%d gauss seidel iterations for computational time of %0.3g secondsn',k,ct_gs)
end
end
0 Kommentare
Siehe auch
Kategorien
Mehr zu Calculus finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!