我试图实现FD在2 d交错网格方案按(Virieux 1986)。我有困难从数值分析公式,这个方案的代码。我没有看到任何交错网格FD的具体的例子。
我希望看到一个例子。如果任何人都可以链接一本书,纸或网页,这将是极大的赞赏。
欢呼,
1回答
\ begingroup美元
\ endgroup美元
这是一个简单的MATLAB脚本为一个简单的FD实现(在时间和空间二阶准确)的声波方程,这是相同的Virieux(1986)的假设下\μ= 0美元和$ \ sigma_ {xx} = \ sigma_ {zz} = - p $压力$ p $。让我知道是否帮助,或许我们可以更新答案。
交错网格建模的主要困惑,你有多个网格(在这里,例如,压力$ p $粒子速度和水平和垂直v_x美元和v_z美元),和计算机实现只允许您在完整的整数索引p [1],p [2],vx [1],……,but in practice they may be offset, i.e., the location ofp [1]不一样的位置vx [1]实际上,他们所抵消\δx / 2美元,只是被相同的整数索引。
% % %设置空间和时间模型——真正的尺寸(单位:米)lx = 1200;ly = 1200;% - - - - - -真正的总运行时间(以秒为单位)lt = 5;% - - - - - -真正的介质参数(m / s和公斤/立方米)vp = 2000;ρ= 1500;% - - - - - -有限差分参数dx = 2;dy = dx;dt = 0.5 e - 3;% - - - - - -相关的离散网格在时间和空间上nx =地板(lx / dx);纽约=地板(ly / dy); nt = floor(lt/dt); x = [0:nx-1]*dx; y = [0:ny-1]*dy; t = [0:nt-1]*dt; % --- Source-time wavelet (supply central frequency fc, in Hertz) ricker = @(fm,t) (1-2*pi^2*fm^2*t.^2) .* exp(-pi^2*fm^2.*t.^2); fc = 20; fs = ricker(fc, t-0.1 ); % Time-delayed Ricker wavelet (plot, if desired!) % --- Source location (positioned at 400, 600) [~,sidx] = min(abs(x - 400)); [~,sidy] = min(abs(x - 600)); %% The computational loop % Initialize fields p = zeros(nx+1,ny+1); vx = zeros(nx+1,ny+1); vy = zeros(nx+1,ny+1); % Stencil update locations X = 1:nx; Y = 1:ny; figure(1); for a = 1:nt % (1) Inject source funtion p(sidx,sidy) = p(sidx,sidy) + dt*fs(a); % (2) Update p p(X+1,Y+1) = p(X+1,Y+1) - (dt*rho*vp^2)*((vx(X+1,Y) - vx(X,Y))/dx ... +(vy(X,Y+1) - vy(X,Y))/dy ); % (3a) Update ux vx(X,Y) = vx(X,Y) - (dt/rho)*(p(X+1,Y+1) - p(X,Y+1))/dx; % (3b) Update uy vy(X,Y) = vy(X,Y) - (dt/rho)*(p(X+1,Y+1) - p(X+1,Y))/dy; % Display every 10th computation if (mod(a,10)==0) imagesc(x,y,p'); axis equal; title(['Timestep = ', num2str(a), 'dt']); xlabel('x [m]'); ylabel('y [m]'); drawnow end end