%% ESS 590D: HEAT AND MASS FLOW MODELLING SPRING 2006 %% %% THIS IS A WAVE-EQUATION SOLVER FOR SH WAVES %% %% d^2(rho phi)/dt^2 = div( Gamma grad phi) %% %% phi = displacement normal to propagation plane %% Gamma = elastic shear modulus %% rho = density %% %% This code is adapted from transient diffusion code %% The only differences are in the transient term. %% Derivation closely follows solution of momentum equations %% to derive Navier Stokes equations. %% In elasticity, displacement is the variable that transmits %% momentum into a volume through shear modulus, rather than %% velocity transitting momentum through viscosity. %% In order to have displacement as the variable inteh transient %% term, we have to write the velocity in %% %% [ (rho v)_P - (rho v)_P^0 ] dx dz %% %% as the _rate of change_ of displacement u, giving %% %% [ (rho u_dot)_P - (rho u_dot)_P^0 ] dx dz %% = [ rho( (u_P - u_P^0)/dt - (u_P^0 - u_P^-1)/dt ] dx dz %% %% where u_P^0 and u_P^-1 are the (known) solution at node P at %% previous and next previous times. %% After dividing through both sides of discretization equation by dt, %% we still have a coefficient a_P_0, except that instead of being %% a_P^0 = rho dx dz/dt %% it is now %% a_P^0 = rho dx dz/dt^2 %% The coefficients a_E, a_W etc are unchanged from the %% diffusion problem. %% The right-hand-side term b of known constants, %% rather than having [ a_P^0 * u_P^0 ] has instead %% [ 2 * a_P^0 * u_P^0 - a_P^0 * u_P^-1 ] %% %% Caution: Although this scheme follows all Patankar's rules, %% results can be inaccurate if: %% - mesh is too coarse to capture shape of wave. %% - time step is too large to ensure that temporal curvature %% (centered at u_P^0) can balance spatial curvature %% (centered on u_P). %% Large temporal mismatch in centers of curvature causes %% dissipation, and waves die out. %% %% %% Boundary Conditions: %% All boundaries are on finite-volume interfaces %% (Patankar Practice B, page 69) %% However, there is an additional "node" X_B located on each % boundary-volume face. %% 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 %% %% %% vector bc_e - specified phi gradient on east boundary %% vector bc_w - specified phi gradient on west boundary %% vector bc_d - specified phi gradient on lower boundary %% vector bc_u - specified phi gradient on upper boundary %% %% There must be at least 1 place where a value of phi is set. %% This is achieved by setting a phi "source" in the upper boundary. %% In the region "source_spot", phi BC is a single cycle of a sine wave %% with amplitude "source_height" and period "source_period" in time. %% By making "source_spot" longer or shorter it is possible to %% simulate everything from plane waves to "point" sources. %% %% Tis code can run 4 different examples, based on the user-set %% value of "fig_num" %% fig_num = 1 plane wave into homogeneous domain %% fig_num = 2 circular (point source) wave in homogeneous domain %% fig_num = 3 plane wave into domain with an internal interface %% fig_num = 4 plane wave into domain with a fluid inclusion %% %% %% %% Ed Waddington May 24 2006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % global source_spot source_period source_height global dt_plot fig_num % % seismic source parameters % source_spot - node index on upper boundary for line source % (will be set after grid has been set) % source_period - period of single sine-wave cycle % source_height - amplitude of displacement source % ----------------------------------------------------- source_period = 1; source_height = 1; % % fig_num determines nature of model % fig_num = 1 plane wave into homogeneous domain % fig_num = 2 circular (point source) wave in homogeneous domain % fig_num = 3 plane wave into domain with an internal interface % fig_num = 4 plane wave into domain with a fluid inclusion % --------------------------------------------------------------- fig_num = 2; % %%%%%%%%%%%%%%%%%%%%%%%%%% % 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 = 2; % x-extent of domain nx_fv = 50; % number of finite volumes in x nx_P = nx_fv + 2; % number of nodes in x % % Lz = 2; % z-extent of domain nz_fv = 60; % number of finite volumes in z nz_P = nz_fv + 2; % number of nodes in z % % % set up matrices X_edges, Z_edges - 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); % 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) ]; % X_P = repmat(x_P_vec, nz_P, 1); % % z_P_vec = [ z_edges_vec(1) ... (z_edges_vec(1:end-1) + diff(z_edges_vec)/2)' ... z_edges_vec(end) ]'; % Z_P = repmat(z_P_vec, 1, nx_P); % % % %% 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; bc_type = 2; bc_e = repmat( [ bc_e_0 bc_type ], nz_P, 1 ); % %% phi gradient at west boundary bc_w_0 = 0; bc_type = 2; bc_w = repmat( [ bc_w_0 bc_type ], nz_P, 1 ); % % %% phi gradient on upper boundary bc_u_0 = 0; bc_type = 2; bc_u = repmat( [ bc_u_0 bc_type ]', 1, nx_P ); % % add displacement source on upper boundary at source_spot if(fig_num == 1) % plane wave in uniform medium source_spot = 1:nx_P; % elseif(fig_num == 2) % "point" source (3 volumes) %% source_spot = fix(nx_P/2) + 1; source_spot = [ (fix(nx_P/2) - 1):(fix(nx_P/2) + 1) ]; % elseif(fig_num == 3) % plane wave hits sloping interface source_spot = 1:nx_P; % elseif(fig_num == 4) % plane wave hits fluid inclusion source_spot = 1:nx_P; end % bc_u(2, source_spot) = 1; % % % %% phi gradient on lower boundary bc_d_0 = 0; bc_type = 2; bc_d = repmat( [ bc_d_0 bc_type ]', 1, nx_P ); % % % % set time-step control parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dt = 0.002; nt = 2000; % dt_plot = 0.1; % set up initial condition % ------------------------ phi_0 = zeros( nz_P, nx_P ); % % eliminate values on corners where there are no nodes bc_e(1,1) = NaN; bc_e(1, end) = NaN; bc_w(1,1) = NaN; bc_w(1, end) = NaN; bc_u(1,1) = NaN; bc_u(1, end) = NaN; bc_d(1,1) = NaN; bc_d(1, end) = NaN; phi_0(1,1) = NaN; phi_0(1,end) = NaN; phi_0(end,1) = NaN; phi_0(end,end) = NaN; % %% 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 interfaces if(fig_num == 1) % plane wave in uniform medium % Gamma_P is set to go; % elseif(fig_num == 2) % "point" source (3 volumes) %% source_spot = fix(nx_P/2) + 1; % Gamma_P is set to go; % elseif(fig_num == 3) % plane wave hits sloping interface % slower wave speed below interface % Z = z_0 * (1 - X) + z_1 * X %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Gamma_1 = 0.2; % z_range = z_P_vec(end) - z_P_vec(1); z_0 = 0.9; z_1 = 0.2; % x_range = x_P_vec(end) - x_P_vec(1); x_0 = 0; x_1 = 1; for ix = 1:nx_P X_bdry = x_P_vec./ x_range; Z_bdry = (z_0 * (1 - X_bdry) + z_1 * X_bdry)'; index = find( Z_P(:,ix)/z_range >= Z_bdry(ix) ); Gamma_P( index, ix) = Gamma_1; end % for ix = 1:nx_P ... % elseif(fig_num == 4) % plane wave hits fluid inclusion at some internal points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % shear modulus of near-fluid inclusion Gamma_1 = 0.0001; x_range = x_P_vec(end) - x_P_vec(1); z_range = z_P_vec(end) - z_P_vec(1); % center and radius of inclusion x_center = x_range/2; z_center = z_range/2; radius = x_range/6; % non-dimensional distance between each node and inclusion center R_2 = ( (X_P -x_center)./radius ).^2 + ... ( (Z_P -z_center)./radius ).^2; % % put Gamma_1 into volumes inside circle index = find( R_2 <= 1); Gamma_P(index) = Gamma_1; % end % if( fig_num == ... % % 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_soln_SH_waves_2D( x_edges_vec, z_edges_vec, ... x_P_vec, z_P_vec, nt, dt, Gamma_P, rho, S_C, S_P, ... bc_u, bc_d, bc_e, bc_w, phi_0 ); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % %