数值模拟的代码问题求解

Question

对网格均匀划分，采用 WENO - 3（即每套模版包含 3 个网格）格式（可选有限体积或有限差分方法）结合空间上的 Lax - Friedrichs 格式，配合 2 阶的 Rung - Kutta 时间积分方法，编制代码（Fortran+matlab两个）计算上述问题在t = 2时的解

Answer 1

让【宝码香车】来帮你解答，本回答参考DeepSeek编写，并整理提供。

如果还有疑问可以点击头像关注私信或评论。如果答案让您满意，请采纳、点赞、关注，非常感谢！

我将为您提供解决无黏Burgers方程的数值模拟代码，采用WENO-3格式和Lax-Friedrichs通量分裂方法，配合二阶Runge-Kutta时间积分。

问题分析

1. 控制方程：无黏Burgers方程 ∂u/∂t + ∂(u²/2)/∂x = 0
2. 初始条件：Riemann问题 u(x,0) = 0.5 (x<2), 0.2 (x≥2)
3. 计算域：x∈[0,4]
4. 数值方法：
   • WENO-3重构（3点模板）
   • Lax-Friedrichs通量分裂
   • 二阶Runge-Kutta时间积分

Fortran 代码实现

program BurgersWENO
    implicit none
    integer, parameter :: Nx = 200, Nt = 1000
    real(8), parameter :: xmin = 0.0d0, xmax = 4.0d0
    real(8), parameter :: tfinal = 2.0d0
    real(8), parameter :: CFL = 0.5d0
    real(8) :: dx, dt, x(Nx), u(Nx), u_new(Nx), u1(Nx), u2(Nx)
    real(8) :: flux(Nx+1), t
    integer :: i, n
    
    ! 初始化网格和初始条件
    dx = (xmax - xmin) / Nx
    do i = 1, Nx
        x(i) = xmin + (i-0.5d0)*dx
        if (x(i) < 2.0d0) then
            u(i) = 0.5d0
        else
            u(i) = 0.2d0
        endif
    enddo
    
    ! 时间推进
    t = 0.0d0
    do n = 1, Nt
        dt = CFL * dx / maxval(abs(u))
        if (t + dt > tfinal) dt = tfinal - t
        
        ! 第一步 RK2
        call compute_flux(u, flux, Nx, dx, dt)
        do i = 2, Nx-1
            u1(i) = u(i) - dt/dx * (flux(i+1) - flux(i))
        enddo
        ! 边界条件
        u1(1) = u1(2)
        u1(Nx) = u1(Nx-1)
        
        ! 第二步 RK2
        call compute_flux(u1, flux, Nx, dx, dt)
        do i = 2, Nx-1
            u2(i) = 0.5d0*(u(i) + u1(i) - dt/dx * (flux(i+1) - flux(i)))
        enddo
        ! 边界条件
        u2(1) = u2(2)
        u2(Nx) = u2(Nx-1)
        
        u = u2
        t = t + dt
        if (t >= tfinal) exit
    enddo
    
    ! 输出结果
    open(unit=10, file='solution_fortran.dat')
    do i = 1, Nx
        write(10,*) x(i), u(i)
    enddo
    close(10)
    
contains


    subroutine compute_flux(u, flux, Nx, dx, dt)
        integer, intent(in) :: Nx
        real(8), intent(in) :: u(Nx), dx, dt
        real(8), intent(out) :: flux(Nx+1)
        real(8) :: alpha, f(Nx), uL, uR
        integer :: i
        
        ! Lax-Friedrichs通量分裂
        alpha = maxval(abs(u))
        do i = 1, Nx
            f(i) = 0.5d0 * u(i)**2
        enddo
        
        ! 计算界面通量 (WENO-3重构)
        do i = 2, Nx
            ! 左侧重构 u_{i-1/2}^+
            call WENO3(u(i-2:i), uL)
            ! 右侧重构 u_{i-1/2}^-
            call WENO3(u(i:i-2:-1), uR)
            flux(i) = 0.5d0*(f_LF(uL,alpha) + f_LF(uR,alpha))
        enddo
        
        ! 边界通量
        flux(1) = f(1)
        flux(Nx+1) = f(Nx)
    end subroutine
    
    subroutine WENO3(u_stencil, u_out)
        real(8), intent(in) :: u_stencil(3)
        real(8), intent(out) :: u_out
        real(8) :: beta1, beta2, w1, w2, eps
        real(8) :: p1, p2
        
        eps = 1.0d-6
        
        ! 计算光滑指示器
        beta1 = (u_stencil(2) - u_stencil(1))**2
        beta2 = (u_stencil(3) - u_stencil(2))**2
        
        ! 计算权重
        w1 = 1.0d0/(eps + beta1)**2
        w2 = 2.0d0/(eps + beta2)**2
        w1 = w1 / (w1 + w2)
        w2 = 1.0d0 - w1
        
        ! 重构多项式
        p1 = 1.5d0*u_stencil(2) - 0.5d0*u_stencil(1)
        p2 = 0.5d0*u_stencil(2) + 0.5d0*u_stencil(3)
        
        u_out = w1*p1 + w2*p2
    end subroutine
    
    function f_LF(u, alpha) result(f)
        real(8), intent(in) :: u, alpha
        real(8) :: f
        f = 0.5d0*(0.5d0*u**2 + alpha*u)
    end function
    
