Added support for boundary conditions of the 2nd type (von Neumann) -…

…- heat fluxes at the walls

heat3D_trans.m

@@ -31,12 +31,12 @@
 end
 
 %% domain
-Lx= 1;                  % plate width (m)
-Ly= 1;                  % plate depth (m)
-Lz= 1;                  % plate height (m)
-Nx=51;                 % nodes in x direction (odd number to have central plane to slice through)
-Ny=51;                 % nodes in y direction (odd number to have central plane to slice through)
-Nz=51;                  % nodes in z direction (odd number to have central plane to slice through)
+Lx= 2;                  %  width (m)
+Ly= 1;                  %  depth (m)
+Lz= 1;                  %  height (m)
+Nx=41;                 % nodes in x direction (odd number to have central plane to slice through)
+Ny=21;                 % nodes in y direction (odd number to have central plane to slice through)
+Nz=21;                  % nodes in z direction (odd number to have central plane to slice through)
 dx=Lx/Nx;               % spacing along x
 dy=Ly/Ny;               % spacing along y
 dz=Lz/Nz;               % spacing along z
@@ -45,77 +45,124 @@
 yc=(Ny+1)/2;
 zc=(Nz+1)/2;
 
-dt = 1/(2*alph*((1/dx^2)+(1/dy^2)+(1/dz^2)));	% stable time step, sec
+dt = 1/(2*alph*((1/dx^2)+(1/dy^2)+(1/dz^2)));	% stable time step, sec for Euler ethod
 tmax=3600; % max time, sec
 
 %% Initial and boundary conditions (1st type, a.k.a Dirichlet)
+% T(y,x,z) % structure of 3D array
+
 T_0= 0;                 % Initial temperature in all nodes 
-T_right   = -20;        % temperature at x=Lx
+T_right  = -20;        % temperature at x=Lx
 T_left   = 20;        % temperature at x=0 
-T_top  = 0;             % temperature at y=0 
-T_bott  = 0;            % temperature at y=Ly  
+
+T_top  = 0;             % temperature at z=Lz 
+T_bott  = 0;            % temperature at z=0  
+
+T_front=0;            % temperature at y=0 
+T_back=0;            % temperature at y=Ly
 
 tol_ss=0.01;            % tolerance for steady state
 
+%% Initial conditions for Steady State
+T=zeros(Ny,Nx,Nz);
+
+T(:,1,:)    =  T_left;
+T(:,Nx,:)   =  T_right;
+
+T(1,:,:)    =  T_front;
+T(Ny,:,:)   =  T_back;
+
+T(:,:,1)    =  T_bott;
+T(:,:,Nz)   =  T_top;
+
+%% Heat fluxes at top and bottom (2nd type -- von Neuman)
+
+% enable calculation of temperature
+% from heat fluxes at the walls (2nd type B.C. at the walls)
+BC2=1; % set to 0 to disable
+
+if BC2==1 % see also lines 130++
+    Q_top   =  0; % Q=0 means thermal insulation
+    Q_bott  =  0;
+
+    Q_left  =  0;
+    Q_right =  0;
+
+    Q_front =  0;
+    Q_back  =  0;
+end
 %% Constants
 
 iter=1;  % counter of iterations
 error  = 1; % initial error for iterations
 
-%% Initial conditions for Steady State
-Tss=zeros(Ny,Nx,Nz);
+tmax=3600; % max time, sec 
+t=0;
+k1=alph*dt; 
 
-Tss(:,1,:)    =T_left;
-Tss(:,Nx,:)   =T_right;
+%% plotting 
+T_min=T_right;
+T_max=T_left;
 
-Tss(1,:,:)    =T_bott;
-Tss(Ny,:,:)   =T_top;
+numisosurf=25; % number of isosurfaces
+isovalues=linspace(T_min,T_max,numisosurf);
 
-%% plotting 
 x=linspace(0,Lx,Nx);
 y=linspace(0,Ly,Ny);
 z=linspace(0,Lz,Nz);
 [X Y Z]=meshgrid(x,y,z);
 
 figure('units','pixels','position',[100 100 1280 720])
-
 
 %% Steady-State Solution
 tic
 i=2:Ny-1; %inner nodes along y
 j=2:Nx-1; %inner nodes along x
-k=2:Nx-1; %inner nodes along z
-Tss_new=Tss;
-T_min=T_right;
-T_max=T_left;
-numisosurf=25; % number of isosurfaces
-isovalues=linspace(T_min,T_max,numisosurf);
-tmax=3600; % max time, sec 
-t=0;
-k1=alph*dt; 
+k=2:Nz-1; %inner nodes along z
+T_new=T;
+
 method='steadystate';
 
 while t<=tmax%error >= tol_ss % by tmax or by tolerance (for Staeady State)
