Solving Stokes equations for zero force boundary condition

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[ And@Flatten@ 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.

Addendum:

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.

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.

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] Table[ 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}], Norm[Mean[vertices]]])]]; (*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]}, gradC[x, 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]*), δΩ], DirichletCondition[ 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:

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.)