end program BurgersWENO

MATLAB 代码实现

function BurgersWENO_Matlab()
    % 参数设置
    Nx = 200;
    xmin = 0;
    xmax = 4;
    tfinal = 2;
    CFL = 0.5;
    
    % 初始化网格和初始条件
    dx = (xmax - xmin) / Nx;
    x = linspace(xmin+dx/2, xmax-dx/2, Nx);
    u = 0.5 * ones(size(x));
    u(x >= 2) = 0.2;
    
    % 时间推进
    t = 0;
    while t < tfinal
        dt = CFL * dx / max(abs(u));
        if t + dt > tfinal
            dt = tfinal - t;
        end
        
        % 第一步 RK2
        flux = compute_flux(u, dx, dt);
        u1 = u;
        u1(2:end-1) = u(2:end-1) - dt/dx * (flux(3:end) - flux(2:end-1));
        u1(1) = u1(2);
        u1(end) = u1(end-1);
        
        % 第二步 RK2
        flux = compute_flux(u1, dx, dt);
        u2 = u;
        u2(2:end-1) = 0.5*(u(2:end-1) + u1(2:end-1) - ...
                      dt/dx * (flux(3:end) - flux(2:end-1));
        u2(1) = u2(2);
        u2(end) = u2(end-1);
        
        u = u2;
        t = t + dt;
    end
    
    % 绘制结果
    plot(x, u, 'b-', 'LineWidth', 2);
    xlabel('x');
    ylabel('u');
    title(['无黏Burgers方程解 @ t = ' num2str(tfinal)]);
    grid on;
    saveas(gcf, 'solution_matlab.png');
end


function flux = compute_flux(u, dx, dt)
    Nx = length(u);
    alpha = max(abs(u));
    f = 0.5 * u.^2;
    flux = zeros(1, Nx+1);
    
    % WENO-3重构计算界面通量
    for i = 2:Nx
        % 左侧重构 u_{i-1/2}^+
        uL = WENO3(u(i-2:i));
        % 右侧重构 u_{i-1/2}^-
        uR = WENO3(u(i:-1:i-2));
        flux(i) = 0.5*(f_LF(uL,alpha) + f_LF(uR,alpha));
    end
    
    % 边界通量
    flux(1) = f(1);
    flux(end) = f(end);
end


function u_out = WENO3(u_stencil)
    eps = 1e-6;
    
    % 计算光滑指示器
    beta1 = (u_stencil(2) - u_stencil(1))^2;
    beta2 = (u_stencil(3) - u_stencil(2))^2;
    
    % 计算权重
    w1 = 1/(eps + beta1)^2;
    w2 = 2/(eps + beta2)^2;
    w1 = w1 / (w1 + w2);
    w2 = 1 - w1;
    
    % 重构多项式
    p1 = 1.5*u_stencil(2) - 0.5*u_stencil(1);
    p2 = 0.5*u_stencil(2) + 0.5*u_stencil(3);
    
    u_out = w1*p1 + w2*p2;
end


function f = f_LF(u, alpha)
    f = 0.5*(0.5*u^2 + alpha*u);
end

方法说明

关键点说明

1. WENO-3重构：使用3点模板进行非线性加权重构，确保在激波附近的高分辨率
2. Lax-Friedrichs通量分裂：将通量分裂为正向和反向分量，保证数值稳定性
3. 二阶Runge-Kutta：采用TVDRK2时间积分方法，保持时间方向二阶精度
4. CFL条件：根据最大特征速度动态调整时间步长，保证计算稳定性

计算结果将显示在t=2时刻的解，可以观察到激波的形成和传播。Fortran代码输出到文件solution_fortran.dat，MATLAB代码会生成图像solution_matlab.png。