From the θ-Scheme to the Stokes Saddle-Point System

Why a matrix formulation in the first place

We want to solve the instationary Stokes equations, but a PDE that carries a time derivative cannot be handed to a stationary linear solver directly. The standard trick is a $θ$ -scheme: discretize in time first, turning the single time-dependent problem into a sequence of stationary problems, one per time step $n = 1, \dots, N$ . Starting from the known initial data $(v_{0}, p_{0})$ — here $v_{0} = 0$ , or in general whatever our measurements give us — each step advances the solution by one $k$ .

It is tempting to hope for a clean propagation rule of the form

(\begin{matrix} v_{n} \\ p_{n} \end{matrix}) = A (\begin{matrix} v_{n - 1} \\ p_{n - 1} \end{matrix}), n = 1, \dots, N,

but this is not quite what comes out, and the reason is worth keeping in mind. The pressure has no time derivative — it is a Lagrange multiplier enforcing incompressibility — so $p_{n - 1}$ never enters the right-hand side; only the velocity is propagated. And rather than a single matrix–vector product, each step requires solving a linear system. What we actually obtain is

(\begin{matrix} V^{'} \\ P^{'} \end{matrix}) = A^{- 1} (\begin{matrix} F (V) \\ O \end{matrix}),

i.e. a fixed operator $A^{- 1}$ applied at every step — “something alike”, just with an inverse and with velocity as the only memory between steps. Deriving $A$ and $F (V)$ explicitly is the goal below.

Strong form

After the $θ$ -scheme has been applied and every $n$ -term collected on the left, every $(n - 1)$ -term on the right, the stationary problem solved at each step reads:

\begin{aligned} \frac{1}{k} v_{n} - ν θ Δ v_{n} + \nabla p_{n} & = f (t_{n}) + \frac{1}{k} v_{n - 1} + ν (1 - θ) Δ v_{n - 1}, \\ \nabla \cdot v_{n} & = 0. \end{aligned}

The right-hand side is entirely known data: the load $f (t_{n})$ plus contributions from the already-computed previous velocity $v_{n - 1}$ .

Weak form

Testing the momentum equation against $u \in V_{0}$ (vanishing on $\partial Ω$ ) and the continuity equation against $q \in Q$ , then integrating the Laplacian and pressure-gradient terms by parts — the boundary terms drop out since $u |_{\partial Ω} = 0$ — we obtain:

\begin{aligned} \frac{1}{k} \int_{Ω} v_{n} \cdot u d x + ν θ \int_{Ω} \nabla v_{n} : \nabla u d x - \int_{Ω} p_{n} (\nabla \cdot u) d x & = \underset{known right-hand side}{\underset{⏟}{\int_{Ω} f (t_{n}) \cdot u d x + \frac{1}{k} \int_{Ω} v_{n - 1} \cdot u d x - ν (1 - θ) \int_{Ω} \nabla v_{n - 1} : \nabla u d x}}, \\ \int_{Ω} q (\nabla \cdot v_{n}) d x & = 0. \end{aligned}

Discretization with shape functions

As before, let $J_{v}$ denote the index set of the local degrees of freedom of the velocity space and $J_{p}$ that of the pressure space. Choose shape functions ${φ_{j}}_{j \in J_{v} \cup J_{p}}$ , where ${φ_{j}}_{j \in J_{v}} \subset V_{h}$ forms a basis of $V_{h}$ and ${φ_{j}}_{j \in J_{p}} \subset Q_{h}$ forms a basis of $Q_{h}$ . Let the discrete solution take the form

v_{n} = \sum_{j \in J_{v}} V_{j}^{'} φ_{j}, p_{n} = - \sum_{j \in J_{p}} P_{j}^{'} φ_{j} .

The minus sign in the pressure ansatz is a deliberate convention: it will flip the $- \int p_{n} (\nabla \cdot u)$ term to a plus, so that the momentum row couples to $+ B^{T}$ and matches the $+ B$ of the continuity row — producing a symmetric block structure instead of an antisymmetric one.

Inserting the ansatz and testing against $u = φ_{i}$ (for $i \in J_{v}$ ) and $q = φ_{i}$ (for $i \in J_{p}$ ) gives:

\begin{aligned} \sum_{j \in J_{v}} V_{j}^{'} \underset{A_{i j}}{\underset{⏟}{(\frac{1}{k} \int_{Ω} φ_{j} \cdot φ_{i} d x + ν θ \int_{Ω} \nabla φ_{j} : \nabla φ_{i} d x)}} + \sum_{j \in J_{p}} P_{j}^{'} \underset{B_{j i} = (B^{T})_{i j}}{\underset{⏟}{\int_{Ω} φ_{j} (\nabla \cdot φ_{i}) d x}} & = \underset{F (V)_{i}}{\underset{⏟}{\int_{Ω} f (t_{n}) \cdot φ_{i} d x + \frac{1}{k} \int_{Ω} v_{n - 1} \cdot φ_{i} d x - ν (1 - θ) \int_{Ω} \nabla v_{n - 1} : \nabla φ_{i} d x}}, \\ \sum_{j \in J_{v}} V_{j}^{'} \underset{B_{i j}}{\underset{⏟}{\int_{Ω} φ_{i} (\nabla \cdot φ_{j}) d x}} & = 0. \end{aligned}

Here $F (V)_{i}$ is the load vector depending on the previously computed velocity $V = V^{n - 1}$ . Note that the same mass and stiffness integrals that build $A$ also build $F (V)$ , so in practice the right-hand side is just $(\frac{1}{k} M - ν (1 - θ) A_{stiff}) V^{n - 1} + F_{ext}$ applied to the previous solution vector.

Matrix formulation

Collecting the blocks, we arrive at the Stokes saddle-point system solved at each time step:

(\begin{matrix} A & B^{T} \\ B & O \end{matrix}) (\begin{matrix} V^{'} \\ P^{'} \end{matrix}) = (\begin{matrix} F (V) \\ O \end{matrix}) .

A couple of remarks to close the loop:

The block $A$ here is the full velocity–velocity block $A = \frac{1}{k} M + ν θ A_{stiff}$ , not the bare stiffness matrix. In the stationary case the mass term $\frac{1}{k} M$ is absent and $A$ reduces to pure stiffness; in the instationary case it “fattens up”, but the shape, symmetry, and zero block stay exactly the same.
The zero block in the bottom-right reflects that pressure has no self-coupling — it only acts as the constraint multiplier — which is precisely why $p_{n - 1}$ never propagates and why each step is a solve of this indefinite system rather than a simple matrix product.

本文采用署名-非商业性使用-相同方式共享 4.0 国际许可协议，转载请注明出处。