Say you want to solve the movement of a particle in a viscous fluid such that the force on the particle is zero. In general, the solution is trivial since one has no-slip boundary conditions, ie the fluid moves at the same velocity as the particle.

But supposed you do have a slip velocity at the surface of the particles due to some microscopic phenomena. (See here for the paper) In particular, if the particles are in a varying field (electric, temperature, concentration of any solute) the velocity of the fluid at the surface of the particle is proportional to the (tangent) gradient of the field in that point.

The question is then, how to solve Stokes equation for those particles with the boundary conditions on the particles such that the force on them is zero and they move with a certain (unknown) velocity.

That was the physical introduction, now the code. First I review the simple boundary conditions for PDE, then present examples of integral boundary conditions, and finally the differential-integral boundary condition required by force=0.

<< NDSolve`FEM`

L = 5;
r = RegionDifference[Rectangle[{-L, -L}, {L, L}], Disk[{0, 0}, 1]];

mesh = ToElementMesh[ r ];

eqn = {
   Laplacian[u[x, y], {x, y}] - D[p[x, y], x] == 0,
   Laplacian[v[x, y], {x, y}] - D[p[x, y], y] == 0,
   D[u[x, y], x] + D[v[x, y], y] == 0};

(*with null boundary conditions in "infinite"*)

