在 Mathematica 中，您可以在多个级别上执行此操作。首先，让我们尝试解决它。我们需要几何、方程和边界条件：

reg = ImplicitRegion[0 <= x <= 2 && 0 <= y <= 0.5 && ! (x >= 1 && y <= 0.1) && ! (x >= 1 && y >= 0.4), {x,y}];
rp = RegionPlot[reg, AspectRatio -> Automatic]

接下来，我们需要方程和边界条件：

stokesFlowOperator = {
Div[{{-1, 0}, {0, -1}}.Grad[u[x, y], {x, y}], {x, y}] + Derivative[1,0][w][x, y], 
Div[{{-1, 0}, {0, -1}}.Grad[v[x, y], {x, y}], {x, y}] + Derivative[0,1][w][x, y], 
Derivative[0, 1][v][x, y] + Derivative[1, 0][u][x, y]};

bc = {DirichletCondition[u[x, y] == 4*0.3*y*(0.5 - y)/(0.41)^2,x == 0.], 
   DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 0 < x < 2], 
   DirichletCondition[w[x, y] == 0., x == 2]};

顶级解决方案：

好的，这就是现在解决我们使用的方程的设置NDSolve。

{xVel, yVel, pressure} = NDSolveValue[{stokesFlowOperator == {0, 0, 0}, bc}, {u, v, w}, {x, y} \[Element] reg, Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, w -> 1}}];

并可视化该区域上的速度矢量场：

Show[rp, StreamPlot[{xVel[x, y], yVel[x, y]}, {x, y} \[Element] reg]]

完毕。

更新

评论中提出的一个问题，这是可视化相同流图的另一种方法：

rmf = RegionMember[reg];
Show[rp, StreamPlot[{xVel[x, y], yVel[x, y]}, {x, 0, 2}, {y, 0, 0.5}, 
  AspectRatio -> Automatic, 
  RegionFunction -> Function[{x, y}, rmf[{x, y}]]]]

这给出了相同的结果。

低级访问：

现在，假设您想深入研究，您当然可以这样做，并且有很多方法可以做到这一点。我展示了一种方式。我们NDSolve用作方程预处理器并从中提取有限元数据。您可以做到这一点意味着您可以访问解决方案的每个阶段。

用作NDSolve方程预处理器：

{ndstate} = 
 NDSolve`ProcessEquations[{stokesFlowOperator == {0, 0, 0}, bc}, {u, 
   v, w}, {x, y} \[Element] reg, 
  Method -> {"FiniteElement", 
    "InterpolationOrder" -> {u -> 2, v -> 2, w -> 1}}]

这将返回一个NDSolve状态对象。我们加载有限元包：

Needs["NDSolve`FEM`"]

并提取和检查有限元数据：

femd = ndstate["FiniteElementData"];
sd = ndstate["SolutionData"][[1]];
femd["Properties"]
{"BoundaryConditionData", "FEMMethodData", "PDECoefficientData", \
"Properties", "Solution"}

所以方程预处理创建了几个数据结构。让我们提取并使用它们：

bcd = femd["BoundaryConditionData"];
methodData = femd["FEMMethodData"];
pded = femd["PDECoefficientData"];

这些是这种低级方法的关键步骤：

dpde = DiscretizePDE[pded, methodData, sd];
dbc = DiscretizeBoundaryConditions[bcd, methodData, sd];

获取系统矩阵：

{load, stiffness, damping, mass} = dpde["SystemMatrices"];

您可以查看它们，修改它们或其他什么。

MatrixPlot[stiffness]

要应用边界条件，请使用：

DeployBoundaryConditions[{load, stiffness, None, None}, dbc]

求解方程：

stationarySolution = LinearSolve[stiffness, load];

创建一个插值函数：

mesh = methodData["ElementMesh"];
offsets = methodData["IncidentOffsets"];
xVel = ElementMeshInterpolation[{mesh}, 
   stationarySolution[[offsets[[1]] + 1 ;; offsets[[2]]]]];

并可视化，例如，等值线图中的 x 速度：

ContourPlot[xVel[x, y], {x, y} \[Element] mesh, 
 AspectRatio -> Automatic, PlotRange -> All, 
 ColorFunction -> ColorData["TemperatureMap"], Contours -> 10]

有关有限元方法的概述，请NDSolve查看此处。

希望这可以帮助。

您当然可以在 Matlab 或类似程序中实现一个简单的算法。或者，您可以使用作为现有有限元软件包之一的现代数值方法的专业实现。例如，这里是deal.II的Stokes教程程序。（免责声明：我编写了这个程序。）我确信 libMesh、FEniCS 等也有类似的程序。

好吧，正如您所说，实现斯托克斯方程相当容易。我建议您自己做，只需阅读 Peyret（1987 年）“流体流动的计算方法”第 7 章的前 20 页，您就可以实现这样一个简单的模型（只需去掉对流项）。