粘弹性介质解析解,什么是弹性介质,matlab

% MATLABScript for obtaning the analytical solution for 2D wave propagation % in a viscoelastic medium, based on the Appendix B of paper of Carcione et % al. (1988)(1), corrected in (2). % % The source used is a Ricker wavelet in the y velocity component. % % A figure of the source temporal function and its inverse Fourier % transform are provided to verify the good working of the method. Also % works for elastic variables by substituting the values of vp and vs by % real constant ones (vp_0 & vs_0). % % Attenuation is implemented by GMB-EK rheology presented in (4) and % thoroughfully explained in (3). % % In order to make a proper comparison with numerical data it is % recommended to either use the same rheology in both the numerical and the % analytical. It is also possible to get an almost constant Q value by % increasing a lot the number of mechanisms used in both the numerical and % the analytical solutions, thus minimizing the effects of the Q fitting in % the comparison. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% By Josep de la Puente, LMU Geophysics, 2005 %%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% home clear all; close all; disp('======================================') disp('COMPUTATION OF VISCOELASTIC ANALYTICAL') disp(' SOLUTION ') disp('======================================') disp(' ((c) Josep de la Puente Alvarez ') disp(' LMU Geophysics 2005)') disp('======================================') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Source parameters (Ricker wavelet): f0=10; % cutoff frequency t0=0.1; % time shift nu=0.5; % weight of exponential (for Carcione's source) epsilon=1.0; % weight of cosinus (for Carcione's source) %Frequency domain: w_max=1001; % Maximum frequency (minimum freq = -w_max). Is better to keep an odd number w_inc=0.1; % Increment in frequency %Position of receiver: x=0.5; y=0.5; %Force of the source: F=1.0; % Constant magnitude. Only alters amplitude of the seismograms %Material parameters: rho=1.0; % Density vp_0=sqrt(3); % P-wave velocity (elastic) vs_0=1; % S_wave velocity (elastic) %Attenuation Q-Solver parameters n=3; % Number of mechanisms we want to have QPval=20; % Desired QP constant value QSval=10; % Desired QS constant value freq=f0; % Central frequency of the absorption band (in Hertz). Good to center it % at the source's central frequency f_ratio=100; % The ratio between the maximum and minimum frequencies of our bandwidth % (Usually between 10^2 and 10^4) %Outputting variables toutmax=1.0; % Maximum value of t desired as output disp('--------------------------------------') disp(' GEOMETRY AND MATERIAL SETTINGS:'); disp('--------------------------------------') disp(strcat('Density= ',num2str(rho))); disp(strcat('VP= ',num2str(vp_0))); disp(strcat('VS= ',num2str(vs_0))); disp('Source at coordinates: (0,0)'); disp(strcat('Receiver at coordinates: (',num2str(x),',',num2str(y),')')); disp(' '); disp('--------------------------------------') disp(' ATTENUATION SETTINGS:'); disp('--------------------------------------') disp(strcat('Number of mechanisms= ',num2str(n))); disp(strcat('Q for P-wave= ',num2str(QPval))); disp(strcat('Q for S-wave= ',num2str(QSval))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% PROBLEM SETUP %%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Derived quantities initialization r=sqrt(x^2+y^2); % Distance to the receiver w0=f0/(2*pi); % Conversion to angular velocity w=[0:w_inc:w_max]; w=2*w-w_max; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% GMB-EK COMPLEX M SOLUTION %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Equation system initialization %% kmax=2*n-1; %Number of equations to solve (the system is overdetermined) AP=zeros(kmax,n); %Initialization of system matrix (P wave) AS=zeros(kmax,n); %Initialization of system matrix (S wave) QP=ones(kmax,1)/QPval; %Desired values of Q for each mechanism inverted (P wave) QS=ones(kmax,1)/QSval; % " " (S wave) YP=zeros(n,1); YS=zeros(n,1); %% Selection of the logarithmically equispaced frequencies wmean=2*pi*freq; wmin_disc=wmean/sqrt(f_ratio); for j=1:kmax w_disc(j)=exp(log(wmin_disc)+(j-1)/(kmax-1)*log(f_ratio)); end %% Filling of the linear system matrix %% for m=1:kmax for j=1:n AP(m,j)=(w_disc(2*j-1).*w_disc(m)+w_disc(2*j-1).^2/QPval)./(w_disc(2*j-1).^2+w_disc(m).^2); end end for m=1:kmax for j=1:n AS(m,j)=(w_disc(2*j-1).*w_disc(m)+w_disc(2*j-1).^2/QSval)./(w_disc(2*j-1).^2+w_disc(m).^2); end end %% Solving of the system %% YP=AP\QP; YS=AS\QS; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%% VISUALIZATION OF Q IN FREQUENCY %%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Getting values for the continuous representation of Q %% wmin=wmean/sqrt(f_ratio); wmax=wmean*sqrt(f_ratio); xfrq=[wmin:(wmax-wmin)/(10000-1):wmax]; %% P-wave Q continuous values numP=0; denP=1; for j=1:n numP=numP+(YP(j)*w_disc(2*j-1)*xfrq(:))./(w_disc(2*j-1)^2+xfrq(:).^2); denP=denP-(YP(j)*w_disc(2*j-1).^2)./(w_disc(2*j-1)^2+xfrq(:).^2); end Q_contP=denP./numP; %% S-wave Q continuous values numS=0; denS=1; for j=1:n numS=numS+(YS(j)*w_disc(2*j-1)*xfrq(:))./(w_disc(2*j-1)^2+xfrq(:).^2); denS=denS-(YS(j)*w_disc(2*j-1).^2)./(w_disc(2*j-1)^2+xfrq(:).^2); end Q_contS=denS./numS; %% Computing fitting quality (RMS and maximum difference) maxPdif=0; maxSdif=0; for j=1:length(Q_contP) tempP=abs(Q_contP(j)-QPval); if tempP >= maxPdif maxPdif=tempP; end tempS=abs(Q_contS(j)-QSval); if tempS >= maxSdif maxSdif=tempS; end end subplot(1,2,1), semilogx(xfrq/(2*pi),Q_contP,xfrq/(2*pi),QPval*ones(length(xfrq),1)), title('Q for the P-wave'),legend('computed','desired') xlabel('frequency (Hz)'),ylabel('Q'),axis([wmin/(2*pi) wmax/(2*pi) 0 QPval*1.25]) subplot(1,2,2), semilogx(xfrq/(2*pi),Q_contS,xfrq/(2*pi),QSval*ones(length(xfrq),1)), title('Q for the S-wave'),legend('computed','desired'), xlabel('frequency (Hz)'),ylabel('Q'),axis([wmin/(2*pi) wmax/(2*pi) 0 QSval*1.25]) disp(strcat('Attenuation bandwidth= [',num2str(wmin/(2*pi)),',',num2str(wmax/(2*pi)),'] Hz')); disp(''); disp(strcat('Maximum QP fitting error= ',num2str(maxPdif/QPval*100),' %')); disp(''); disp(strcat('Maximum QS fitting error= ',num2str(maxSdif/QSval*100),' %')); disp(' '); disp('--------------------------------------') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%% VELOCITIES COMPUTATION %%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Translating into Y_kappa and Y_mu %% YK=YP*0; YM=YP*0; YK(:)=(vp_0^2*YP(:)-4/3*vs_0^2*YS(:))/(vp_0^2-4/3*vs_0^2); YM(:)=YS(:); % Complex modulus mu_0=vs_0^2*rho; lam_0=vp_0^2*rho-2*mu_0; MUK=(lam_0+2/3*mu_0); % Elastic bulk modulus MUM=2*mu_0; % Elastic shear modulus MK=MUK*ones(length(w),1); MM=MUM*ones(length(w),1); lambda=zeros(length(w),1); mu=zeros(