fulldynd.m

这是一个模拟第3类模式地震波的matlab脚本。 This a collection of Matlab scripts that solve the antiplane (mode III)

%fulldynD.m	2D Spectral BIEM for fault dynamics
%		Antiplane mode, DISPLACEMENT formulation
%		References: Morrisey and Geubelle 1997
%			    Lapusta et al. 2000
%
% INPUTS: see section "SIMULATION PARAMETERS"
% OUTPUTS: see section "STORE RESULTS"
% TO DO: adaptive timestep
%        eliminate replicate
% USES: C3.m, friction.m
%
% Jean-Paul Ampuero	[e-mail]	May 2003

%-----------------------------------------------
% UNITS
MPa = 1e6;
year = 356*24*60*60 ;

%-----------------------------------------------
% SIMULATION PARAMETERS

%-- physical properties --
L = 300. ; 	% fault length (simulated period)
CS = 3600.; 	% S wave velocity
RHO = 2800.;	% density
MU = RHO*CS^2; 	% shear modulus

%-- background stress --
TAUINIT1 = 60*MPa ;	% initial shear stress 
TAUINIT2 = 100 *MPa ;	% normal stress (+ compressive)
TAURATE1 = 0.2e8*MPa/year; % tectonic shear loading rate

%-- heterogeneity of initial shear stress --
HET.L = 100. ;	% length scale
HET.A = 5*MPa;	% amplitude

%-- numerical settings --
N = 128 ; 	% number of space elements, grid size
TMAX = 1300; 	% number of time steps

%-- friction --
FRIC.Dc = 8e-3;	% friction critical slip Dc
FRIC.MUs = 0.65;	% static friction coefficient
FRIC.MUd = 0.6;	% dynamic friction coefficient

%-- output settings --
TOUT = 10;	% time stride for outputs
XOUT = (N/4:3*N/4);	% space stride for outputs
XSIS = [N/2 3*N/8];	% coordinates of full time outputs
ECHO = 2;	% output info

%-- expert numerical settings --
RTCUT = 2 ;	% kernel_time-cut / fault_travel_time for mode_1, typically = 2
QW = 4 ; 	% kernel_time-cut Nyquist_mode / mode_1 , typically = 4
CFL = 0.5 ;	% stability number (Courant), typically = 0.5
RDELAY = 0.5 ; 	% kernel delay / dt , typically = 0.5

%-----------------------------------------------
% INITIALIZE

dx = L/N ;
dt = CFL*dx/CS ; 
x = (-N/2+1:N/2).' *dx ;
t = (1:TMAX).' *dt;

imped1 = MU/(2*CS) ;
dload1 = TAURATE1*dt ;
load1 = TAUINIT1 + HET.A*exp(-x.^2/HET.L^2);
stress2 = ones(N,1)*TAUINIT2 ;

% barrier:
NWEAK = N/2 ; % size of the central weak zone
RBAR = 10 ; % relative strength of the barrier
stress2(1:N/2-NWEAK/2) = RBAR*TAUINIT2 ;
stress2(N/2+NWEAK/2+1:end) = RBAR*TAUINIT2 ;

FRIC.W = (FRIC.MUs-FRIC.MUd) ./ FRIC.Dc ;

NK=N/2; % number of strictly positive wavenumbers
NKnyq = N/2+1 ; % Nyquist wavenumber position in classical FFT storage
delay = RDELAY*dt ;
TCUT1 = floor(RTCUT*N/CFL);

%--
% Kernel pre-integrated with trapezoidal rule, including delay
% Mode 0 has null contribution, don't need to store it,
% only strictly positive wavenumbers.
% For negative wavenumbers we will use symetries of real FFT
% Also implemented: KERNEL CUT-OFF (TCUT), 
% although this should be done only in the velocity formulation.
if ECHO > 1, disp('Computing the kernel'), end
dk = 2*pi/L ;
TCUT = ceil(  TCUT1*( 1+(QW-1)/(NK-1)*(0:NK-1) ) ./ (1:NK)  ); 
for ik=1:NK
  tc = TCUT(ik) ;
  k = dk*ik;
  times = (0:tc).'*dt+delay;
  arg = k*CS*times ;
  CIII = C3(arg) ;
  factor = -0.5*MU*CS*k^2 *dt/2 ;
  kernel(ik).f = factor*( CIII(1:tc) + CIII(2:tc+1) );
 % The following makes the DISPLACEMENT formulation with cut-off
 % EXACTLY equivalent to the VELOCITY formulation with cut-off
  WIII = 1- k*CS*( 0.5*CIII(1) +sum(CIII(2:end)) +0.5*CIII(end) )*dt ;
  kernel(ik).f(tc) =  kernel(ik).f(tc) -MU/2*k*WIII ;
