Mathematica中耦合PDE有限元分析的法向电流边界条件设置问题

核心思路：转化为标量电势PDE #

你的方程组可以通过代入转化为仅关于电势phi的二阶椭圆型PDE，大幅简化边界条件设置：

1. 电流密度表达式：$J_i = \sum_j c_{ij} \partial_j \phi$
2. 代入连续性方程$\sum_i \partial_i J_i = 0$，展开得到：
   $$\partial_x(c_{xx}\partial_x\phi + c_{xy}\partial_y\phi) + \partial_y(c_{yx}\partial_x\phi + c_{yy}\partial_y\phi) = 0$$

边界条件的具体实现 #

1. 绝缘边界（$\hat{n} \cdot J = 0$） #

该条件等价于法向电流通量为零，对应phi的边界通量条件：$\hat{n} \cdot (c \cdot \nabla\phi) = 0$。在Mathematica中直接用NeumannValue[0, 边界条件]实现，它专门用于指定通量型边界条件。

2. 源极/漏极边界 #

用DirichletCondition直接指定源极和漏极的电势值（比如源极设为$V_0$，漏极设为0）。

3. 复杂几何的边界区分 #

对于复杂形状，通过离散边界并识别特定区域区分源/漏和绝缘边界：

完整代码示例 #

(* 定义复杂几何区域 *)
reg = Polygon[{{0, 0}, {2, 0}, {2, 1}, {1, 2}, {0, 1}}];

(* 定义各向异性电导系数（可改为空间依赖函数） *)
cxx = 1; cxy = 0.1; cyx = 0.1; cyy = 1;
c = {{cxx, cxy}, {cyx, cyy}};

(* 推导电流密度和PDE *)
Jx = c[[1, 1]] D[phi[x, y], x] + c[[1, 2]] D[phi[x, y], y];
Jy = c[[2, 1]] D[phi[x, y], x] + c[[2, 2]] D[phi[x, y], y];
pde = D[Jx, x] + D[Jy, y] == 0;

(* 离散边界并提取源/漏段 *)
bnd = BoundaryDiscretizeRegion[reg];
edges = MeshPrimitives[bnd, 1];
sourceEdge = Select[edges, RegionMember[#, {0, 0.5}] &]; (* 左侧边为源 *)
drainEdge = Select[edges, RegionMember[#, {2, 0.5}] &]; (* 右侧边为漏 *)

(* 设置边界条件 *)
V0 = 1; (* 源极电压 *)
bc = {
  DirichletCondition[phi[x, y] == V0, RegionMember[sourceEdge, {x, y}]],
  DirichletCondition[phi[x, y] == 0, RegionMember[drainEdge, {x, y}]],
  NeumannValue[0, !RegionMember[sourceEdge, {x, y}] && !RegionMember[drainEdge, {x, y}]]
};

(* 求解PDE *)
sol = NDSolve[{pde, bc}, phi, {x, y} ∈ reg];

(* 可视化结果 *)
ContourPlot[phi[x, y] /. sol[[1]], {x, y} ∈ reg, ColorFunction -> "TemperatureMap", PlotLegends -> Automatic]
StreamPlot[{Jx, Jy} /. sol[[1]], {x, y} ∈ reg, StreamStyle -> Blue, PlotRange -> All]

关键注意事项 #

• 如果电导系数$c_{ij}$是空间变化的（比如随x,y变化），只需将其定义为函数（如cxx[x,y] = Exp[-x^2]），代码逻辑无需修改。
• 对于导入的外部几何（如STL文件），可以用Import["model.stl"]加载区域，再用BoundaryDiscretizeRegion处理边界，同样适用上述方法。
• 若边界段难以通过坐标识别，可结合MeshCoordinates和MeshCells手动指定边界单元的索引。

内容的提问来源于stack exchange，提问作者Rain