How to implement point source or volume source in finite element implementations

Question

I'm trying to do a simple implementation to study the advection-diffusion-reaction dynamics in a straight pipe.

I have points positioned along the length of the pipe (blue dots in the image above).

I came across a blog (link) that explains how a point source is implemented. However, I am not really sure how to adapt it in this case.

The equation that captures the concentration dynamics is the following

The R in the last term of RHS corresponds to sink/source and the units is generally defined in mol/m^3 s. I want to define R = kc where the units of k is per sec. But I want to define R only at the points colored blue (i.e only as a point source/sink and not as a line source or sink).

I would like to ask for suggestions on how this implementation can be done.

EDIT: I would like to compare the FEM result obtained for 1D straight pipe simulation of convection-diffusion dynamics with point sinks located along the length of the pipe (problem posted above) against the analytical solution.

Suggestions on the analytical expression against which the numerical result can be compared will be really helpful.

Answer 1

Seen that you also have an account over at mathematica stackexchange I am going to show an implementation using Mathematica. This is certainly not the only way you can do it but hopefully gets you started.

We start by creating a 1D mesh

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

The mesh has points included at 1/3 and 2/3. This is where we will position the point sources later. I also refined the mesh, because I am going to add a convection term which will make this PDE convection dominant (large Peclet number).

We set up the dependent variable, the mass concentration $c$, time variable $t$ and the spatial independent variable $x$. Also some material properties are defined.

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

For the point sources the idea is to use a regularized delta function. The Dirac delta function poses a problem in numerical simulations with FEM implementation that is not prepared to take it as an input. This is because the Dirac delta function is singular at the source location $X_s$ and needs special treatment that differs from how FEM treats other terms. This can be avoided with an approximation to the Dirac delta function. The process of approximating the Dirac delta function is called regularization.

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}}]

More information about this can be found in the Heat Transfer tutorial. The process for mass transport is the same. So if you also want to volume sources you can proceed in the same way as mentioned in that tutorial.

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

The regularization function needs a parameter gamma, which is based on the minimal element diameter. Let's construct the point source at the first included point and visualize it.

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

You see that this is an idealized approximation to the delta distribution. Note that

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

The Comsol approach uses a weak formulation. That is a different approach. I am trying to show that (point) sources can be set up in different ways.

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 is the point source strength. The unit would be [$mol/(m^3 s)$], because the regularization function is dimension less the source has the same units as the overall PDE components, see here. The first point source is active all the time while the second only starts at a specific time but is stronger.

We set up the PDE and solve it.

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

Visualize / explore the result.

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

More information can be found in the Mass Transport tutorial and the MassTransportPDEComponent.

I hope this is useful to you.

Update for older version:

The code I posted above works for version 12.3. Since the OP has 12.0 the following is the same code, suitable for version 12.0. In essence the MassTransportPDEComponent is not available in 12.0. But that's not a problem as one can write the PDE manually. Here is the complete code:

This first part is the same as before:

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;

This part is where things change. I added a parameters section and added a manual form of the 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];

This should give you the same result in version 12.0.

Update: Cloud version

Another alternative is to use the wolfram cloud. Here is a cloud version of the notebook.

Answer 2

Contrary to what @user21's answer, I don't think that you need to do anything special for point loads. Let's see why.

A point load can be represented as a Dirac delta "function". So, in your case it would be something like

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

where $\rho$ is the intensity of the source and $x_i$ is the position.

If we use a weighted residuals for the finite element method you end up with the following integral for "body" sources

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

being $w(x)$ the weight function. If $R(x) = \rho \delta(x - x_i)$,

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

Thus, it suffices to evaluate your weight function in the point $x_i$. For "common" finite elements, we use interpolating functions. This implies the following

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

where $r_i$ represents the position of the point in the reference element.

Notice the following:

• For points that coincide with nodes (your case) the source contribute to exactly one entry in the $b$ vector. But for points inside elements it contributes to all the nodes related to that element, one for each function $\phi_j$.

• For linear elements $\phi_j(r_i) > 0$, but that's not the case for higher-order elements.

Regarding the question if one should or should not use point sources and its physical interpretation... It is a different matter, more related to the physics of the problem, I would argue. Respecting the units that you mention at the end of your question, I think that you don't provide enough information to assess it.