end
clear CIII times arg WIII
%--

% fields
u = zeros(N,1);
v = zeros(N,1);
for ik=1:NK
  fu(ik).f = zeros(TCUT(ik),1);
end  
f1 = zeros(N,1);

% counters
ito = 1;
techo = ceil(TMAX/10);
ite = 1;
itcyc(1:NK) = 1;

%-----------------------------------------------
% BEGIN TIME LOOP
if ECHO > 1, disp('Begin time loop'), end

for it = 1:TMAX ,
  
% Some progress info  
  if ECHO > 1 & ite == techo
    ite = 1;
    disp(sprintf('Step %i/%i, Vmax = %0.4g, Umax = %0.4g' ...
                ,it,TMAX,max(v),max(u)))
  else
    ite = ite+1;
  end  

% Update strength
  frimu = friction(u,FRIC);
%  frimu(1:N) = 100*FRIC.MUs ; frimu(64-31:64+32) = FRIC.MUd; % crack test
  strength = stress2.*frimu ; % normal stress >0 is compressive
  
% Update external load  
  load1 = load1 + dload1 ;

% Trial state: assuming no further slip
  stress1 = load1 + f1 ;

% Solve friction
  v = max( stress1-strength ,0)/imped1 ;
  stress1 = min(stress1,strength) ;

%--- Begin STORE RESULTS ---
 % (x,t) fields, coarsely sampled
  if it == ito*TOUT
    usto(:,ito) = u(XOUT) ;
    vsto(:,ito) = v(XOUT) ;
    s1sto(:,ito) = stress1(XOUT) ;
    ststo(:,ito) = strength(XOUT) ;
    ito = ito+1;
  end  
 % fields(t) at selected locations, full time sampling 
  vsis(it,:) = v(XSIS)';
  usis(it,:) = u(XSIS)';
  s1sis(it,:) = stress1(XSIS)';
 % macroscopic fields, full time sampling
  um(it) = sum(u)/N ;
  vm(it) = sum(v)/N ;
  incrack = find(v>0);
  lm(it) = length(incrack); 
  sm(it) = sum(stress1)/N ;
%  smcrack(it) = sum( stress1(incrack) )*dx;
%--- End STORE RESULTS ---

  if it==TMAX, break, end

% Update slip for NEXT time step
  u = u+dt*v;
  fftu = fft(u);

  for ik=1:NK  

   % Update slip spectrum for NEXT time step
   % "fu" stores the slip spectrum history: fu(MODE).f(TIME)
   % The storage is CYCLIC in the time dimension:
   % don't need to keep slip history longer than TCUT(mode).
   % "itcyc" points, on entry, to the position of the current time step 
   %         in "fu(mode).f(:)"   = mod(it,TCUT(mode))+1
    tcut = TCUT(ik);
    itc = itcyc(ik)+1 ; if itc > tcut, itc = 1;, end 
    itcyc(ik) = itc; % now points to the next time step
    fu(ik).f(itc) = fftu(ik+1);

   % Update stress functional for NEXT time step
   % An alternative implementation (Lapusta et al. 2000)
   % avoids delay by using a predicted value for the current slip.
   % "ittab" points to the slip history from it+1 back to it-TCUT+2
    it2 = min( it+1, tcut );
    ittab = [ (itc:-1:1) (tcut:-1:itc+1) ];
    ff1(ik+1) = sum( kernel(ik).f(1:it2) .* fu(ik).f(ittab(1:it2)) ) ;

  end
  
  ff1(1) = 0 ; % mode 0 is null
  ff1(NKnyq+1:N) = conj( ff1(NK:-1:2) ); % real signal: complete the spectrum
  f1 = real(ifft(ff1.')) ;

end % END TIME LOOP

%-----------------------------------------------
% INFO

W = TAUINIT2*FRIC.W ; 	% slip weakening rate
Sm = W / imped1 ; 	% growth rate
Lc = 1.158*MU/W ;	% nucleation size

disp(sprintf('Space resolution Lc/dx = %0.4g',Lc/dx));
disp(sprintf( 'Time resolution 1/Sm*dt = %0.4g',1/(Sm*dt)) );