%% ESS 590: HEAT AND MASS FLOW MODELLING SPRING 2008 %% %% Finite-Volume solution for Steady-State phi(x,z) %% diffusion with advection %% %% EXPLORE FALSE DIFFUSION %% %% 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 %% updated May 22, 2008 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % global c_x c_z c_prime % velocity parameterization % % % False Diffusion Test % bring in 2 streams at different phi on boundary % When Gamma is very small, the 2 streams should stay distinct % ------------------------------------------------------------ % Gamma_j = [ 0.1 0.1 0.01 0.001*[1 1 1 1] ]; Gamma_j = [ 0.1 0.1 0.01 0.0001*[1 1 1 1] ]; % % set up numbers of volumes in each direction % for each of several runs % decreasing volume size should decrease false diffusion % -------------------------------------------------------- % N_vol = [ 9 9 9 9 9 25 100 ]; N_vol = [ 10 10 10 10 10 50 1000 ]; % % set up velocity patterns % 1,2 velocity parallel to grid % same behavior with u=0 or w=0 % 3,4 velocity parallel to grid (w=0) % decreasing Gamma % 5,6,7 velocity at 45 degrees to grid, decreasing dx,dz % ------------------------------------------------------- c_x_vec = [ 0 sqrt(2)*[1 1 1] 1 1 1 ]; c_z_vec = [ sqrt(2) 0*[1 1 1] 1 1 1 ]; % % set values of phi in 2 incoming streams bc_phi_0 = 50; % deg C bc_phi_1 = 100; % deg C for i_count = 1:length(N_vol) % loop over different models pause % % set parameters for velocity calculations % u = c_x + c_prime * x % w = c_z - c_prime * z % -------------------------- c_x = c_x_vec(i_count); c_z = c_z_vec(i_count); 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 = N_vol(i_count); % number of finite volumes in x nx_P = nx_fv + 2; % number of nodes in x % % Lz = 10; % z-extent of domain nz_fv = N_vol(i_count); % 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 % --------------------------------------------------- % if(i_count == 1) % % Flow in z only % ---------------- %% phi gradients 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 gradients at west boundary bc_w_0 = 0; % deg m^-1 bc_type = 2; % gradient BC bc_w = repmat( [ bc_w_0 bc_type ], nz_P, 1 ); % % %% temperature phi_u on upper boundary % west half cold, east half hot n_cold = fix(nx_P/2) + 1; n_hot = nx_P - n_cold; % bc_type = 1; % specified-value BC %% bc_u = repmat( [ bc_phi_0 bc_type ]', 1, nx_P ); % bc_u = [ repmat( [ bc_phi_0 bc_type ]', 1, n_cold ) ... repmat( [ bc_phi_1 bc_type ]', 1, n_hot ) ]; % % %% 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 ); % % elseif( (i_count > 1) & (i_count < 5) ) % % % Flow in x only % ---------------- % 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 n_cold = fix(nx_P/2) + 1; n_hot = nx_P - n_cold; % bc_type = 1; % fixed phi BC % bc_w = [ repmat( [ bc_phi_0 bc_type ]', 1, n_cold ) ... repmat( [ bc_phi_1 bc_type ]', 1, n_hot ) ]'; % % %% phi gradient on upper boundary bc_u_0 = 0; % deg C m^-1 bc_type = 2; % specified-flux BC bc_u = repmat( [ bc_u_0 bc_type ]', 1, nx_P ); % % %% phi gradient on lower boundary bc_d_0 = 0; % deg C m^-1 bc_type = 2; % specified-flux BC bc_d = repmat( [ bc_d_0 bc_type ]', 1, nx_P ); % % else % % % flow is at 45 degres on 2-D domain % ------------------------------------ % 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 = bc_phi_1; % 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 = bc_phi_0; % deg C bc_type = 1; % specified-value BC bc_u = repmat( [ bc_u_0 bc_type ]', 1, nx_P ); %% bc_u_1 = 100; %% bc_u(1, 1:fix(nx_P/2)) = bc_u_1; % % %% 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 ); % end % if(i_count == 1) ... % % %% material parameters %%%%%%%%%%%%%%%%%%%%%%% % set up thermal-property matrix - uniform conductivity Gamma_0 = Gamma_j(i_count); 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(i_count) % % set up colors for plotting red = [ 1 0 0 ]; green = [ 0 1 0 ]; blue = [ 0 0 1 ]; magenta = [ 1 0 1 ]; black = [ 0 0 0 ]; color = [ red' blue' magenta' green' black']; % % % set up marker symbols for plotting marker = [ 'o' 'x' '*' '+' 'none ' ]; % % 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; if(i_count == length(N_vol) ) 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), ... 'linestyle', 'none' ) else 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) ) end %% shading interp colorbar hold on; grid on; box on set(gca,'zlim', [ bc_phi_0, bc_phi_1 ] ) xlabel('Distance (m)', 'FontSize', 12) ylabel('Depth (m)', 'FontSize', 12) zlabel('Temperature ( ^oC)', 'FontSize', 12) % % show parameters of current model display( ' ' ); display( ['Gamma = ' num2str(Gamma_0) ' dx = ' num2str(dx) ] ); display( [ 'Velocity (v_x, v_z) = (' num2str(c_x) ', ' num2str(c_z) ')' ] ); % P_x = rho_0*c_x*dx/Gamma_0; P_z = rho_0*c_z*dz/Gamma_0; display( [ 'Grid Peclet numbers (P_x, P_z) = (' num2str(P_x) ', ' num2str(P_z) ')' ]); % end % for_i_count = 1:length(N_vol) ... % % % % %