有时间项和灌注源项的 Stokes 方程，用 Talor-Hood 方法实现

stokes2D_2.py

# import matplotlib.pyplot as plt
from dolfin import *

mesh = UnitSquareMesh(100, 100)
T = 1.0
Nt = 100
dt = T/Nt
nu = 0.01

# Define function spaces
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)

# No-slip boundary condition for velocity
# Boundaries
def right(x, on_boundary): return x[0] > (1.0 - DOLFIN_EPS)
def left(x, on_boundary): return x[0] < DOLFIN_EPS
def bottom(x, on_boundary): return x[1] < DOLFIN_EPS
def top(x, on_boundary): return x[1] > (1.0 - DOLFIN_EPS)

# No-slip boundary condition for velocity
noslip = Constant((0.0, 0.0))
bc3 = DirichletBC(W.sub(0), noslip, right)
bc1 = DirichletBC(W.sub(0), noslip, left)
bc2 = DirichletBC(W.sub(0), noslip, bottom)

# Inflow boundary condition for velocity
topflow = Expression("1", degree=1)
bc0 = DirichletBC(W.sub(0).sub(0), topflow, top)

no_pressure = Constant(0.0)
bc4 = DirichletBC(W.sub(1), no_pressure, top)

# Collect boundary conditions
bcs = [bc0, bc1, bc2, bc3, bc4]

# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

f = Constant((0, 0))
s = Expression('-100*x[1]*(x[1]-1)*x[0]*(x[0]-1)', degree=5, t=0)
k = Constant(1/dt)
w_ = Function(W)
(un, _) = Function(W).split(True)


F = k*inner((u - un), v) * dx + (nu*inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx - inner(f, v)*dx - s*q*dx
a = lhs(F)
L = rhs(F)

t = 0
ufile = File("result/velocity.pvd")
pfile = File("result/pressure.pvd")
for it in range(Nt):
    t += dt
    s.t = t
    print("t = ", t)
    # Assemble matrix and vector
    A = assemble(a)
    b = assemble(L)
    # Compute solution
    solve(a == L, w_, bcs)
    (u_, p_) = w_.split(True)
    ufile << u_
    pfile << p_
    un.assign(u_)

Nx = 100
Ny = 100
i = 5

# for i in range(Nx):
for j in range(Ny):
    # print(i/Nx+0.5/Nx, j/Ny)
    result = "{:.16e}".format(u_(i/Nx+0.5/Nx, j/Ny)[1])
    print(result)