%% ESS 590D: HEAT AND MASS FLOW MODELLING SPRING 2006 %% %% Finite Volume solution for 1-D or 2_D steady phi(x,z) %% on a rectangle %% compared with Finite Difference solution %% %% 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 COMPARES FDM AND FVM SOLUTIONS TO ANALYTICAL SOLUTION %% BASED ON POLYNOMIALS FOR CONDUCTIVITY AND SOURCE TERM %% %% 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 varies with position (x,z) %% %% Source %% Source is independent of phi, but S_C varies with position %% %% 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 Control Volumes %% fvm = 0; % use Finite Differences % % set up control 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 % % 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 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) 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) ); % if(nx <= 2) % (2 half-volumes) % % set up inverse-polynomial conductivity % -------------------------------------- %% z_K = z_edges_vec; %% K_z = repmat(K_0, size(z_K) ); z_K = [ 0 0.35 0.7 1 2 3 4 5 6.3 7 8 9 9.5 10 ]; K_z = [ 2 1.7 1.4 1 3 1 1 2 2 3 1 2.5 1.5 0.5 ]; % % get polynomial through K_z(z_K) (check it doesn't go negative!) P_K_inv = polyfit( z_K, 1./K_z, (length(z_K)-1) ); % get polynomial values at volume centers K_P = 1./polyval( P_K_inv, Z_P ); end % if(nx <= 2) % %% source term %%%%%%%%%%%%%%%%%%%%% % S = S_C + S_P * phi_P % spatially uniform source term, no dependence on solution phi S_C_0 = 0; % W m^-3 s^-1 %% S_P_0 = 10; % W m^-3 s^-1 deg^-1 S_P_0 = 0; % 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 % %% if(nx == 1) if(nx <= 2) %% if(nx <= 0) % set up polynomial source %% z_S = z_edges_vec; %% S_z = repmat(S_P_0, size(z_S), 1 ); z_S = [ 0 0.5 1 2 3 4 5 6 7 8.5 9.7 10 ]; S_z = 0.1 * [ 2 1.5 1 3 1 1 2 2.1 1.5 3 3.6 0 ]; % get polynomial fitting S_z(z_S) %% P_S = polyfit( z_S, S_z, 2 ); P_S = polyfit( z_S, S_z, (length(z_S)-1) ); % % get polynomial source term S_C at control volume centers S_C = polyval( P_S, Z_P ); end % if(nx == 1) % % 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 ) ... % if(nx <= 2) %% find analytical solution in 1-D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% zz = linspace(z_P_vec(1), z_P_vec(end), 200 ); % integrate polynomial for source term S_C P_S_int = polyint( P_S ); % % multiply S_int and 1/K to get int_x1^x [S(x')] dx' / K(x) S_int_over_K = polyval(P_S_int, zz).* polyval(P_K_inv, zz); % % find polynomial representing S_int_over_K n_poly = length(P_S_int) + length(P_K_inv) - 1; P_S_int_over_K = polyfit( zz, S_int_over_K, n_poly); % % now integrate it again P_S_int_over_K_int = polyint( P_S_int_over_K ); % % integrate polynomial for inverse conductivity P_K_inv_int = polyint( P_K_inv ); % phi_analytic = phi_u_0 - ( polyval(P_S_int_over_K_int, zz ) ... - polyval(P_S_int_over_K_int, zz(1) ) ) ... + ( polyval(P_K_inv_int, zz ) ... - polyval(P_K_inv_int, zz(1) ) ) ... .* ( polyval(P_S_int, zz(end) ) - q_d_0 ); end % if(nx <= 2) ... %% Plot results %%%%%%%%%%%%%%%%% % if(fvm == 1) figure(1) % Finite-Volume solution else figure(2) % Finite-Difference solution end % if( nx <= 2 ) % 1-D, plot phi(z) solution as line graph % --------------------------------------- plot(Z_P(:,1), phi(:,1), 'bo-' ) hold on; grid on; box on % analytical solution plot(zz, phi_analytic, 'r-') xlabel('Depth (m)') ylabel('Temperature (deg C)') if(fvm == 1) title('Finite Volume Method') legend( 'blue - fvm', 'Red - analytical', 'location', 'northwest') else title('Finite Difference Method') legend( 'blue - fdm', 'Red - analytical', 'location', 'northwest') 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 <= 2 ) ... % if(nx <= 2) % plot conductivity (K_z) % ----------------------- figure(3) hold on; grid on; box on plot( Z_P(:,1), 1./polyval( P_K_inv, Z_P(:,1) ), 'ro' ) plot( zz, 1./polyval( P_K_inv, zz ), 'b-' ) ylabel('Conductivity (W m^{-1} deg^{-1})') xlabel('Depth (m)') title('Polynomial Conductivity K_P(z)') % plot source (S_C) % ----------------- figure(4) hold on; grid on; box on plot( Z_P(:,1), S_C(:,1), 'gx' ) plot( Z_P, polyval( P_S, Z_P(:,1) ), 'ro' ) plot( zz, polyval( P_S, zz ), 'b-' ) ylabel('Source S_C (W m^{-3})') xlabel('Depth (m)') title('Polynomial Source Term S_C(z)') end % if(nx <= 2) ... % % % % %