+
+%% von Neuman B.C. (fluxes at walls)
+% grad T as central difference    
+if BC2==1
+
+% grad T via central difference -- Neuman B.C. at the walls (constant fluxes or grad(T) )
+% uncomment the respective walls
+
+        T(1,:,:)  =  T(3,:,:)     +Q_front*2*dy;     %  flux on front
+        T(Ny,:,:) =  T(Ny-2,:,:)  +Q_back*2*dy;      %  flux on back
+
+%       T(:,1,:)  =  T(:,3,:)    +Q_left*2*dx;       %  flux on left
+%       T(:,Nx,:) =  T(:,Nx-2,:) +Q_right*2*dx;      %  flux on right
+
+%       T(:,:,1)  =  T(:,:,3)    +Q_bott*2*dz;       %  flux on bottom
+%       T(:,:,Nz) =  T(:,:,Nz-2) +Q_top*2*dz;        %  flux on top
+end    
 
+
+%%     
   switch method
         case 'steadystate'
 
    % Solving Steady-State Laplace's eqn on 3D FD stencil
-    Tss_new(i,j,k)=(dz^2*(Tss(i,j,k-1)+Tss(i,j,k+1))+dy^2*(Tss(i+1,j,k)+Tss(i-1,j,k))+dx^2*(Tss(i,j+1,k)+Tss(i,j-1,k)))/(2*(dx^2+dy^2+dz^2));    
+    T_new(i,j,k)=(dz^2*(T(i,j,k-1)+T(i,j,k+1))+dy^2*(T(i+1,j,k)+T(i-1,j,k))+dx^2*(T(i,j+1,k)+T(i,j-1,k)))/(2*(dx^2+dy^2+dz^2));    
 
         case 'Euler'
 
     % FTCS transient -- forward Euler
-   Tss_new(i,j,k) =Tss(i,j,k)+k1*(((Tss(i-1,j,k)-2*Tss(i,j,k)+Tss(i+1,j,k))/dy^2)+((Tss(i,j-1,k)-2*Tss(i,j,k)+Tss(i,j+1,k))/dx^2)+((Tss(i,j,k-1)-2*Tss(i,j,k)+Tss(i,j,k+1))/dz^2));
-   
+   T_new(i,j,k) =T(i,j,k)+k1*(((T(i-1,j,k)-2*T(i,j,k)+T(i+1,j,k))/dy^2)+((T(i,j-1,k)-2*T(i,j,k)+T(i,j+1,k))/dx^2)+((T(i,j,k-1)-2*T(i,j,k)+T(i,j,k+1))/dz^2));
+
   end
     iter=iter+1;
     t=iter*dt; % pointless for Steady State
-    error=max(max(max(abs(Tss_new-Tss))));    
-    Tss=Tss_new;
+%     error=max(max(max(abs(T_new-T))));    
+  T(i,j,k)=T_new(i,j,k);
 
-if mod(iter,5)==0
-    plot3D(X,Y,Z,Lx,Ly,Lz,Tss,isovalues,iter,t); % e.g. plot every 5 iter
+if mod(iter,2)==0  % e.g. plot every 5 iter
+    plot3D(X,Y,Z,Lx,Ly,Lz,dx,dy,dz,T,isovalues,iter,t); 
 end
 
 if video >0 

plot3D.m

@@ -1,7 +1,15 @@
 
-function plot3D(X,Y,Z,Lx,Ly,Lz,T,isovalues,iter, t);
+function plot3D(X,Y,Z,Lx,Ly,Lz,dx,dy,dz,T,isovalues,iter, t);
 clf
 
+numisosurf=25; % number of isosurfaces
+
+%% recalculate scales (dsable this section to keep the static scale)
+ T_max = max(T,[],'all');
+ T_min = min(T,[],'all');
+ isovalues=linspace(T_min,T_max,numisosurf);
+
+ %%
 num=numel(isovalues);
  for i=1:num     
      p(i)=patch(isosurface(X,Y,Z,T,isovalues(i)));
@@ -19,10 +27,11 @@
 
 camproj perspective
 camlight; lighting gouraud; alpha(0.75);
-axis([0 Lx Ly/2 Ly 0 Lz]); % plot only half along Y to see inside
+axis([0 Lx 0 Ly 0 Lz]); % plot only half along Y to see inside
 axis manual
-view(iter*2,20);
-%view(3); switch on to have static view
+%view(iter,20); % orbiting camera
+%view(180,0); %switch on to have static view
+view(3) % iso view
 title(['Steady State Solution, iter = ' num2str(iter), ' time = ' num2str(t) ,' s']);
 drawnow