bcs = {DirichletCondition[u[x, y] == 0, 
    x == L || x == -L || y == L || y == -L],
   DirichletCondition[v[x, y] == 0, 
    x == L || x == -L || y == L || y == -L],
   DirichletCondition[p[x, y] == 0., x == L || x == -L]

One can solve for Dirichlet condition on the particle like this

pde = {eqn, bcs, 
   DirichletCondition[u[x, y] == 1, Norm[{x, y}]^2 <= 1.], 
   DirichletCondition[v[x, y] == 1, Norm[{x, y}]^2 <= 1.]};

{xVel, yVel, pressure} = 
  NDSolveValue[pde, {u, v, p}, {x, y} ∈ mesh2, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

And it works.

Now, say you want to impose that the integral of the function over the particle is constant, something like this:

   pde = {eqn, bcs,NIntegrate[u[x, y] + v[x, y], {x, y} ∈ Circle[]]==1};

{xVel, yVel, pressure} = 
  NDSolveValue[pde, {u, v, p}, {x, y} ∈ mesh2, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

This is what I understand for an integral boundary condition and this won't work. For what I could find online seems that mathematica cannot deal with this atm.

Now, the condition of force=0 on the particle is slightly more complicated. First we need to remember that the force acting on a surface due to the fluid is $F_i = \sigma_{ik}dA_k$. For the sake of simplicity let's say that the force is only in x direction, so for the case of our particle we have that

force = NIntegrate[(-pressure[Cos[ϕ], Sin[ϕ]] + 
     Subscript[σ, xx][Cos[ϕ], Sin[ϕ]]) Cos[ϕ] + 
  Subscript[σ, xy][Cos[ϕ], 
    Sin[ϕ]] Sin[ϕ], {ϕ, 0, 2 π}]

is the force acting on the particle, with the stress tensor components given by

Subscript[σ, yy][x_, y_] = Evaluate[2 D[yVel[x, y], {y}]];
Subscript[σ, xx][x_, y_] = Evaluate[2 D[xVel[x, y], {x}]];
Subscript[σ, xy][x_, y_] = 
 Evaluate[D[xVel[x, y], {y}] + D[yVel[x, y], {x}]]

That's an integro-differential boundary condition and it's what I need to solve. I have a way to to do this which consists of iteratively, given the slip velocity on the particle, I start with a velocity, solve Stokes equation for that boundary condition, and compute the force. If the force is negative I need a larger velocity and if it's positive a smaller velocity. Eventually I reach the velocity where the force on the particle is zero.

However, what I would like to do is to pass the boundary condition to NDSolve just like this:

 pde = {eqn, bcs, force==0};
{xVel, yVel, pressure} = 
  NDSolveValue[pde, {u, v, p}, {x, y} ∈ mesh2, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]

too much to ask?

  • $\begingroup$ Could you clarify this: "...impose that the integral of the function over the particle is constant..." Do you mean the integral over the sum of the velocity components (u + v) over the boundary of the Disk (i.e. a Circle? The Disk is what you have subtracted from the region why do you want to integrate over it now? $\endgroup$
    – user21
    Commented Apr 4, 2016 at 10:02
  • $\begingroup$ Do I understand correctly that you have two questions, one is about how to implement NIntegrate[u[x,y]+v[x,y], Element[{x,y}, Circle[]]]==1 as a boundary condition. And the second is to compute some force from the resulting interpolation functions? Or is this force feed back (i.e. non-linear) into the boundary integral constraint? $\endgroup$
    – user21
    Commented Apr 4, 2016 at 10:06
  • $\begingroup$ The first integral is just an example of an integral boundary condition. What I do want to implement is the integral of the force, which is a function of the spatial derivatives of the fluid velocity and the pressure. I thought that dividing the problem in smaller steps would be easier to understand. Sorry for the typos. $\endgroup$ Commented Apr 4, 2016 at 13:21

1 Answer 1


I'm answering my own question since it may eventually help somebody else who is facing a similar problem.

In this case since the equations are linear on the particles' velocities one can find "easily" the one that makes the force zero.

This is how I'm doing it.

First, let's put the problem in abstract terms and for only one dimension for the sake of simplicity.

We have that the force on the particle $i$ is linear in each of the particles' velocities.

$$ F_i = \sum_j m_{ij} U_j + b_i$$

with $j$ running over all the particles. Evaluating the force for the "homogeneous" solution, ie. $\vec{U} = 0$ on has that $b_i = F_i(\vec{0})$. For the two particles case in turn

\begin{eqnarray} F_1 = m_{11} U_1 + m_{12} U_2 + b_1~,\\ F_2 = m_{21} U_1 + m_{22} U_2 + b_2~. \end{eqnarray} To find the constants $m_{ij}$ we evaluate the force for different $U_i$ linearly independent, even better by an orthonormal base in the velocity space. $m_{ij}$ are obtained then by solving this system of equations. Since there are four unknowns and four equations the system is solvable.

Having found the dependence of the force on the particles' velocities one can impose now that the force on both particles is zero: $$ M.U = -b$$ whose solution for the velocity is given in terms of the force evaluated at the base vectors. In Mathematica:

Clear[m, f]
n = 3
e0 = 0 Range@n;
eaux = e0;
eaux[[1]] = 1;

e = Permutations[eaux]
sol = Solve[ 
    Table[f[i][e[[j]]] == m[i, j] + f[i][e0], {j, n}, {i, n}], 
  Flatten@Table[m[i, j], {i, n}, {j, n}]]

LinearSolve[Table[First@(m[i, j] /. sol), {i, n}, {j, n}], 
  Table[f[i][e0], {i, n}]] // MatrixForm

Where f[i][e[[i]]] is the force acting on the particle i obtained by using the corresponding velocities in the Stokes equation from the question.

In the general case where the equations are not linear on the particles' velocities one needs to make a minimization routine that finds the configuration that makes the force zero, say for example with a steepest descent method. This would be even more computationally expensive.


As requested by user21, here is my code.

First, I need to explain that I'm solving the phoretic interaction between colloids. The particles are platinum spheres in a solution of peroxide water. The surface of the particles reacts with the peroxide water and produces hydrogen. This generates a concentration field for the hydrogen having as sources the particles. Gradients on the fields will induce forces on the particles such that there is a solute-mediated interaction between them. Here's the paper you want to read if you are interested: review of phoresis

So in my case I first compute the concentration field due to the particles. With the gradient of the concentration I can compute the slip velocity of the fluid on the surface of the particles, and when a velocity is imposed on the particles, I can solve stokes equations for the system. What I'm really interested on is how the phoretic force depends on the distance between the particles. Currently, people assume this is a pair-wise interaction that decays as $1/r^2$ in 3D, however it's slightly more complicated than that.

Let's take the case of the two equal particles so the minimization problem is unidimensional. Without loss of generality let's assume the particles repel each other. A diagram of the situation may be useful. enter image description here

If the particles are hold at null velocity they feel a repulsive force, black arrows on top figure. If they move with a certain velocity they feel a drag opposite to their movement, middle panel. For a certain velocity $U_0$ the force on the particles is zero and this move apart from each other, lower panel.

enter image description here

This is the force for the $U=0$ (blue) and $U=1$ (red) as a function of the distance between the particles. The force decays in my case because the gradients of the concentration field decay with distance (on the very long range the solution to the concentration equation is almost the sum of two independent sources). Furthermore, the red curve changes it sign since for any positive velocity there is a distance such that this velocity is larger than the critical velocity. Or said another way, $U=1$ is the critical velocity (the velocity for which the force on the particle is zero) for $\delta r = 4$.

I hope you are not lost. I make all this introduction because this is not the general kind of boundary condition that one solves, and I want to give the rationale behind it.

So here the code to compute force1 as a function of r. You also need to compute force0 by replacing U = {0,0}

<< NDSolve`FEM`
L = 40;
force1 = Range@20;
ClearAll[concentration, vSlip]
 x0 = 1.5 + i/5;

 cords = {{-x0, 0}, {x0, 0}};

 δΩ =  
  Apply[Or, Norm[{x, y} - #]^2 == 1.0 & /@ cords];
 Ω = Apply[And, Norm[{x, y} - #]^2 > 1 & /@ cords];

 (*Create the mesh and visualise*)

 (*The same mesh will be used to both solve the diffussion equation \
and the stokes eq. *)

 r = RegionDifference[Rectangle[{-L, -L}, {L, L}], 
   RegionUnion[Disk[#, 1] & /@ cords]]; 
 mesh2 = ToElementMesh[ r, "MaxBoundaryCellMeasure" -> .1, 
   "MaxCellMeasure" -> 1.5, MeshQualityGoal -> "Maximal", 
   MeshRefinementFunction -> 
    Function[{vertices, area}, 
     area > 0.0125 (0.1 + 
         Min[Norm[Abs@Mean[vertices] - {x0, 0}], 

 (*First we solve the diffusion equation in the system*)

 concentration = NDSolveValue[{\!\(
\*SubsuperscriptBox[\(∇\), \({x, y}\), \(2\)]\(u[x, y]\)\) == 
     NeumannValue[1, (Abs[x] - x0)^2 + y^2 == 1], 
    DirichletCondition[u[x, y] == 0, Abs[x] >= L || Abs[y] >= L]}, 
   u, {x, y} ∈ mesh2, Method -> "FiniteElement"];

 (*once the perturbation due to the particle is taken into account, \
one can use this to compute vslip*)

 gradC[x_, y_] := Evaluate[D[concentration[x, y], {{x, y}}]]; 

 (*Second, define stoke's problem*)

 eqn = {
   Laplacian[u[x, y], {x, y}] - D[p[x, y], x] == 0,
   Laplacian[v[x, y], {x, y}] - D[p[x, y], y] == 0,
   D[u[x, y], x] + D[v[x, y], y] == 0};

 (*with null boundary conditions in "infinite"*)

 bcs = {u[x, -L] == 0, u[x, L] == 0,
   u[-L, y] == 0, u[L, y] == 0,
   v[x, -L] == 0, v[x, L] == 0,
   v[-L, y] == 0, v[L, y] == 0,
   p[L, y] == 0, p[-L, y] == 0};

 (*deine the slip velocity*)

 vSlip[x_, y_] := 
  Module[{ϕ = ArcTan[If[x > 0, x - x0, x + x0], y]}, 
     y] - {Cos[ϕ], 
      Sin[ϕ]} ({Cos[ϕ], Sin[ϕ]}.gradC[x, y])];

 (*To use the same notation as the paper, we call the velocity of the \
particle U, this is our variable to minimise later. The ansatz is U = \
\*SubscriptBox[\(∇\), \(∞\)]C\)*)

 U = {1., 0.};

 (*solve diffusion equation wtih fixed boundary conditions and \
evaluate vslip as the gradient on the surface of the particle*)

 pde = {eqn, bcs, 
   DirichletCondition[u[x, y] == vSlip[x, y][[1]] + U[[1]](*Sign[
    x]*), δΩ], 
    v[x, y] == vSlip[x, y][[2]] + U[[2]], δΩ]};

 {xVel, yVel, pressure} = 
  NDSolveValue[pde, {u, v, p}, {x, y} ∈ mesh2, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

 Subscript[σ, xx][x_, y_] = Evaluate[2 D[xVel[x, y], {x}]];
 Subscript[σ, xy][x_, y_] = 
  Evaluate[D[xVel[x, y], {y}] + D[yVel[x, y], {x}]];

 force1[[i]] = 
  Quiet@NIntegrate[(-pressure[x0 + Cos[ϕ], Sin[ϕ]] + 
        Subscript[σ, xx][x0 + Cos[ϕ], 
         Sin[ϕ]]) Cos[ϕ] + 
     Subscript[σ, xy][x0 + Cos[ϕ], 
       Sin[ϕ]] Sin[ϕ], {ϕ, 0, 2 π}];

 , {i, 1}]

When you have both force1 and force0 you can obtain by u = force0/(force1 - force0) and it looks something like:enter image description here

In blue force0, red is the critical velocity that makes $f=0$ on the particle and green is the currently used model. Force and velocity are in the same plot since we are in the low reynolds regime.

(sorry for the plots with different axis but this is still work in progress.)

  • $\begingroup$ Is there a chance that you post your code and how it works together with the Stokes equation? If I have a solution that you think is correct I could try write a FEM boundary condition for that. $\endgroup$
    – user21
    Commented May 3, 2016 at 11:37
  • $\begingroup$ there it is, I hope it's clear enough. If not please just ask. $\endgroup$ Commented May 3, 2016 at 14:46
  • $\begingroup$ Great, thanks! I'll have a look and see what can be done. Please be patient I'll not be able to do it right now. I need to finish something else first. $\endgroup$
    – user21
    Commented May 3, 2016 at 21:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.