%% ESS 590D: HEAT AND MASS FLOW MODELLING SPRING 2006 %% %% Control Volume solution for 1-D or 2_D steady phi(x,z) %% on a rectangle %% %% div( K grad phi) + S = 0 %% or %% del^2( phi ) + (grad K) dot (grad phi ) + S = 0 %% %% For example, %% phi = temperature (deg) %% K = thermal conductivity (W m^-1 deg^-1) %% S = sources (W m^-3) %% + S_C + S_P phi %% %% THIS CODE SETS UP POSITIVE SOURCE FEEDBACK S_P %% %% Boundary Conditions follow Patankar's approach - half-volumes on %% boundaries %% - Specified quantity phi_u(x) on upper surface %% - specified fluxes q_e, q_w, and q_d through other 3 sides of domain %% %% Material Properties %% K can vary with position (x,z) %% %% In the spirit of matlab's vector capabilities, the program %% avoids loops over spatial indices. %% Instead, matrices are formed for all spatially varying parameters %% and variables. %% %% Ed Waddington April 2006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Choose Finite Volumes or Finite Differences % -------------------------------------------- fvm = 1; % use Finite Volumes %% fvm = 0; % use Finite Differences % % set up Finite volumes %%%%%%%%%%%%%%%%%%%%%%%%%% % This code sets up uniform-sized volumes % (could be generalized) % Lx = 1; % x-extent of domain % nx = 2; % number of control volumes in x dx = Lx/max((nx-1),1 ); % x-widths of volumes % %% Lz = 10; % z-extent of domain Lz = 4 * pi; % z-extent of domain % % nz = 6; % number of control volumes in z % nz = 11; % number of control volumes in z % nz = 21; % number of control volumes in z nz = 41; % number of control volumes in z % nz = 81; % number of control volumes in z % nz = 161; % number of control volumes in z dz = Lz/max((nz-1),1 ); % z-heights of volumes %% set up matrices X_P, Z_P - defining (x,z) coordinates % of finite-volume centers x_P_vec = dx * (0:(nx-1) ); z_P_vec = dz * (0:(nz-1) )'; X_P = repmat(x_P_vec,nz,1); Z_P = repmat(z_P_vec,1,nx); % % % use Finite-Volume method % if(fvm == 1) X_edges = [ X_P(:,1) (X_P(:, 1:end-1) + diff(X_P')'/2) ... X_P(:,end) ]; % Z_edges = [ Z_P(1,:)' (Z_P(1:end-1, :) + diff(Z_P)/2)' ... Z_P(end,:)' ]'; % else % use Finite-Difference method % (don't ever really need X_edges, Z_edges with fdm) X_edges = X_P; Z_edges = Z_P; end % %% Boundary Conditions %%%%%%%%%%%%%%%%%%%%%% %% flux through east, west and down boundaries % q_d is negative because heat flows toward decreasing depths %% q_d_0 = -0.10; % W m^-2 q_d_0 = 0.0; % W m^-2 q_e_0 = 0.0; % W m^-2 q_w_0 = 0.0; % W m^-2 % q_d = repmat( q_d_0, 1, nx ); q_w = repmat( q_w_0, nz, 1 ); q_e = repmat( q_e_0, nz, 1 ); % % %% temperature phi_u on upper boundary %% phi_u_0 = -20; % deg C phi_u_0 = 10; % deg C phi_u = repmat( phi_u_0, 1, nx ); % % %% material parameters %%%%%%%%%%%%%%%%%%%%%%%%% % set up thermal-property matrix - uniform conductivity %% K_0 = 2.0; % W m^-1 deg^-1 K_0 = 1; % W m^-1 deg^-1 K_P = repmat( K_0, size(X_P) ); % % %% source term %%%%%%%%%%%%%%%%%%%%% % S = S_C + S_P * phi_P % spatially uniform source term, no dependence on solution phi %% S_C_0 = 10; % W m^-3 s^-1 %% S_P_0 = 0; % W m^-3 s^-1 deg^-1 % % source term depends on solution phi (positive S_P) S_C_0 = 0; % W m^-3 s^-1 S_P_0 = 1; % W m^-3 s^-1 deg^-1 % % S_C = repmat( S_C_0, size(X_P) ); % constant source term S_P = repmat( S_P_0, size(X_P) ); % phi_P-dependent source % S_P_0 two_K_over_dz_squared = 2*K_0/dz^2 % if (fvm == 1) % % This is it! solve for phi using Finite Volumes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi = fvm_soln_SS_bdry_nodes( X_edges, Z_edges, X_P, Z_P, K_P, ... S_C, S_P, phi_u, q_d, q_e, q_w ); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % else % % % This is it! solve for phi using Finite Differences %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi = fdm_soln_SS( X_edges, Z_edges, X_P, Z_P, K_P, S_C, S_P, ... phi_u, q_d, q_e, q_w ); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % end % if (fvm == 1 ) ... % %% Plot results %%%%%%%%%%%%%%%%% % if(fvm == 1) figure(1) % Control-Volume solution else figure(2) % Finite-Difference solution (under construction) end % if( nx <= 2 ) % 1-D, plot phi(z) solution as line graph % --------------------------------------- %% plot(Z_P(:,1), phi(:,1), 'bo-' ) plot(Z_P(:,1), phi(:,1), 'ro-' ) hold on; grid on; box on %% plot analytical solution zz = linspace(0, Lz); % if( S_P_0 >= 0) beta = sqrt(S_P_0/K_0); %% (assumes q_d=0, S_P_0=1, and Lz is multiple of 2 pi)) phi_zz = phi_u(1,1) * cos(beta * zz); else beta = sqrt(-S_P_0/K_0); b = phi_u(1,1)/(1 + exp(2*Lz) ); a = phi_u(1,1) - b; phi_zz = a*exp(-beta.*zz) + b*exp(beta.*zz); end % plot(zz,phi_zz, 'g') % xlabel('Depth (m)') ylabel('Temperature (deg C)') if(fvm == 1) title('Finite Volume Method') else title('Finite Difference Method') end % else % 2-D, plot phi(x,z) solution as perspective surface % -------------------------------------------------- surf(X_P, Z_P, phi ) colorbar hold on; grid on; box on xlabel('Distance (m)') ylabel('Depth (m)') zlabel('Temperature (deg C)') if(fvm == 1) title('Control Volume Method') else title('Finite Difference Method') end end % if( nx <= 1 ) ... % % % % %