%% ESS 590D: HEAT AND MASS FLOW MODELLING SPRING 2006 %% %% Finite-Volume solution for Steady-State phi(x,z) %% diffusion with advection %% %% div( rho u phi ) = div( Gamma grad phi) + S %% %% For example, %% phi = 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 %% %% This code uses PATANKAR POWER-LAW SCHEME for advective fluxes %% %% Externally specified velocity field (u,w) must satisfy %% mass conservation %% %% div( rho u) = 0 %% %% Boundary Conditions: %% All boundaries are finite-volume interfaces %% (Patankar Practice B, page 69) %% However, there is an additional "point" X_B located on boundary. %% Its value and coefficients are determined by phi value or gradient %% boundary conditions there. %% Adjacent interior point is X_I %% (X_I - X_B) = dX_B %% %% phi_B = phi_0 %% a_B = 1, all other a_J = 0 %% b = phi_0 %% %% d phi/dx = dphi_dx_0 %% use gradient between X_B and X_I (interior point) %% d phi/dx|B = Gamma_I (phi_I - phi_B)/dx_B %% or %% phi_B = dphi_dx_0 * dx_B + phi_I %% %% For each boundary there is a matrix bc_x %% with either 2 rows (Up,Down) or 2 columns (East,West). %% e.g. East %% bc_e is [nz_P x 2] where nz_P is number of nodes on boundary %% - first column contains value of BC at that node on boundary %% - second column specifies type of BC %% - 1 -> specified phi %% - 2 -> specified phi gradient across boundary %% %% Similarly bc_w, bc_u, bc_d on other boundaries %% %% fvm_SS_adv_diff_2D.m will correctly apply any combination of BC %% that are passed to it in matrices bc_e etc %% %% %% Material Properties %% Gamma(x,z) and rho((x,z) can vary 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 May 7 2006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % global c_x c_z c_prime % velocity parameterization % % % set parameters for velocity calculations % u = c_x + c_prime * x % w = c_z - c_prime * z % -------------------------- c_x = 1; c_z = 1; c_prime = 0; % % %%%%%%%%%%%%%%%%%%%%%%%%%% % set up finite volumes %%%%%%%%%%%%%%%%%%%%%%%%%% % This code sets up uniform-sized volumes (could be generalized) % There are points on all boundaries, which are also % volume interfaces (Patankar Practice B, page 69) % Lx = 10; % x-extent of domain nx_fv = 20; % number of finite volumes in x nx_P = nx_fv + 2; % number of nodes in x % % Lz = 10; % z-extent of domain nz_fv = 20; % number of finite volumes in z nz_P = nz_fv + 2; % number of nodes in z % % % set up vectors x_edges_vec, z_edges_vec - % defining (x,z) coordinates of finite-volume edges % % set x-widths of volumes % (this bit of code sets up uniform edge spacing) % ------------------------------------------------- dx = Lx/nx_fv; % z-widths of volumes x_edges_vec = dx * (0:nx_fv); % %% X_edges = repmat(x_edges_vec, nz_fv+1, 1); % % % set x-widths of volumes % (this bit of code sets up uniform edge spacing) % ------------------------------------------------- dz = Lz/nz_fv; % z-widths of volumes z_edges_vec = dz * (0:nz_fv)'; % %% Z_edges = repmat(z_edges_vec, 1, nx_fv+1); % % set nodes at finite-volume centers and on all boundaries x_P_vec = [ x_edges_vec(1) ... (x_edges_vec(1:end-1) + diff(x_edges_vec)/2) ... x_edges_vec(end) ]; % % z_P_vec = [ z_edges_vec(1) ... (z_edges_vec(1:end-1) + diff(z_edges_vec)/2)' ... z_edges_vec(end) ]'; % % %%%%%%%%%%%%%%%%%%%%%%%%% %% Boundary Conditions %%%%%%%%%%%%%%%%%%%%%%%%% %% % Boundary conditions along each boundary are carried in an array % labelled bc_x. For each boundary point, % - first entry is numerical value of BC, % - second entry indicates type of BC % 1 for specified value % 2 for specified gradient % --------------------------------------------------- % % phi gradient at east boundary bc_e_0 = 0; % deg m^-1 bc_type = 2; % gradient BC bc_e = repmat( [ bc_e_0 bc_type ], nz_P, 1 ); % % %% phi at west boundary bc_w_0 = 25; % deg m^-1 bc_type = 1; % gradient BC bc_w = repmat( [ bc_w_0 bc_type ], nz_P, 1 ); % % %% phi on upper boundary bc_u_0 = 50; % deg C bc_type = 1; % specified-value BC bc_u = repmat( [ bc_u_0 bc_type ]', 1, nx_P ); % % %% phi gradient on lower boundary bc_d_0 = 0; % deg C bc_type = 2; % specified-value BC bc_d = repmat( [ bc_d_0 bc_type ]', 1, nx_P ); % % %% material parameters %%%%%%%%%%%%%%%%%%%%%%% % set up thermal-property matrix - uniform conductivity Gamma_0 = 1; % W m^-1 deg^-1 Gamma_P = repmat( Gamma_0, nz_P, nx_P ); % % % set up density matrix rho_0 = 1; % kg m^-3 rho = repmat( rho_0, nz_P, nx_P ); % % %% source term %%%%%%%%%%%%%%%%%%%%% % S = S_C + S_P * phi_P % ---------------------- 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, nz_P, nx_P ); % constant source term S_P = repmat( S_P_0, nz_P, nx_P ); % phi_P-dependent source % % % This is it! solve for phi using Finite Volumes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% phi = fvm_SS_adv_diff_2D( x_edges_vec, z_edges_vec, ... x_P_vec, z_P_vec, Gamma_P, rho, S_C, S_P, ... bc_u, bc_d, bc_e, bc_w ); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % %%%%%%%%%%%%%%%%% %% Plot results %%%%%%%%%%%%%%%%% % figure(1) % % % recover nodal coordinates X_P = repmat(x_P_vec, nz_P, 1); Z_P = repmat(z_P_vec, 1, nx_P); % % 2-D, plot phi(x,z) solution as perspective surface % -------------------------------------------------- %% surf(X_P, Z_P, phi ) i_skip = 1; % decimate grids for plotting? % surf(X_P(1:i_skip:end, 1:i_skip:end), ... Z_P(1:i_skip:end, 1:i_skip:end), ... phi(1:i_skip:end, 1:i_skip:end) ) % colorbar hold on; grid on; box on xlabel('Distance (m)') ylabel('Depth (m)') zlabel('Temperature (deg C)') % Gamma_0__c_x__c_z__dx = [ Gamma_0 c_x c_z dx ] % % % % %