精华内容
下载资源
问答
  • 解二维抛物型方程初边值问题的交替方向隐式方法.详细地研究了各种情形的求解二维抛物型方程初边值问题的交替方向隐式差分法。这 种方法可以将二维隐式方法归结为求解三对角线性方程组。与一维情况类似 ,可以继续...
  • 提出了数值求解三维波动方程的2种精度分别为0(τ2 +h4)和0(τ4+h4)的交替方向隐式(ADI)格式,并且通过Fourier方法证明了格式的稳定性。该方法在沿每个空间方向上只涉及三个网格基架点,因此可以重复采用TDMA算法,从而...
  • 交替方向隐式欧拉方法研究二维带有弱奇异核的偏积分微分方程的数值解,在空间方向上采用二阶差商,时间方向上使用向后欧拉方法,积分项用一阶卷积求积逼近,该方法具备了交替方向存储量少,计算量低的特点.
  • ADI交替方向隐式法求抛物型PDE的fortran程序
  • 在本文中,我们使用Adomian分解法(ADM),有限差分法和交替方向隐式法来评估上述方法的优缺点。 为此,我们对使用这些方法构造的不同解决方案进行了数值模拟,并比较了所研究误差的情况。
  • 在该方法中,表征器件全矢量特性的偏振耦合项被引入光束传播方程,并使用新型的交替方向隐式进行高效处理。数值模拟的结果表明,相比于传统的迭代方法,该方法在保持较高精度的同时,计算效率有了显著的提高。
  • 此函数使用交替方向隐式 (ADI) 方法求解均匀介质中的三维 Pennes 生物热传递 (BHT) 方程。 该代码是为组织中的高强度聚焦超声 (HIFU) 治疗而开发的,但它也可以应用于其他加热问题。 如果需要,该解决方案会考虑组织...
  • 第2节讨论了热力学模拟物理。第3节概述了热传导的数值公式。第4节介绍了通过ADI方法进行的热模拟。实施和实验结果在第5节中提出
  • 基于单轴各向异性完全匹配层边界条件的交替方向隐式有限差分法GPR全波场数值模拟,冯德山,陈承申,交替方向隐式差分(ADI-FDTD)法突破了Courand-Friedrich-Levy(CFL)条件的约束,具有无条件稳定的特点;而单轴各向...
  • 二维ADI(交替方向隐式求解)

    千次阅读 2018-09-22 09:40:07
    数值解就是在我关注的点上,通过差分计算方法得到的值。 精确解是带入关注的点到表达式。 u(i,j)是矩阵 x=0:h:1是向量。 x(i)便可以取出来一个值。 关于时间和整个循环。先执行内部的循环,再...

    1问题描述

    用ADI法求解二维抛物线问题

    2.物理模型:

    3.数学描述:

    4.算法描述:

    1. 对于时间层n,增加虚拟时间层n+1\2;
    2. 对uxx在时间层n+1\2上进行前向差分,对uyy在时间层上n上进行后向差分;
    3. 对uxx在时间层n+1\2上进行向后差分,对uyy在时间层n+1上进行向前差分;
    4. 已知u在时间层n上的值,按行求解u在事件层n+1\2层的值;
    5. 将u在n+1\2上的值代入到第二个方程中,形成b,按列求解u在时间层n+1上的差分解,即u(i,j,k)。

    5.程序伪代码

             (1)区域离散;

               x=i*h,y=j*h,t=tao*ht。或者x=0:h:1,y=0:h:1,t=0:ht:t

              注意在编程的时候,一定要把0:m写成1:m+1,就是每个坐标自动加1。在纸上写的时候,和在电脑上写的时候是不一样的。         

           (2)初始值和边界值离散;

               u(x,y,t)=u(i*h,j*h,tao*ht);

               其实就是固定一个值,然后让i,j进行变化,也叫扫描过程。

               并且给出初始的值。u(x,y,n)的值是给定的,存起来然后调用。

             (3)方程离散;

                差分法离散成两个方程,写成Ax=b的形式。

                并且给出u的第一行和最后一行的值,按照边界条件给出,给出b的第一行和最后一行的值。

             (4)迭代求解;

                在时间层n,

               已知u(i,j,n),求解u(i,j,n+1\2);

                固定列数,按行进行扫描,对每一列都有一个方程组,Ax=b;求解后按列存入u(i,j,n+1\2)矩阵中。

                固定行数,按列进行扫描,对每一行都有一个方程组,Ax=b;求解后按行存入u(i,j,n+1)矩阵中。

             (5)后处理和作图(可视化过程)

    plot(x,y,t,u(i,j,k)

    surf(x,y,t,u(i,j,k))

    编程思路:

    首先离散空间区域和时间区域,

    然后固定时间层,再固定列数,对行进行扫描,

    求解一个Ax=b的方程组,

    先把b写出来,然后把A写出来,然后计算x=b\A;

    再将x1的值送到b2中,解出来x2,装入U中。

    数值解就是在我关注的点上,通过差分计算方法得到的值。

    精确解是带入关注的点到表达式。

    u(i,j)是矩阵

    x=0:h:1是向量。

    x(i)便可以取出来一个值。

    关于时间和整个循环。先执行内部的循环,再执行第二个循环,然后是时间的循环,解下一个方程,最后需要得到u2的.

    u2只不过是一个中间值。不需要知道u2在边界上的值的。

    值随着时间变化,并且用新的u1的值来更新b的值,从而求得新的u1的解。

    需要知道u1在边界上的值。

    这个是完善u1这个表格。

    展开全文
  • 煤层气开采过程中通常采用...交替方向格式可以把二维问题转化成一维问题,对x,y两个方向的迭代矩阵均为三对角矩阵,结构相同,易于编程计算。所建立的差分格式具有计算量小,稳定性好等优点,数值试验的结果表明效果良好。
  • 针对变系数空间分数阶电报方程,利用Grtinwald-Letnikov分数阶导数的定义,在交替方向法的基础上构造了一种修正交替方向隐式差分格式.通过Fourier分析和Lax等价定理证明了所提出的格式是绝对稳定、相容和无条件收敛...
  • 研究了三维常系数反应扩散方程的紧交替方向隐式差分格式.首先综合运用降阶法和降维法导出了紧差分格式,并给出了差分格式截断误差的表达式;其次引进过渡层变量,给出了紧交替方向隐式差分格式算法;利用Fouier稳定性...
  • 本文介绍了一种旨在减少3-D高阶交替方向隐式有限差分时域(ADI-FDTD)方法数值离散的新方法。 首先,我们用人工各向异性介电材料修改了3-D高阶ADI-FDTD方法的数值公式,并通过分析得出了新的数值色散关系。 另外,给...
  • 研究了二维变系数非齐次抛物型方程的紧交替方向隐式差分格式,首先运用算子方法导出了紧差分格式,给出了差分格式的截断误差,接着讨论了差分格式的稳定性和收敛性,最后给出了数值例子,数值结果和理论分析是吻合的。
  • 三维导热问题的ADI-TDMA算法

    千次阅读 2020-04-30 17:06:03
    本文在上一篇文章中关于二维非稳态导热问题的交替方向隐式方法的基础上,设计了一个三维算例,并通过Fortran程序,演示适用于三维模型的Brian ADI格式的交替方向隐式方法的求解过程。

    三维导热问题的ADI-TDMA算法

    基本原理

    上一篇文章,讲述了二维导热问题采用交替方向隐式方法[1]的求解过程。对于二维问题,Peaceman-Rachford ADI格式是无条件稳定的。然而,该格式对于三维问题是条件稳定的,其稳定条件为:

    aΔt[1Δx2+1Δy2+1Δz2]1.5 a\Delta t[\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}+\frac{1}{\Delta z^2}]\le 1.5

    可见该格式对时间步长Δt\Delta t的限制很严格。相比之下,Brian ADI格式则是绝对稳定的,故可以选取较大的时间步长,同时每一时间步的计算量较小,可大幅加快计算速度。

    下面介绍Brian ADI格式下三维非稳态导热方程的求解过程,求解过程中涉及两个温度中间量:UUVV

    基于有限体积法的全隐离散方程可表示为:
    aP0(TPTP0)=aE(TETP)aW(TPTW)+aF(TFTP)aB(TPTB)+aN(TNTP)aS(TPTS)+SΔxΔyΔz a_P^0(T_P-T_P^0)=a_E(T_E-T_P)-a_W(T_P-T_W)+a_F(T_F-T_P)-a_B(T_P-T_B)+a_N(T_N-T_P)-a_S(T_P-T_S)+S\Delta x\Delta y\Delta z

    (1)沿xx方向求解:

    先取Δt/2\Delta t/2的时间步长,对x方向做隐式处理,yz方向取温度场T0T^0的值,求解温度场的中间量UU,则全隐离散方程可改成:

    aP0(UPTP0)=aE(UEUP)aW(UPUW)+aF(TF0TP0)aB(TP0TB0)+aN(TN0TP0)aS(TP0TS0)+SΔxΔyΔz a_P^0(U_P-T_P^0)=a_E(U_E-U_P)-a_W(U_P-U_W)+a_F(T_F^0-T_P^0)-a_B(T_P^0-T_B^0)+a_N(T_N^0-T_P^0)-a_S(T_P^0-T_S^0)+S\Delta x\Delta y\Delta z

    整理后得到:

    (aP0+aE+aW)UP=aEUE+aWUW+(aP0aFaBaNaS)TP0+aFTF0+aBTB0+aNTN0+aSTS0+SΔxΔyΔz (a_P^0+a_E+a_W)U_P=a_EU_E+a_WU_W+(a_P^0-a_F-a_B-a_N-a_S)T_P^0+a_FT_F^0+a_BT_B^0+a_NT_N^0+a_ST_S^0+S\Delta x\Delta y\Delta z

    在上式中,除了UPU_PUEU_EUWU_W,其他量都是已知的。令aP=aP0+aE+aWa_P=a_P^0+a_E+a_WD=(aP0aFaBaNaS)TP0+aFTF0+aBTB0+aNTN0+aSTS0+SΔxΔyΔzD=(a_P^0-a_F-a_B-a_N-a_S)T_P^0+a_FT_F^0+a_BT_B^0+a_NT_N^0+a_ST_S^0+S\Delta x\Delta y\Delta z,则上式可以简化为:

    aPUP=aEUE+aWUW+D a_PU_P=a_EU_E+a_WU_W+D

    采用沿xx方向的TDMA可以求得温度场中间量UU

    (2)沿yy方向求解:

    yy方向做隐式处理,xx方向取求得的中间量UU的值,zz方向取温度场T0T^0的值,求解温度场的第二个中间量VV,则全隐离散方程可改成:

    aP0(VPTP0)=aE(UEUP)aW(UPUW)+aF(VFVP)aB(VPVB)+aN(TN0TP0)aS(TP0TS0)+SΔxΔyΔz a_P^0(V_P-T_P^0)=a_E(U_E-U_P)-a_W(U_P-U_W)+a_F(V_F-V_P)-a_B(V_P-V_B)+a_N(T_N^0-T_P^0)-a_S(T_P^0-T_S^0)+S\Delta x\Delta y\Delta z

    整理后得到:

    (aP0+aF+aB)VP=aFVF+aBVB+(aP0aNaS)TP0+aNTN0+aSTS0(aE+aW)UP+aEUE+aWUW+SΔxΔyΔz (a_P^0+a_F+a_B)V_P=a_FV_F+a_BV_B+(a_P^0-a_N-a_S)T_P^0+a_NT_N^0+a_ST_S^0-(a_E+a_W)U_P+a_EU_E+a_WU_W+S\Delta x\Delta y\Delta z

    在上式中,除了VPV_PVFV_FVBV_B,其他量都是已知的。令aP=aP0+aF+aBa_P=a_P^0+a_F+a_BD=(aP0aNaS)TP0+aNTN0+aSTS0(aE+aW)UP+aEUE+aWUW+SΔxΔyΔzD=(a_P^0-a_N-a_S)T_P^0+a_NT_N^0+a_ST_S^0-(a_E+a_W)U_P+a_EU_E+a_WU_W+S\Delta x\Delta y\Delta z,则上式可以简化为:

    aPVP=aFVF+aBVB+D a_PV_P=a_FV_F+a_BV_B+D

    采用沿yy方向的TDMA可以求得温度场中间量VV

    (3)沿zz方向求解:

    zz方向做隐式处理,xx方向取仍取步骤(1)中求得的中间量UU的值,yy方向取步骤(2)求得的中间量VV的值,求解下一时刻温度场T1T^1,则全隐离散方程可改成:

    aP0(TP1VP)=aE(UEUP)aW(UPUW)+aF(VFVP)aB(VPVB)+aN(TN1TP1)aS(TP1TS1)+SΔxΔyΔz a_P^0(T_P^1-V_P)=a_E(U_E-U_P)-a_W(U_P-U_W)+a_F(V_F-V_P)-a_B(V_P-V_B)+a_N(T_N^1-T_P^1)-a_S(T_P^1-T_S^1)+S\Delta x\Delta y\Delta z

    整理后得到:

    (aP0+aN+aS)TP1=aNTN1+aSTS1+(aP0aFaB)VP0+aFVF+aBVB(aE+aW)UP+aEUE+aWUW+SΔxΔyΔz (a_P^0+a_N+a_S)T_P^1=a_NT_N^1+a_ST_S^1+(a_P^0-a_F-a_B)V_P^0+a_FV_F+a_BV_B-(a_E+a_W)U_P+a_EU_E+a_WU_W+S\Delta x\Delta y\Delta z

    在上式中,除了TP1T_P^1TF1T_F^1TB1T_B^1,其他量都是已知的。令aP=aP0+aN+aSa_P=a_P^0+a_N+a_SD=(aP0aFaB)VP0+aFVF+aBVB(aE+aW)UP+aEUE+aWUW+SΔxΔyΔzD=(a_P^0-a_F-a_B)V_P^0+a_FV_F+a_BV_B-(a_E+a_W)U_P+a_EU_E+a_WU_W+S\Delta x\Delta y\Delta z,则上式可以简化为:

    aPTP1=aNTN1+aSTS1+D a_PT_P^1=a_NT_N^1+a_ST_S^1+D

    采用沿zz方向的TDMA可以求得下一时刻的温度场T1T^1

    边界节点的处理方式与上述过程相似,主要区别在于该节点的离散方程中未知项少了一项,例如对于顶部对流边界,aNa_N的取值得改为αΔxΔy\alpha\Delta x\Delta yTNT_N由未知的待求解量改为已知量TfT_f

    求解程序

    下面设计一个算例,并通过Fortran程序给出其求解过程:

    有一个长600mm,宽400mm,高200mm的矩形计算域。计算域内材料的密度为2500kg/m3,热容为1372J/(kg·K)。含内热源,热生成率恒定为1050J/(m3·s)。导热系数为温度的多项式函数:

    λ=1.13588×1015T5+3.25358×1012T43.25305×109T3+1.32926×106T29.27637×105T+1.04478 \lambda=-1.13588\times 10^{-15}T^5+3.25358\times 10^{-12}T^4-3.25305\times 10^{-9}T^3+1.32926\times 10^{-6}T^2-9.27637\times 10^{-5}T+1.04478

    计算域的初始温度为1500℃,外表面与环境进行对流换热,其中顶部对流给热系数α1\alpha_1为200W/(m2·K),底部、左右侧和前后侧的对流给热系数α2\alpha_2为80W/(m2·K),环境温度TfT_f为20℃。模拟的物理时间长度为3h。

    采用均匀结构化网格,每个网格节点的尺寸为2mm×2mm×2mm。由于该算例的模型沿yOzyOz面、zOxzOx是对称的,实际计算域可以只取左后半部分,即选取的计算域尺寸为300mm×200mm×200mm,网格节点数量为151×101×101。计算域右侧面、前侧面作为对称边界处理,右侧面上TE=TWT_E=T_WaE=aWa_E=a_W;前侧面上TF=TBT_F=T_BaF=aBa_F=a_B

    在选取时间步长时,由于非稳态导热随着时间的进行,会趋向于稳态,刚开始时温度变化率会很大,随后则会大幅减小,本算例选取时间步长的准则是使选定的时间步长内最大温度变化幅度不超过10℃,这样刚开始时间步长会很小,随后则会变得很大。

    计算过程中选取了8个采样点,并设定每隔5分钟输出一次这8个采样点的温度值。

    选用的坐标系如下:
    三维坐标系
    求解程序如下,模拟用到的参数在网格材料边界条件初始条件求解控制模块中通过常量定义出来,以便于参数的调整:

    module mesh
        implicit none
        real, parameter :: length = 300e-3, width = 200e-3, height = 200e-3
        real, parameter :: delta_X = 2e-3, delta_Y = 2e-3, delta_Z = 2e-3
    end module mesh
    
    module material
        implicit none
        real, parameter :: rho = 2500, Cp = 1372, S = 1050
    end module material
    
    module boundary_condition
        implicit none
        real, parameter :: alpha1 = 300, alpha2 = 80
        real, parameter :: Tf = 20
    end module boundary_condition
    
    module initial_condition
        implicit none
        real, parameter :: T0 = 1500
    end module initial_condition
    
    module solution_control
        implicit none
        real, parameter :: time_length = 180 * 60
        real, parameter :: sampling_interval = 5 * 60
        real, parameter :: max_deltaT_per_step = 10
        real, dimension(8) :: sample_points_X = (/40e-3, 300e-3, 40e-3, 300e-3, 40e-3, 300e-3, 40e-3, 300e-3/)
        real, dimension(8) :: sample_points_Y = (/100e-3, 100e-3, 170e-3, 170e-3, 100e-3, 100e-3, 170e-3, 170e-3/)
        real, dimension(8) :: sample_points_Z = (/100e-3, 100e-3, 100e-3, 100e-3, 170e-3, 170e-3, 170e-3, 170e-3/)
        integer :: spIndex(8,3)
    end module solution_control
    
    
    program TDMA
        use mesh
        use material
        use solution_control
        implicit none
        integer :: Nx, Ny, Nz
        
        !确定X、Y方向的网格节点数量
        Nx = length / delta_X + 1
        Ny = width / delta_Y + 1
        Nz = height / delta_Z + 1
    
        !根据采样点位置获得采样点在网格节点中的索引
        spIndex(:,1) = floor(sample_points_X / delta_X + 1)
        spIndex(:,2) = floor(sample_points_Y / delta_Y + 1)
        spIndex(:,3) = floor(sample_points_Z / delta_Z + 1)
    
        !开始计算
        call calculate(Nx, Ny, Nz)
    
    end program TDMA
    
    
    subroutine calculate(Nx, Ny, Nz)
        use mesh
        use material
        use boundary_condition
        use initial_condition
        use solution_control
        implicit none
        integer, intent(in) :: Nx, Ny, Nz
        integer :: i
        real :: current_time, current_sampling_time, delta_Tau
        real, dimension(0:Nx+1,0:Ny+1,0:Nz+1) :: T, U, V, T1, delta_T
        real, dimension(Nx,Ny,Nz) :: lambda, aP0, aE, aW, aF, aB, aN, aS, aP, Sv, D, Dx, P, Q
        real :: kx(Ny,Nz), ky(Nx,Nz),kz(Nx,Ny)
        real :: max_delta_T, max_rate
        real :: tStart, tEnd
    
        !初始化系数矩阵A
        aP0 = 0
        aE = 0
        aW = 0
        aF = 0
        aB = 0
        aN = 0
        aS = 0
    
        !初始化温度矩阵
        T = Tf
        T(1:Nx,1:Ny,1:Nz) = T0
        U = Tf
        V = Tf
        T1 = Tf
    
        !初始化源项矩阵
        Sv = S * delta_X * delta_Y * delta_Z
        Sv(1,:,:) = Sv(1,:,:) / 2
        Sv(:,1,:) = Sv(:,1,:) / 2
        Sv(:,:,1) = Sv(:,:,1) / 2
        Sv(:,:,Ny) = Sv(:,:,Ny) / 2
    
        !初始化输出
        write(*,"(A4,$)") "Time"
        do i = 1, 7
            write(*,"(A9,I1,$)") "SP", i
        end do
        write(*,"(A10)") "SP8"
    
        !按时间步进行温度场计算
        current_time = 0
        current_sampling_time = sampling_interval
        delta_Tau = 0.0001  !考虑到降温刚开始时的温度变化率很大,初始delta_Tau取一个很小的值
        call CPU_TIME(tStart)
        do while (current_time < time_length)
            !更新当前时间
            current_time = current_time + delta_Tau
    
            !通过多项式函数计算材料导热系数
            lambda = -1.13588E-15 * T(1:Nx,1:Ny,1:Nz)**5 + 3.25358E-12 * T(1:Nx,1:Ny,1:Nz)**4 &
                   - 3.25305E-09 * T(1:Nx,1:Ny,1:Nz)**3 + 1.32926E-06 * T(1:Nx,1:Ny,1:Nz)**2 &
                   - 9.27637E-05 * T(1:Nx,1:Ny,1:Nz) + 1.04478
    
            !计算系数矩阵A
            !(1)内部节点的计算
            aP0 = rho * Cp * delta_X * delta_Y * delta_Z / (delta_Tau / 2)
            aE(1:Nx-1,:,:) = 2 / (1 / lambda(1:Nx-1,:,:) + 1 / lambda(2:Nx,:,:)) / delta_X * delta_Y * delta_Z
            aW(2:Nx,:,:) = 2 / (1 / lambda(2:Nx,:,:) + 1 / lambda(1:Nx-1,:,:)) / delta_X * delta_Y * delta_Z
            aF(:,1:Ny-1,:) = 2 / (1 / lambda(:,1:Ny-1,:) + 1 / lambda(:,2:Ny,:)) / delta_Y * delta_X * delta_Z
            aB(:,2:Ny,:) = 2 / (1 / lambda(:,2:Ny,:) + 1 / lambda(:,1:Ny-1,:)) / delta_Y * delta_X * delta_Z
            aN(:,:,1:Nz-1) = 2 / (1 / lambda(:,:,1:Nz-1) + 1 / lambda(:,:,2:Nz)) / delta_Z * delta_X * delta_Y
            aS(:,:,2:Nz) = 2 / (1 / lambda(:,:,2:Nz) + 1 / lambda(:,:,1:Nz-1)) / delta_Z * delta_X * delta_Y
            !(2)左侧、上下侧、后侧对流边界节点的处理
            aW(1,:,:) = alpha2 * delta_Y * delta_Z
            aN(:,:,Nz) = alpha1 * delta_X * delta_Y
            aS(:,:,1) = alpha2 * delta_X * delta_Y
            aB(:,1,:) = alpha2 * delta_X * delta_Z
            !(3)右侧、前侧对称边界节点的处理
            aW(Nx,:,:) = aW(Nx,:,:) * 2
            aF(:,Ny,:) = aF(:,Ny,:) * 2
            !(4)边界节点的减半处理
            !(4-1)上侧aP0、aE、aW、aF、aB减半
            aP0(:,:,Nz) = aP0(:,:,Nz) / 2
            aE(:,:,Nz) = aE(:,:,Nz) / 2
            aW(:,:,Nz) = aW(:,:,Nz) / 2
            aF(:,:,Nz) = aF(:,:,Nz) / 2
            aB(:,:,Nz) = aB(:,:,Nz) / 2
            !(4-2)下侧aP0、aE、aW、aF、aB减半
            aP0(:,:,1) = aP0(:,:,1) / 2
            aE(:,:,1) = aE(:,:,1) / 2
            aW(:,:,1) = aW(:,:,1) / 2
            aF(:,:,1) = aF(:,:,1) / 2
            aB(:,:,1) = aB(:,:,1) / 2
            !(4-3)左侧aP0、aN、aS、aF、aB减半
            aP0(1,:,:) = aP0(1,:,:) / 2
            aN(1,:,:) = aN(1,:,:) / 2
            aS(1,:,:) = aS(1,:,:) / 2
            aF(1,:,:) = aF(1,:,:) / 2
            aB(1,:,:) = aB(1,:,:) / 2
            !(4-4)后侧aE、aW、aN、aS减半
            aP0(:,1,:) = aP0(:,1,:) / 2
            aE(:,1,:) = aE(:,1,:) / 2
            aW(:,1,:) = aW(:,1,:) / 2
            aN(:,1,:) = aN(:,1,:) / 2
            aS(:,1,:) = aS(:,1,:) / 2
    
            !沿x方向使用TDMA隐式求解
            aP = aP0 + aE + aW
            !(1)定义TDMA中的常数矩阵
            D = (aP0 - aF - aB - aN - aS) * T(1:Nx,1:Ny,1:Nz)+ aF * T(1:Nx,2:Ny+1,1:Nz) + aB * T(1:Nx,0:Ny-1,1:Nz) &
              + aN * T(1:Nx,1:Ny,2:Nz+1) + aS * T(1:Nx,1:Ny,0:Nz-1) + Sv
            !(2)左侧边界的特殊处理
            D(1,:,:) = D(1,:,:) + aW(1,:,:) * Tf
            !!!!!!
            P(1,:,:) = aE(1,:,:) / aP(1,:,:)
            Q(1,:,:) = D(1,:,:) / aP(1,:,:)
            do i = 2, Nx
                kx = aP(i,:,:) - aW(i,:,:) * P(i-1,:,:)
                P(i,:,:) = aE(i,:,:) / kx
                Q(i,:,:) = (D(i,:,:) + aW(i,:,:) * Q(i-1,:,:)) / kx
            end do
            U(Nx,1:Ny,1:Nz) = Q(Nx,:,:)
            do i = Nx-1, 1, -1
                U(i,1:Ny,1:Nz) = P(i,:,:) * U(i+1,1:Ny,1:Nz) + Q(i,:,:)
            end do
            !(3)由于沿y、z方向定义常数矩阵时,均会用到x方向的计算量U加和,将其定义为Dx
            Dx = -(aE + aW) * U(1:Nx,1:Ny,1:Nz) + aE * U(2:Nx+1,1:Ny,1:Nz) + aW * U(0:Nx-1,1:Ny,1:Nz)
    
            !沿y方向使用TDMA隐式求解
            aP = aP0 + aF + aB
            !(1)定义TDMA中的常数矩阵
            D = (aP0 - aN - aS) * T(1:Nx,1:Ny,1:Nz) + aN * T(1:Nx,1:Ny,2:Nz+1) + aS * T(1:Nx,1:Ny,0:Nz-1) + Dx + Sv
            !(2)后侧边界的特殊处理
            D(:,1,:) = D(:,1,:) + aB(:,1,:) * Tf
            !!!!!!
            P(:,1,:) = aF(:,1,:) / aP(:,1,:)
            Q(:,1,:) = D(:,1,:) / aP(:,1,:)
            do i = 2, Ny
                ky = aP(:,i,:) - aB(:,i,:) * P(:,i-1,:)
                P(:,i,:) = aF(:,i,:) / ky
                Q(:,i,:) = (D(:,i,:) + aB(:,i,:) * Q(:,i-1,:)) / ky
            end do
            V(1:Nx,Ny,1:Nz) = Q(:,Ny,:)
            do i = Ny-1, 1, -1
                V(1:Nx,i,1:Nz) = P(:,i,:) * V(1:Nx,i+1,1:Nz) + Q(:,i,:)
            end do
    
            !沿z方向使用TDMA隐式求解
            aP = aP0 + aN + aS
            !(1)定义TDMA中的常数矩阵
            D = (aP0 - aF - aB) * V(1:Nx,1:Ny,1:Nz) + aF * V(1:Nx,2:Ny+1,1:Nz) + aB * V(1:Nx,0:Ny-1,1:Nz) + Dx + Sv
            !(2)上下侧边界的特殊处理
            D(:,:,1) = D(:,:,1) + aS(:,:,1) * Tf
            D(:,:,Nz) = D(:,:,Nz) + aN(:,:,Nz) * Tf
            !!!!!!
            P(:,:,1) = aN(:,:,1) / aP(:,:,1)
            Q(:,:,1) = D(:,:,1) / aP(:,:,1)
            do i = 2, Nz
                kz = aP(:,:,i) - aS(:,:,i) * P(:,:,i-1)
                P(:,:,i) = aN(:,:,i) / kz
                Q(:,:,i) = (D(:,:,i) + aS(:,:,i) * Q(:,:,i-1)) / kz
            end do
            T1(1:Nx,1:Ny,Nz) = Q(:,:,Nz)
            do i = Nz-1, 1, -1
                T1(1:Nx,1:Ny,i) = P(:,:,i) * T1(1:Nx,1:Ny,i+1) + Q(:,:,i)
            end do
    
            !当达到采样时间时,输出采样点温度
            if (current_time >= current_sampling_time) then
                !输出6个采样点的温度值
                write(*,"(I4,$)") floor(current_time / 60)
                do i = 1, 7
                    write(*,"(f10.2,$)") T1(spIndex(i,1), spIndex(i,2), spIndex(i,3))
                end do
                write(*,"(f10.2)") T1(spIndex(6,1), spIndex(6,2), spIndex(6,3))
                !
                !更新下一次需要采样的时间
                current_sampling_time = current_sampling_time + sampling_interval
            end if
    
            !计算最大冷却速率
            delta_T = abs(T1 - T)
            max_delta_T = maxval(delta_T)
            max_rate = max_delta_T / delta_Tau
    
            !根据对每个时间步内的最大温度变化率限制,以及下一次采样时间确定下一个时间步长
            delta_Tau = min(max_deltaT_per_step / max_rate, current_sampling_time - current_time)
    
            !更新温度矩阵T
            T = T1
    
        end do
        call CPU_TIME(tEnd)
        write(*,"('Calculate complete, time consumption: ', f8.4, ' sec')") (tEnd - tStart)
    
    end subroutine calculate
    

    参考文献

    [1] 陶文铨. 数值传热学[M]. 第二版. 西安: 西安交通大学出版社, 2001: 99-104.

    展开全文
  • 简单的二维抛物线方程例子 采用二维交替方向隐格式求解 并且附有matlab程序 适合借鉴
  • 该文基于SMAC(simplified marker and cell...压力Poisson方程用Tschebyscheff SLOR方法交替方向迭代求解。Reynolds时均动量方程、k方程和ε方程对流项均采用Chakravarthy-Osher TVD格式离散,该格式不但有助于提高数值
  • 二维导热问题的ADI-TDMA算法

    千次阅读 2020-04-29 12:50:15
    交替方向隐式方法是CFD中用于求解对流-扩散问题的一种非常高效的求解方法,它将二维/三维隐式求解问题转化为各方向的一维隐式求解问题处理,并采用TDMA求解,比直接求解多维问题的线性方程组更高效,且时间步长选取...

    二维导热问题的ADI-TDMA算法

    基本原理

    对于多维非稳态导热问题,其时间项的离散格式有两种:显式格式隐式格式。在显式格式中,下一时间步的温度场由上一时间步的已知量求出,不需要求解代数方程组,但其时间步的选取受稳定性条件的限制,故通常只能取很小的时间步。隐式格式不受稳定性限制,但其需要对五对角矩阵(二维)或七对角矩阵(三维)进行求解,故每一时间步的计算量较大。

    相比之下,交替方向隐式方法[1](alternating direction implicit, ADI)是在显式与隐式之间折中的一种方法。对于二维问题,先将xx方向的温度值按隐式处理,yy方向的温度值按显式处理,这样就将二维问题转化为一系列并列的一维问题,可以直接使用三对角阵算法(TDMA)进行计算;然后再倒过来,将yy方向隐式,xx方向显式处理。这种方法被称为Peaceman-Rachford ADI格式。可以证明,该方法对于二维问题是绝对稳定的,故没有稳定性限制,可以取较大的时间步长,同时每一时间步的计算量较小,可大幅加快计算速度。

    下面介绍Peaceman-Rachford ADI格式下二维非稳态导热方程的求解过程。

    基于有限体积法的全隐离散方程可表示为:
    aP0(TPTP0)=aE(TETP)aW(TPTW)+aN(TNTP)aS(TPTS)+SΔxΔy a_P^0(T_P-T_P^0)=a_E(T_E-T_P)-a_W(T_P-T_W)+a_N(T_N-T_P)-a_S(T_P-T_S)+S\Delta x\Delta y

    (1)沿xx方向求解:

    先取Δt/2\Delta t/2的时间步长,对x方向做隐式处理,求解温度场的中间量UU,则全隐离散方程可改成:
    aP0(UPTP0)=aE(UEUP)aW(UPUW)+aN(TN0TP0)aS(TP0TS0)+SΔxΔy a_P^0(U_P-T_P^0)=a_E(U_E-U_P)-a_W(U_P-U_W)+a_N(T_N^0-T_P^0)-a_S(T_P^0-T_S^0)+S\Delta x\Delta y

    整理后得到:

    (aP0+aE+aW)UP=aEUE+aWUW+(aP0aNaS)TP0+aNTN0+aSTS0+SΔxΔy (a_P^0+a_E+a_W)U_P=a_EU_E+a_WU_W+(a_P^0-a_N-a_S)T_P^0+a_NT_N^0+a_ST_S^0+S\Delta x\Delta y

    在上式中,除了UPU_PUEU_EUWU_W,其他量都是已知的。令aP=aP0+aE+aWa_P=a_P^0+a_E+a_WD=(aP0aNaS)TP0+aNTN0+aSTS0+SΔxΔyD=(a_P^0-a_N-a_S)T_P^0+a_NT_N^0+a_ST_S^0+S\Delta x\Delta y,则上式可以简化为:

    aPUP=aEUE+aWUW+D a_PU_P=a_EU_E+a_WU_W+D

    通过使用沿xx方向的TDMA可以求得温度场中间量UU

    (2)沿yy方向求解:

    在求得中间量UU的基础上,再取Δt/2\Delta t/2的时间步长,对y方向做隐式处理,求解下一时刻温度场T1T^1,则全隐离散方程可改成:
    aP0(TPTP0)=aE(TE1TP1)aW(TP1TW1)+aN(UNUP)aS(UPUS)+SΔxΔy a_P^0(T_P-T_P^0)=a_E(T_E^1-T_P^1)-a_W(T_P^1-T_W^1)+a_N(U_N-U_P)-a_S(U_P-U_S)+S\Delta x\Delta y

    整理后得到:

    (aP0+aN+aS)TP1=aNTN1+aSUS1+(aP0aEaW)UP+aEUE+aEUW+SΔxΔy (a_P^0+a_N+a_S)T_P^1=a_NT_N^1+a_SU_S^1+(a_P^0-a_E-a_W)U_P+a_EU_E+a_EU_W+S\Delta x\Delta y

    在上式中,除了Tp1T_p^1TN1T_N^1TS1T_S^1,其他量都是已知的。令aP=aP0+aN+aSa_P=a_P^0+a_N+a_SD=(aP0aEaW)UP+aEUE+aEUW+SΔxΔyD=(a_P^0-a_E-a_W)U_P+a_EU_E+a_EU_W+S\Delta x\Delta y,则上式可以简化为:

    aPTP1=aNTN1+aNTN1+D a_PT_P^1=a_NT_N^1+a_NT_N^1+D

    再使用沿yy方向的TDMA求得下一时刻的温度场T1T^1

    边界节点的处理方式与上述过程相似,主要区别在于该节点的离散方程中未知项少了一项,例如对于顶部对流边界,aNa_N的取值得改为αΔx\alpha\Delta xTNT_N由未知的待求解量改为已知量TfT_f

    求解程序

    下面设计一个算例,并通过Fortran程序给出其求解过程:

    有一个宽400mm,高200mm的矩形计算域。计算域内材料的密度为2500kg/m3,热容为1372J/(kg·K)。含内热源,热生成率恒定为1050J/(m3·s)。导热系数为温度的多项式函数:

    λ=1.13588×1015T5+3.25358×1012T43.25305×109T3+1.32926×106T29.27637×105T+1.04478 \lambda=-1.13588\times 10^{-15}T^5+3.25358\times 10^{-12}T^4-3.25305\times 10^{-9}T^3+1.32926\times 10^{-6}T^2-9.27637\times 10^{-5}T+1.04478

    计算域的初始温度为1500℃,四周与环境进行对流换热,其中顶部对流给热系数α1\alpha_1为200W/(m2·K),底部和左右侧的对流给热系数α2\alpha_2为80W/(m2·K),环境温度TfT_f为20℃。模拟的物理时间长度为3h。

    采用均匀结构化网格,每个网格节点的尺寸为2mm×2mm。由于该算例的模型沿yy轴是对称的,实际计算域可以只取左半部分,即选取的计算域尺寸为200mm×200mm,网格节点数量为101×101。计算域右侧作为对称边界处理,该边界上TE=TWT_E=T_WaE=aWa_E=a_W

    由于Peaceman-Rachford ADI格式的时间步长选取不受稳定性限制,其时间步长的选取比较灵活。由于非稳态导热随着时间的进行,会趋向于稳态,刚开始时温度变化率会很大,随后则会大幅减小,本算例选取时间步长的准则是使选定的时间步长内最大温度变化幅度不超过10℃,这样刚开始时间步长会很小,随后则会变得很大。

    计算过程中选取了6个采样点,并设定每隔5分钟输出一次这6个采样点的温度值。

    求解程序如下,模拟用到的参数在网格材料边界条件初始条件求解控制模块中通过常量定义出来,以便于参数的调整:

    module mesh
        implicit none
        real, parameter :: width = 200e-3, height = 200e-3
        real, parameter :: delta_X = 2e-3, delta_Y = 2e-3
    end module mesh
    
    module material
        implicit none
        real, parameter :: rho = 2500, Cp = 1372, S = 1050
    end module material
    
    module boundary_condition
        implicit none
        real, parameter :: alpha1 = 300, alpha2 = 80
        real, parameter :: Tf = 20
    end module boundary_condition
    
    module initial_condition
        implicit none
        real, parameter :: T0 = 1500
    end module initial_condition
    
    module solution_control
        implicit none
        real, parameter :: time_length = 180 * 60
        real, parameter :: sampling_interval = 5 * 60
        real, parameter :: max_deltaT_per_step = 10
        real, dimension(6) :: sample_points_X = (/30e-3, 100e-3, 200e-3, 30e-3, 100e-3, 200e-3/)
        real, dimension(6) :: sample_points_Y = (/100e-3, 100e-3, 100e-3, 170e-3, 170e-3, 170e-3/)
        integer :: spIndex(6,2)
    end module solution_control
    
    
    program TDMA
        use mesh
        use material
        use solution_control
        implicit none
        integer :: Nx, Ny
        
        !确定X、Y方向的网格节点数量
        Nx = width / delta_X + 1
        Ny = height / delta_Y + 1
    
        !根据采样点位置获得采样点在网格节点中的索引
        spIndex(:,1) = floor(sample_points_X / delta_X + 1)
        spIndex(:,2) = floor(sample_points_Y / delta_Y + 1)
    
        !开始计算
        call calculate(Nx, Ny)
    
    end program TDMA
    
    
    subroutine calculate(Nx, Ny)
        use mesh
        use material
        use boundary_condition
        use initial_condition
        use solution_control
        implicit none
        integer, intent(in) :: Nx, Ny
        integer :: i
        real :: current_time, current_sampling_time, delta_Tau
        real, dimension(0:Nx+1,0:Ny+1) :: T, U, T1, delta_T
        real, dimension(Nx,Ny) :: lambda, aP0, aE, aW, aN, aS, aP, Sv, D, P, Q
        real :: kx(Ny), ky(Nx)
        real :: max_delta_T, max_rate
        real :: tStart, tEnd
    
        !初始化系数矩阵A
        aP0 = 0
        aE = 0
        aW = 0
        aN = 0
        aS = 0
    
        !初始化温度矩阵
        T = Tf
        T(1:Nx,1:Ny) = T0
        U = Tf
        T1 = Tf
    
        !初始化源项矩阵
        Sv = S * delta_X * delta_Y
        Sv(1,:) = Sv(1,:) / 2
        Sv(:,1) = Sv(:,1) / 2
        Sv(:,Ny) = Sv(:,Ny) / 2
    
        !初始化输出
        write(*,"(A4,$)") "Time"
        do i = 1, 5
            write(*,"(A9,I1,$)") "SP", i
        end do
        write(*,"(A10)") "SP6"
    
        !按时间步进行温度场计算
        current_time = 0
        current_sampling_time = sampling_interval
        delta_Tau = 0.0001  !考虑到降温刚开始时的温度变化率很大,初始delta_Tau取一个很小的值
        call CPU_TIME(tStart)
        do while (current_time < time_length)
            !更新当前时间
            current_time = current_time + delta_Tau
    
            !通过多项式函数计算材料导热系数
            lambda = -1.13588E-15 * T(1:Nx,1:Ny)**5 + 3.25358E-12 * T(1:Nx,1:Ny)**4 - 3.25305E-09 * T(1:Nx,1:Ny)**3 &
                   + 1.32926E-06 * T(1:Nx,1:Ny)**2 - 9.27637E-05 * T(1:Nx,1:Ny) + 1.04478
    
            !计算系数矩阵A
            !(1)内部节点的计算
            aP0 = rho * Cp * delta_X * delta_Y / (delta_Tau / 2)
            aE(1:Nx-1,:) = 2 / (1 / lambda(1:Nx-1,:) + 1 / lambda(2:Nx,:)) / delta_X * delta_Y
            aW(2:Nx,:) = 2 / (1 / lambda(2:Nx,:) + 1 / lambda(1:Nx-1,:)) / delta_X * delta_Y
            aN(:,1:Ny-1) = 2 / (1 / lambda(:,1:Ny-1) + 1 / lambda(:,2:Ny)) / delta_Y * delta_X
            aS(:,2:Ny) = 2 / (1 / lambda(:,2:Ny) + 1 / lambda(:,1:Ny-1)) / delta_Y * delta_X
            !(2)左侧、上下侧对流边界节点的处理
            aW(1,:) = alpha2 * delta_Y
            aN(:,Ny) = alpha1 * delta_X
            aS(:,1) = alpha2 * delta_X
            !(3)右侧对称边界的处理
            aW(Nx,:) = aW(Nx,:) * 2
            !(4)边界节点的减半处理
            !(4-1)上侧aP0、aE、aW减半
            aP0(:,Ny) = aP0(:,Ny) / 2
            aE(:,Ny) = aE(:,Ny) / 2
            aW(:,Ny) = aW(:,Ny) / 2
            !(4-2)下侧aP0、aE、aW减半
            aP0(:,1) = aP0(:,1) / 2
            aE(:,1) = aE(:,1) / 2
            aW(:,1) = aW(:,1) / 2
            !(4-3)左侧aP0、aN、aS减半
            aP0(1,:) = aP0(1,:) / 2
            aN(1,:) = aN(1,:) / 2
            aS(1,:) = aS(1,:) / 2
    
            !沿x方向使用TDMA隐式求解
            aP = aP0 + aE + aW
            !(1)定义TDMA中的常数矩阵
            D = (aP0 - aN - aS) * T(1:Nx,1:Ny) + aN * T(1:Nx,2:Ny+1) + aS * T(1:Nx,0:Ny-1) + Sv
            !(2)左侧边界的特殊处理
            D(1,:) = D(1,:) + aW(1,:) * Tf
            !!!!!!
            P(1,:) = aE(1,:) / aP(1,:)
            Q(1,:) = D(1,:) / aP(1,:)
            do i = 2, Nx
                kx = aP(i,:) - aW(i,:) * P(i-1,:)
                P(i,:) = aE(i,:) / kx
                Q(i,:) = (D(i,:) + aW(i,:) * Q(i-1,:)) / kx
            end do
            U(Nx,1:Ny) = Q(Nx,:)
            do i = Nx-1, 1, -1
                U(i,1:Ny) = P(i,:) * U(i+1,1:Ny) + Q(i,:)
            end do
    
            !沿y方向使用TDMA隐式求解
            aP = aP0 + aN + aS
            !(1)定义TDMA中的常数矩阵
            D = (aP0 - aE - aW) * U(1:Nx,1:Ny) + aE * U(2:Nx+1,1:Ny) + aW * U(0:Nx-1,1:Ny) + Sv
            !(2)上下侧边界的特殊处理
            D(:,1) = D(:,1) + aS(:,1) * Tf
            D(:,Ny) = D(:,Ny) + aN(:,Ny) * Tf
            !!!!!!
            P(:,1) = aN(:,1) / aP(:,1)
            Q(:,1) = D(:,1) / aP(:,1)
            do i = 2, Ny
                ky = aP(:,i) - aS(:,i) * P(:,i-1)
                P(:,i) = aN(:,i) / ky
                Q(:,i) = (D(:,i) + aS(:,i) * Q(:,i-1)) / ky
            end do
            T1(1:Nx,Ny) = Q(:,Ny)
            do i = Ny-1, 1, -1
                T1(1:Nx,i) = P(:,i) * T1(1:Nx,i+1) + Q(:,i)
            end do
    
            !当达到采样时间时,输出采样点温度
            if (current_time >= current_sampling_time) then
                !输出6个采样点的温度值
                write(*,"(I4,$)") floor(current_time / 60)
                do i = 1, 5
                    write(*,"(f10.2,$)") T1(spIndex(i,1), spIndex(i,2))
                end do
                write(*,"(f10.2)") T1(spIndex(6,1), spIndex(6,2))
                !
                !更新下一次需要采样的时间
                current_sampling_time = current_sampling_time + sampling_interval
            end if
    
            !计算最大冷却速率
            delta_T = abs(T1 - T)
            max_delta_T = maxval(delta_T)
            max_rate = max_delta_T / delta_Tau
    
            !根据对每个时间步内的最大温度变化率限制,以及下一次采样时间确定下一个时间步长
            delta_Tau = min(max_deltaT_per_step / max_rate, current_sampling_time - current_time)
    
            !更新温度矩阵T
            T = T1
    
        end do
        call CPU_TIME(tEnd)
        write(*,"('Calculate complete, time consumption: ', f8.4, ' sec')") (tEnd - tStart)
    
    end subroutine calculate
    

    参考文献

    [1] 陶文铨. 数值传热学[M]. 第二版. 西安: 西安交通大学出版社, 2001: 99-104.

    展开全文
  • ADI(交替方向隐格式)求解二维抛物方程(含matlab程序)其中包括了ADI算法的解析步骤和计算例子,最后附上MATLAB程序以供参考
  • 交替方向隐式方法是为了对二维区域内偏微分方程的差分求解而发展起来的一类方法。二维区域的偏微分方程,当使用隐式差分方法格式求解时,形成的五对角方程组如果用高斯赛德尔迭代法、线松弛迭代等方法求解,工作量...
  • 压力 Poisson方程用 Tschebyscheff SLOR方法交替方向迭代求解。N-S方程数值离散时对流项采用了 Chakravaythy-Osher TVD格式。用该方法计算后台阶流场的结果与经典的实验结果相当吻合,表明该方法是可靠的,在合适的...
  • 里面主要是三种二维ADI算法的matlab实现程序以及代码解释,没有具体算法步骤的分析。
  • 这组模拟使用交替方向隐式方法来求解 FD BPM 中出现的抛物线波动方程。 zip 文件包括以下程序: FDBPM3D_free_space.m - 3D 传播的动画和视频。 FDBPM3D_free_space_slice.m - 它显示了传播脉冲的切片。 FDBPM3D_...
  • 为解决大型稀疏矩阵的求解问题,采用一种改进的交替方向隐式方法和SSOR法提出两种不同的预处理器,并对Householder-GMRES(m)算法进行左端预处理,形成两种新算法,对算法的收敛性进行分析,给出数值算例验证新算法的可行...
  • 通过数值方法使用交替方向隐式方法(ADI)来处理所得的微分方程,这是有限差分方法之一,并且在两种情况下:非稳态和稳态。 并且我们完成了对以下各项的影响的研究:普朗特数,施密特数,格拉茨霍夫数和辐射参数。
  • 提出了一种求解二维波动方程的高精度紧致差分方法, 该方法首先利用紧交替方向隐式差分格式, 其截断误差为O(τ2 h4), 分别在粗网格和细网格上对原方程进行求解, 然后利用Richardson外推计算一次, 进一步提高精度, ...
  • 为了研究煤与瓦斯突出后瓦斯浓度的分布规律,基于流体力学理论,叙述了巷道内突出后瓦斯运移的数学模型,采用交替方向隐式格式方法(ADI方法)对瓦斯运移的对流扩散数学方程进行离散化处理,并在给定的物理模型以及...
  • 在这篇研究文章中,描述了两个有限差分隐式数值方案,以近似以耦合形式存在的二维修正反应扩散费舍尔系统的数值解。 有限差分隐式方案显示了计算算法的无条件稳定和二阶准确性质,并且通过具有已知解析解的示例完成...
  • 为解决微型无线传感器节点设计过程中无法分析射频微机电(RFMEMS)器件微小可动结构的电磁特性问题,提出一种电磁建模方法——变形时域交替方向隐式有限差分(DADI-FDTD)法。采用有限差分技术,引入时间变量实现一...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,996
精华内容 798
关键字:

交替方向隐式方法