令

则方程(1)的变分形式为：求(u，p)∈(V，M)，满足

其中$a\left( u, v \right)=\int_{\mathit{\Omega }}{\nabla u\nabla v\text{d}x\text{d}y, b\left( v, p \right)}=\int_{\mathit{\Omega }}{p\nabla v\text{d}x\text{d}y}, \left( f, v \right)=\int_{\mathit{\Omega }}{fv\text{d}x\text{d}y}$.由文献[1]知空间(V，M)满足LBB条件，方程(2)的解存在且唯一.

为简单起见，设J_h是区域Ω中单元最大直径为h的三角形(或四边形)正则剖分簇，即$\overline{\mathit{\Omega }}=\bigcup\limits_{K\in {{J}_{h}}}{K}$.定义有限元空间

其中R₁(K)表示定义在单元K上的P₁元或Q₁元.由于空间对(V_h，M_h)不满足离散的LBB条件，所以我们利用文献[10-11]中的方法，引入相应的稳定项来构造如下稳定化求解格式：

求(u_h，p_h)∈(V_h，M_h)，满足

其中R_hu₀是初值u₀(x)的投影，将在下文(6)式中定义.

其中Π：L²(Ω)→R₀是标准的L²投影算子，满足性质

其中R₀表示单元K上的分片常数.这里及下文中的C均为与h无关的正常数，且取值可以不同.

为了得到有限元解(u_h，p_h)的最优误差估计结果，引入如下投影算子：

其中

投影算子(R_h，Q_h)是适定的.

下面我们将对稳定化格式(3)进行误差分析，并得到如下结论：

定理1  设(u，p)和(u_h，p_h)分别为方程(2)和(3)的解，(u，p)∈((H²(Ω))²∩V，H¹(Ω)∩M)，则有误差估计

其中σ(t)=max{1，${{t}^{-\frac{1}{2}}}$}.

证  令

在方程(2)中取(v，q)=(v_h，q_h)，并与(3)式相减，得

利用投影定义(6)式可得

在(10)式中取(v_h，q_h)=(θ，φ)，可得

进而有

对(11)式两端从0到t进行积分，注意到θ(0)=0，可得

再利用Gronwall不等式得

最后结合三角不等式以及文献[11]中引理4.1的结论可得(7)式.

(10) 式两边同时对t求导，得

在(13)式中取(v_h，q_h)=(θ_t，φ_t)，并利用Young不等式，得

将(14)式的两端乘以t，再同时从0到t进行积分，得

进一步利用Gronwall不等式，可得

由三角不等式和文献[11]中的引理4.1可以推出(8)式.

下面对压力p的误差进行分析，得出如下结论：

定理2  在定理1的条件下，有

证  由文献[12]中的定理3.1，可得

再利用三角不等式、定理1和文献[11]中引理4.1的结论可以推出(17)式.

为了验证理论分析的正确性，我们给出下列数值实验.数值实验1对本文使用的局部压力投影稳定化方法进行了数值验证，测试收敛阶是否符合理论分析结果；数值实验2给出了几种满足离散的LBB条件的常用低阶混合有限元方法的计算结果，以便于进行对比分析.

在方程(1)中取区域Ω=[0, 1]×[0, 1]，其精确解为

右端项可以通过计算得到.为比较方便，我们统一使用向后欧拉格式进行计算，计算中采用直角三角形或正方形网格剖分.

数值实验1  基于本文稳定化方法，分别利用同阶P₁²-P₁元和Q₁²-Q₁元进行数值求解.在t=1时的计算结果见表 1、表 2.

由表 1、表 2可以看到速度的H¹模和压力的L²模的误差收敛阶均超过1阶，这也进一步验证了理论分析的正确性.

数值实验2  为了便于和本文方法对比，我们同时也给出了几个常用的满足离散的LBB条件的低阶混合元方法求解方程(1)的结果.其中，表 3给出的是协调双线性Q₁元配常数Q₀元对Q₁²-Q₀的计算结果，表 4是非协调旋转Q₁元^[18]配常数Q₀元对Q₁^rot-Q₁^rot-Q₀的计算结果，表 5是文献[5-6]中使用的CR型非协调三角形元配常数P₀元对P₁^nc-P₁^nc-P₀的计算结果，表 6是文献[7]中使用的CR型非协调矩形元配常数Q₀元对MQ₁^rot-Q₀的计算结果.

通过对结果进行分析，我们发现上述方法关于速度u的能量模的收敛阶都能达到最优.相对而言，利用同阶Q₁²-Q₁元对计算压力p时，会有高半阶的超收敛结果.关于数值实验2中所列的单元在求解流体问题时各自所具有的优势可查阅相关文献，本文不再详述.