%% ESS 590D: HEAT AND MASS FLOW MODELLING SPRING 2006 %% %% Finite Volume solution for 1-D or 2_D transient phi(x,z) %% on a rectangle %% compared with Finite Difference solution %% %% div( Gamma grad phi) + S = d(rho phi)/dt %% %% For example, %% phi = internal energy (J kg^-1) %% = int_0^T c(T') dT' %% T = temperature (deg) %% Gamma = K/c %% K = thermal conductivity (W m^-1 deg^-1) %% c = specific heat capacity (J kg^-1 deg^-1) %% S = sources (W m^-3) %% + S_C + S_P phi %% rho = density %% %% THIS CODE (what's special?) %% %% 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 S = S_C + S_P * T %% %% 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. %% %% %% This main program sets and delivers: %% - Mesh geometry (X_P, Z_P, X_edges, Z_edges) %% - time-marching info (nt, dt) %% - constitutive properties (K_P) %% - sources (S_C, S_P) %% - initial conditions (phi_0) %% - boundary conditions (phi_u, q_d, q_e, q_w) %% %% %% Ed Waddington April 2006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %% set up finite volumes %%%%%%%%%%%%%%%%%%%%%%%%%% % This code sets up uniform-sized volumes in both space and time % (could be generalized) % % Space dimensions % ----------------- Lx = 20; % x-extent of domain % nx = 21; % number of control volumes in x dx = Lx/max((nx-1),1 ); % x-widths of volumes % max function prevents division by zero % in 1-D (i.e. when nx=1) % Lz = 20; % z-extent of domain % nz = 21; % 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); % % %% set up matrices X_edges, Z_edges - defining (x,z) coordinates %% of finite-volume interfaces 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,:)' ]'; % % % Time dimension % --------------- sec_per_year = 365.25*24*3600; %% Lt = 5; % Total duration of time integration %% nt = 101; % number of future points Lt = 1; % Total duration of time integration nt = 51; % number of future points dt_yr = Lt/max((nt-1),1); % time step in years dt = dt_yr * sec_per_year; % % %% Boundary Conditions %%%%%%%%%%%%%%%%%%%%%%%%% %% flux through east, west and down boundaries % q_d is negative because heat flows toward decreasing depths q_d_0 = -0.50; % W m^-2 q_e_0 = 0.0; % W m^-2 %% q_w_0 = -0.50; % W m^-2 q_w_0 = 0.0; % W m^-2 % % This code sets of simple steady B.C.s % put "interesting" transient B.C.s into these matrices q_d = repmat( q_d_0, nt, nx ); q_w = repmat( q_w_0, nz, nt ); q_e = repmat( q_e_0, nz, nt ); % % %% temperature phi_u on upper boundary %% put your favorite pattern here % ------------------------------------- phi_u_0 = -20; % deg C % % constant and uniform surface temperature %% phi_u = repmat( phi_u_0, nt, nx ); % % seasonally varying but spatially uniform phi_range = 5; time_vec = dt_yr * (0:(nt-1)); phi_seasonal = phi_u_0 + phi_range * sin( 2*pi * time_vec ); phi_u = repmat( phi_seasonal', 1, nx ); % % %% material parameters %%%%%%%%%%%%%%%%%%%%%%% % set up thermal-property matrix - e.g. uniform conductivity K_0 = 2.0; % W m^-1 deg^-1 K_P = repmat( K_0, size(X_P) ); % % % set up density matrix - e.g. uniform on domain rho_0 = 900; % kg m^-3 (e.g. ice) rho_P = repmat( rho_0, size(X_P) ); % % % set up specific heat capacity - uniform on domain c_0 = 2000; % J kg^-1 deg^-1 c_P = repmat( c_0, size(X_P) ); % % % set up generic Gamma matrix - for heat flow, Gamma = K/c %% Gamma_P = K_P./c_P; rho_c_P = rho_P.* c_P; % %% source term %%%%%%%%%%%%%%%%%%%%% % S = S_C + S_P * phi_P % If you want transient source terms, you need to build that information % into a transient vector S(t), or (storage expensive) 3-D arrays % S_C(x,z,t) and S_P(x,z,t) % spatially uniform source term, no dependence on solution phi S_C_0 = 0; % W m^-3 s^-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 % % % % initial condition % ----------------- %% phi_0 = zeros( size(X_P) ); %% phi_0 = repmat( phi_u, nz,1 )... %% - [ Z_P(1,:)' cumsum(diff(Z_P))']'.* repmat(q_d, nz,1)./ K_P; % % For example, set constant initial value, equal to BC on upper boundary phi_0 = repmat( phi_u(1,:), nz, 1 ); % % % This is it! solve for phi(x,z,t) using Finite Volumes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi = fvm_soln_transient_bdry_nodes( X_edges, Z_edges, X_P, Z_P, ... nt, dt, rho_c_P, K_P, S_C, S_P, ... phi_0, phi_u, q_d, q_e, q_w ); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % dt rho_c_dz_2_over_2K = rho_c_P(1,1)*(Z_P(2,1)-Z_P(1,1)).^2/(2 * K_P(1,1)) % %% Plot results %%%%%%%%%%%%%%%%% % if(nx <= 2) % %% 1-D line graphs %%%%%%%%%%%%%%%%%% % % plot conductivity (K_z) % ----------------------- figure(1) hold on; grid on; box on plot( Z_P(:,1), K_P(:,1), 'ro' ) ylabel('Conductivity (W m^{-1} deg^{-1})') xlabel('Depth (m)') title('Conductivity K_P') % plot source (S_C) % ----------------- figure(2) hold on; grid on; box on plot( Z_P(:,1), S_C(:,1), 'rx' ) ylabel('Source S_C (W m^{-3})') xlabel('Depth (m)') title('Source Term S_C') % % plot source (S_C) % ----------------- figure(3) hold on; grid on; box on plot( Z_P(:,1), S_P(:,1), 'bx' ) ylabel('Source S_P (W m^{-3} deg^{-1})') xlabel('Depth (m)') title('Source Term S_P') % % 1-D, plot phi(z) solution as line graph % --------------------------------------- figure(4) plot(Z_P(:,1), phi(:,1), 'bo-' ) hold on; grid on; box on xlabel('Depth (m)') ylabel('Temperature (deg C)') title('Finite Volume Method') % % % else % %% plot in 2-D %%%%%%%%%%%%%%%%%%% % % plot conductivity (K_z) % ----------------------- figure(1) hold on; grid on; box on surf(X_P, Z_P, K_P ) colorbar zlabel('Conductivity (W m^{-1} deg^{-1})') xlabel('Depth (m)') title('Conductivity K_P') % plot source (S_C) % ----------------- figure(2) hold on; grid on; box on surf(X_P, Z_P, S_C ) colorbar zlabel('Source S_C (W m^{-3})') xlabel('Distance (m)') ylabel('Depth (m)') title('Source Term S_C') % % plot source (S_P) % ----------------- figure(3) hold on; surf(X_P, Z_P, S_P ) colorbar zlabel('Source S_P (W m^{-3} deg^{-1})') xlabel('Distance (m)') ylabel('Depth (m)') title('Source Term S_P') % % 2-D, plot phi(x,z) solution as perspective surface % -------------------------------------------------- figure(4) hold on surf(X_P, Z_P, phi ) colorbar xlabel('Distance (m)') ylabel('Depth (m)') title('Control Volume Method') end % if( nx <= 2 ) ... % % % %