| 查看: 1706 | 回复: 9 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
基于有限差分的波场模拟，pml边界，震源为单力源已有2人参与
|
||
|
本人最近在做波场模拟的程序设计，在波场分离的时候遇到下图问题，震源附近发生突变，有没有大佬能够指出问题所在？ 模型大小500x500x500，vp为3200km/s，vs为vp/1.732，密度为2250kg/m^3，震源位于中心，采用雷克子波，中心频率为50hz。  水平.jpg  3D.jpg |
回复此楼» 本帖@通知
@在
» 猜你喜欢
报！！！某书最新消息说月底出结果
已经有13人回复
-
今年的陕西省自然科学基金啥时候出结果？
已经有10人回复
-
地球化学论文润色/翻译怎么收费?
已经有70人回复
-
有大佬投过Earth System Science Data吗
已经有0人回复
-
江苏科技大学2024年人才招聘公告
已经有0人回复
-
厦门工学院2024年学科带头人和专任教师招聘启事
已经有0人回复
-
Elsevier审稿
已经有11人回复
-
面上评委可能变为7个了，是好是坏呢？
已经有20人回复
baobiao007
木虫 (职业作家)
中国特色
- 应助: 201 (大学生)
- 金币: 6482.7
- 散金: 557
- 红花: 40
- 帖子: 3050
- 在线: 1009.9小时
- 虫号: 505962
- 注册: 2008-02-18
- 专业: 应用地球物理学
6楼2019-04-19 19:21:44
baobiao007
木虫 (职业作家)
中国特色
- 应助: 201 (大学生)
- 金币: 6482.7
- 散金: 557
- 红花: 40
- 帖子: 3050
- 在线: 1009.9小时
- 虫号: 505962
- 注册: 2008-02-18
- 专业: 应用地球物理学
【答案】应助回帖
★ ★ ★ ★ ★
感谢参与，应助指数 +1
Ego悄吟: 金币+5, ★★★★★最佳答案 2019-04-21 09:30:08
感谢参与，应助指数 +1
Ego悄吟: 金币+5, ★★★★★最佳答案 2019-04-21 09:30:08
|
送你段例子代码吧，直接运行即可 % Finite-differences in time domain(FDTD) acoustic wave propagation in 2D % medium close all; % Output every ... time steps IT_DISPLAY = 10; %% MODEL % Model dimensions nx = 401; nz = 401; dx = 10; % [m] dz = 10; % [m] % Elastic parameters vp = 3300.0 * ones(nz, nx); % velocity of compressional waves, [m/s] %% TIME STEPPING t_total = .55; % [sec] recording duration dt = 0.7 * min(dx,dz)/max(vp(  ); % min grid space / max velocitynt = round(t_total/dt); % number of time steps t = [0:nt]*dt; CFL = max(vp( ) * dt / min(dx,dz);%% SOURCE f0 = 10.0; % dominant frequency of the wavelet t0 = 1.20 / f0; % half Ricker wavelet excitation time factor = 1e10; % amplitude coefficient angle_force = 90.0; % spatial orientation of source % (not relevant for acoustic case) jsrc = round(nz/2); % source location along OZ isrc = round(nx/2); % source location along OX a = pi*pi*f0*f0; dt2 = dt^2; source_term = factor * exp(-a*(t-t0).^2); % Gaussian force_x = sin(angle_force * pi / 180) * source_term * dt2 / (dx * dz); min_wavelengh = 0.5*min(vp(vp>330))/f0; % shortest wavelength bounded by velocity in the air %% SUMMARY fprintf('#################################################\n'); fprintf('2D acoustic FDTD wave propagation in isotripic \nmedium in displacement formulation\n'); fprintf('#################################################\n'); fprintf('Model:\n\t%d x %d\tgrid nz x nx\n\t%.1e x %.1e\t[m] dz x dx\n',nz, nx, dz,dx); fprintf('\t%.1e x %.1e\t[m] model size\n',nx*dx, nz*dz); fprintf('\t%.1e...%.1e\t[m/s] vp\n', min(vp( ), max(vp());fprintf('Time:\n\t%.1e\t[sec] total\n\t%.1e\tdt\n\t%d\ttime steps\n',t_total,dt,nt); fprintf('Source:\n\t%.1e\t[Hz] dominant frequency\n\t%.1f\t[sec] index time\n',f0,t0); fprintf('Other:\n\t%.1f\tCFL number\n', CFL); fprintf('\t%.2f\t[m] shortest wavelength\n\t%d, %d\t points-per-wavelength OX, OZ\n', min_wavelengh, floor(min_wavelengh/dx), floor(min_wavelengh/dz)); fprintf('#################################################\n'); %% ALLOCATE MEMORY FOR WAVEFIELD p3 = zeros(nz+2,nx+2); % Wavefields at t p2 = zeros(nz+2,nx+2); % Wavefields at t-1 p1 = zeros(nz+2,nx+2); % Wavefields at t-2 % Coefficients for derivatives co_dxx = 1/dx^2; co_dzz = 1/dz^2; %% Loop over TIME tic; for it = 1:nt p3 = zeros(size(p2)); % Second-order derivatives dp_dxx = co_dxx * (p2(2:end-1,1:end-2) - 2*p2(2:end-1,2:end-1) + p2(2:end-1,3:end)); dp_dzz = co_dzz * (p2(1:end-2,2:end-1) - 2*p2(2:end-1,2:end-1) + p2(3:end,2:end-1)); % U(t) = 2*U(t-1) - U(t-2) + G dt2/rho; p3(2:end-1,2:end-1) = 2.0*p2(2:end-1,2:end-1) - p1(2:end-1,2:end-1) + (vp.^2).*(dp_dxx + dp_dzz).*dt2; % Add source term p3(jsrc, isrc) = p3(jsrc, isrc) + force_x(it); % Exchange data between t-2 (1), t-1 (2) and t (3) and apply ABS p1 = p2 .* weights; p2 = p3 .* weights; % Output if mod(it,IT_DISPLAY) == 0 fprintf('Time step: %d \t %.4f s\n',it, single(t(it))); imagesc(p3); colorbar; axis equal tight; colormap jet; title(['Step = ',num2str(it),'/',num2str(nt),', Time: ',sprintf('%.4f',t(it)),' sec']); drawnow; end end toc; disp('End'); |
8楼2019-04-20 10:36:59
baobiao007
木虫 (职业作家)
中国特色
- 应助: 201 (大学生)
- 金币: 6482.7
- 散金: 557
- 红花: 40
- 帖子: 3050
- 在线: 1009.9小时
- 虫号: 505962
- 注册: 2008-02-18
- 专业: 应用地球物理学
9楼2019-04-20 10:38:46
80