与@user21 的回答相反，我认为您不需要为点负载做任何特别的事情。让我们看看为什么。

点载荷可以表示为狄拉克增量“函数”。所以，在你的情况下，它会像

$$R = \rho \delta(x - x_i)\, ,$$

其中是源的强度，是位置。$\rho$$x_i$

如果我们对有限元方法使用加权残差，您最终会得到以下“体”源的积分

$$b = \int_\Omega R(x) w(x)\, ,$$

为权重函数。如果，$w(x)$$R(x) = \rho \delta(x - x_i)$

$$b = \int_\Omega \rho \delta(x - x_i) w(x) = \rho w(x_i)\, .$$

因此，在点评估您的权重函数就足够了。对于“普通”有限元，我们使用插值函数。这意味着以下$x_i$

$$b = \begin{cases} \rho\quad \text{for nodes}\, ,\\ \rho \phi_j(r_i)\quad \text{otherwise}\, , \end{cases}$$

其中表示参考元素中点的位置。$r_i$

请注意以下事项：

• 对于与节点重合的点（您的情况），源仅对向量中的一个条目做出贡献。但是对于元素内的点，它对与该元素相关的所有节点都有贡献，每个函数一个。$b$$\phi_j$

• 对于线性元素，但对于高阶元素并非如此。$\phi_j(r_i) > 0$

关于一个人是否应该使用点源及其物理解释的问题......我认为这是另一回事，与问题的物理学更相关。尊重您在问题末尾提到的单位，我认为您没有提供足够的信息来评估它。

看到您在mathematica stackexchange上也有一个帐户，我将展示一个使用Mathematica 的实现。这当然不是你能做到的唯一方法，但希望能让你开始。

我们首先创建一个 1D 网格

Needs["NDSolve`FEM`"]
region = Line[{{0}, {1}}];
includePoints = {{1/3}, {2/3}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints, 
  "MaxCellMeasure" -> 0.0008]

网格在 1/3 和 2/3 处包含点。这是我们稍后将放置点源的地方。我还改进了网格，因为我要添加一个对流项，这将使这个 PDE 对流占主导地位（大 Peclet 数）。

我们设置了因变量、质量浓度、时间变量和空间自变量。还定义了一些材料属性。$c$$t$$x$

vars = {c[t, x], t, {x}};
pars = <|"DiffusionCoefficient" -> 1, 
   "MassConvectionVelocity" -> {100}|>;

对于点源，想法是使用正则化 delta 函数。Dirac delta 函数在使用 FEM 实现的数值模拟中提出了一个问题，该函数不准备将其作为输入。这是因为 Dirac delta 函数在源位置处是奇异的，并且需要与 FEM 处理其他项的方式不同的特殊处理。这可以通过对 Dirac delta 函数的近似来避免。逼近狄拉克δ函数的过程称为正则化。$X_s$

RegularizedDeltaPoint[g_, X_List, Xs_List] := 
 Piecewise[{{Times @@ Thread[1/(4 g) (1 + Cos[\[Pi]/(2 g) (X - Xs)])],
     And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}]

有关这方面的更多信息，请参阅传热教程。大众运输的过程是相同的。因此，如果您还想对源进行批量处理，您可以按照该教程中提到的相同方式进行操作。

Subscript[h, mesh] = Sqrt[Min[mesh["MeshElementMeasure"]]];
Subscript[gamma, reg] = Subscript[h, mesh]/2;

正则化函数需要一个基于最小元素直径的参数 gamma。让我们在第一个包含点处构建点源并将其可视化。

temp = RegularizedDeltaPoint[Subscript[gamma, reg], {x}, 
  includePoints[[1]]];
Plot[temp, {x, 0, 1}]

您会看到这是对 delta 分布的理想化近似。注意

Integrate[temp, {x, 0, 1}]
1.

Comsol 方法使用了一个弱公式。那是一种不同的方法。我试图表明（点）源可以以不同的方式设置。

Qp = 10;
pars["MassSource"] = 
  RegularizedDeltaPoint[Subscript[gamma, reg], {x}, 
     includePoints[[1]]]*Qp +
   RegularizedDeltaPoint[Subscript[gamma, reg], {x}, 
     includePoints[[2]]]*If[t > 1/2*10^-3, 5*Qp, 0];

Qp是点源强度。单位将是 [ ]，因为正则化函数的维度少，源具有与整体 PDE 组件相同的单位，请参见此处。第一个点源始终处于活动状态，而第二个点源仅在特定时间开始但更强。$mol/(m^3 s)$

我们设置 PDE 并解决它。

pde = {MassTransportPDEComponent[vars, pars] == 0, c[0, x] == 0};
tEnd = 10^-3;
cfun = NDSolveValue[pde, c, {t, 0, tEnd}, {x} \[Element] mesh];

可视化/探索结果。

Manipulate[
 Plot[cfun[t, x], {x} \[Element] region, 
  PlotRange -> {{0, 1}, {0, 1/2}}], {t, 0, tEnd}]

更多信息可以在MassTransport教程和MassTransportPDEComponent中找到。

我希望这对你有用。

旧版本更新：

我在上面发布的代码适用于 12.3 版。由于 OP 有 12.0 以下是相同的代码，适用于 12.0 版本。本质上MassTransportPDEComponent12.0 中不可用。但这不是问题，因为可以手动编写 PDE。这是完整的代码：

第一部分与以前相同：

RegularizedDeltaPoint[g_, X_List, Xs_List] := 
  Piecewise[{{Times @@ 
      Thread[1/(4 g) (1 + Cos[\[Pi]/(2 g) (X - Xs)])], 
     And @@ Thread[RealAbs[X - Xs] <= 2 g]}, {0, True}}];

Needs["NDSolve`FEM`"]
region = Line[{{0}, {1}}];
includePoints = {{1/3}, {2/3}};
mesh = ToElementMesh[region, "IncludePoints" -> includePoints, 
   "MaxCellMeasure" -> 0.0008];

Subscript[h, mesh] = Sqrt[Min[mesh["MeshElementMeasure"]]];
Subscript[gamma, reg] = Subscript[h, mesh]/2;

这部分是事情发生变化的地方。我添加了一个参数部分并添加了 PDE 的手动形式。

parameters = {kappa -> {{1}}, v1 -> 100, 
   gamma -> Subscript[gamma, reg], Qp -> 10};
pde = {Derivative[1, 0][c][t, x] + 
      Inactive[Div][(-kappa) . 
        Inactive[Grad][c[t, x], {x}], {x}] + {v1} . 
       Inactive[Grad][c[t, x], {x}] - 
      Qp*RegularizedDeltaPoint[gamma, {x}, {1/3}] - 
      If[t > 1/2000, 5*Qp, 0]*
       RegularizedDeltaPoint[gamma, {x}, {2/3}] == 0, 
    c[0, x] == 0} /. parameters;

tEnd = 10^-3;
cfun = NDSolveValue[pde, c, {t, 0, tEnd}, {x} \[Element] mesh];

这应该在版本 12.0 中为您提供相同的结果。

更新：云版

另一种选择是使用 wolfram 云。这是笔记本的云版本。