当前位置: 首页 > 地学 >基于有限差分的波场模拟,pml边界,震源为单力源

基于有限差分的波场模拟,pml边界,震源为单力源

作者 Ego悄吟
来源: 小木虫 450 9 举报帖子
+关注

本人最近在做波场模拟的程序设计,在波场分离的时候遇到下图问题,震源附近发生突变,有没有大佬能够指出问题所在?
模型大小500x500x500,vp为3200km/s,vs为vp/1.732,密度为2250kg/m^3,震源位于中心,采用雷克子波,中心频率为50hz。

基于有限差分的波场模拟,pml边界,震源为单力源
水平.jpg


基于有限差分的波场模拟,pml边界,震源为单力源-1
3D.jpg 返回小木虫查看更多

今日热帖
  • 精华评论
  • diwushuxue

    网格多大,你这跟波场分离没关系,是你正演没做好。有可能是频散过于严重

  • Ego悄吟

    引用回帖:
    2楼: Originally posted by diwushuxue at 2019-04-17 11:39:28
    网格多大,你这跟波场分离没关系,是你正演没做好。有可能是频散过于严重

    您好,看到您的回帖之后,我回去划细网格了,之前网格大小为10m,计算步长为1.60e-3s,之后网格大小为2.5m,计算步长为0.8e-3s,可是依然出现这个情况。
    基于有限差分的波场模拟,pml边界,震源为单力源-2
    3D.jpg


    基于有限差分的波场模拟,pml边界,震源为单力源-3
    水平.jpg

  • diwushuxue

    引用回帖:
    3楼: Originally posted by Ego悄吟 at 2019-04-18 09:21:33
    您好,看到您的回帖之后,我回去划细网格了,之前网格大小为10m,计算步长为1.60e-3s,之后网格大小为2.5m,计算步长为0.8e-3s,可是依然出现这个情况。

    3D.jpg

    水平.jpg
    ...

    还是正演问题,你空间网格10m,时间采样1ms,主频30Hz。再试一下!

  • Ego悄吟

    引用回帖:
    4楼: Originally posted by diwushuxue at 2019-04-18 17:34:15
    还是正演问题,你空间网格10m,时间采样1ms,主频30Hz。再试一下!...

    您好,根据您给的建议,我回去试了试,结果还是出现之前的问题。
    基于有限差分的波场模拟,pml边界,震源为单力源-4
    3D.jpg


    基于有限差分的波场模拟,pml边界,震源为单力源-5
    水平.jpg

  • baobiao007

    明显是正演的递推过程没写对。。

  • baobiao007

    送你段例子代码吧,直接运行即可


    % 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 velocity
    nt = 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');,

猜你喜欢
下载小木虫APP
与700万科研达人随时交流
  • 二维码
  • IOS
  • 安卓