Page 103 - 2021年第52卷第8期
P. 103

将式(11)代入式(10)左端,将式(8)代入式(10)右端,可得:
                     γ             γ
                                                                         )
                                                                                               )
                                                    )
                                     é
                     1              1
                                                                                                    ù
                        V Ψ dxdy =  å  Ψ Ψ V ( x ,y dxdy +   Ψ Ψ V ( x ,y dxdy +   Ψ Ψ V ( x ,y dxdy (12)
                    å x    I        ë   I  i  x  i  i        I  j  x  j  j         I  k  x  k  k   û
                     1             1
                       γ  1                γ  1  é  ∂Ψ                 ∂Ψ                ∂Ψ        ù
                                                             i 
                                                                                j 
                     - å  K x ∂H  Ψ dxdy = - å ê  x  I  i dxdyH +  K Ψ I  j  dxdyH +  K Ψ I  k dxdyH k ú ú  (13)
                                             ê K Ψ
                                  I
                                                                   x
                                                                                      x
                       1      ∂x           1  ë      ∂x                ∂x                 ∂x       û
                   将式(12)、式(13)结合式(5),可得 Zone 1 上关于达西流速 V 的方程组的具体形式,求解可以获
                                                                          x
               得 Zone 1 区域上的 V 值,其中包括了在 Zone 1 的边界上,即在 Zone 1 一侧的介质交界面 AB 上的 V                            x
                                 x
               值。结合 Jump 函数式(3),即可获 Zone 2 一侧的介质交界面 AB 上的 V 值,这将作为 Zone 2 上的子问
                                                                              x
               题的第一类边界条件。在 Zone 2 上,运用类似上述达西流速的模拟过程,即可获得 Zone 2 区域的粗
               尺度 V 值。在所有子问题求解完毕后,可以得到 Ω 上的 V 值。在Ω进行类似上述 V 的达西流速的模
                    x
                                                                                           x
                                                                    x
               拟过程,即可以得到Ω上的粗尺度 V 值。
                                               y
                   上述过程中,MSFEM-D 运用了多尺度基函数模拟地下水流和达西流速,得到了水头和达西流速
               的粗尺度网格上的解。因而,MSFEM-D 的水头和达西流速的总方程组的阶数均远低于精细剖分的传
               统方法,具其有极高的计算效率。
               3  数值模拟实验
                   采用如下缩写方式:V 表示 x 方向上的达西流速;V 表示 y 方向上的达西流速;FEM-F 表示精细
                                                                  y
                                       x
               剖分的有限元法;Model-Yeh-F 表示精细剖分的 Yeh 的伽辽金有限元模型;Model -Zhou-F 表示精细
               剖分的 Zhou 等的 S 有限元模型;MSFEM-D 表示多尺度有限元与区域分解法的组合方法。
                               1
                   本节模拟了含水平介质交界面的渗透系数非均质各向异性的地下水二维稳定流问题以及含水平
               与倾斜两个介质交界面的地下水稳定流问题的水头和达西流速,并将 MSFEM-D 所模拟的水头与
                                                                                       [3]
                                                                    [8]
               FEM-F 比较,MSFEM-D 所模拟的达西流速与 Model-Yeh-F 以及 Model-Zhou-F 两种现有达西流速
               算法进行比较。为了比较计算消耗与精度,FEM-F、Model-Yeh-F、Model-Zhou-F 的网格数目均等
               于 MSFEM-D 的细网格数目,且所有数值方法均采用了线性基函数的边界条件。同时,所有数值方法
               均采用 C++编写,均采用 Cholesky 分解法求解总方程的矩阵,没有应用并行计算技术。在测量 CPU
               时间时,各方法均在同一计算机上独立运行。
               3.1  含水平交界面的渗透系数非均质各向异性的二维地下水稳定流问题                                 本例基于 Cainelli 等在 2012
                     [4]
               年工作 中的数值算例所设置,具有解析解。本例的地下水水流和达西流速的控制方程分别为,二维
                                æ  ∂ æ   ∂H  ö  ∂ æ  ∂H  ö    ö
               地下水稳定流方程 ç ç      -   ç K x  ÷ -   ç K y  ÷ = W  ÷ ÷  和达西方程(式(9))。研究区域Ω为[0,1]×[0,
                                è  ∂x è  ∂x  ø  ∂y  è  ∂y  ø  ø
               1]。研究区内含有两种不同的各向同性含水介质(图 2),分别位于 Zone 1 和 Zone 2 两个区域,它们的
               交界面为 AB(y=0.5)。渗透系数呈非均质各向异性变化,Zone 1、Zone 2 的渗透系数分别为:K =(1+
                                                                                                     1x
               x)(1+y),K =10(1+x)(1+y)和 K =10(1+x)(1+y),K =100(1+x)(1+y)。本例中,地下水水流方程中
                         1y
                                                               2y
                                            2x
               的源汇项 W 和研究区的边界的第一类边界条件由如下水头 H 的解析解给出:
                                                             )
                                                                 3
                                                                   4
                                          ìZone 1:H = sin(10xy - x y + 1
                                          ï
                                          ï
                                                              )
                                          í             é(2y + 9 x  ù  3æ 2y + 9 ö 4                  (14)
                                          ï Zone 2:H = sinê ê   ú - x ç    ÷ + 1
                                                                ú
                                          ï
                                          î             ë   2   û    è  20  ø
                   此外,x 方向上的达西流速 V 的解析解为:
                                             x
                                                               )
                                    ì Zone 1:V = -K [ 10ycos(10xy - 3x y  4 ]
                                                                    2
                                    ï         x    x
                                    ï
                                                                     )
                                    í               ì      )   é(2y + 9 x ù  2 æ 2y + 9 ö  4 ù        (15)
                                                    ï(2y + 9
                                    ïZone 2:V = -K í        cos ê ê    ú - 3x ç     ÷  ú ú
                                                                       ú
                                    ï         x    2x  ï  2    ë   2   û     è  20  ø
                                    î               î                                û
                                                                                               — 983  —
   98   99   100   101   102   103   104   105   106   107   108