A sphere heads east from Los Angeles ( centre of hydrodynamic 'mass' problem)

[vputz]

A sphere heads east from Los Angeles... ("centre of hydrodynamic 'mass'" problem)

Homework Statement

This is actually from research, but it's one of those "this should be an easy problem and I'm wasting too much time making it too difficult" problems, so I thought it fit better under math/calculus--at its heart it is a math-y problem rather than a physics-y one.

The problem is this: I'm doing numeric simulations of small hydrodynamic swimmers composed of pairs of spheres (dumbbells); for notation, each sphere has position [tex]x_1[/tex] and [tex]x_2[/tex], and radius [tex]a_1[/tex] and [tex]a_2[/tex]. Each dumbbell has a time-varying force applied in opposite directions to each sphere.

So as you can imagine, you basically have something that expands and contracts in a periodic fashion. We model the effect of force on a sphere by Stokes drag, ie

[tex]v=\frac{1}{6\pi\eta a}F[/tex]

(where [tex]\eta[/tex] is viscosity and F is the force on the sphere) and hydrodynamic interactions between spheres with the Oseen tensor

[tex]\frac{1}{8\pi\eta r}(\delta_{\alpha\beta} + \hat{r_\alpha}\hat{r_\beta})[/tex]

(actually we use the Rotne-Prager tensor, but an Oseen approximation may be good enough). In one dimensions, it's fairly simple (see sec 2).

So finally, here's the problem: What is the "position" of a dumbbell? I want some function of the position of the two spheres ([tex]x_1[/tex] and [tex]x_2[/tex]) that does not vary over the duration of a cycle. For symmetric dumbbells, the forces are equal and we can say

[tex]x_{cm} = \frac{x_1+x_2}{2}[/tex]

...but this naive version fails with asymmetric spheres and the quantity varies over time. A naive update based on the Stokes force is

[tex]x_{cm} = \frac{ a_1 x_1 + a_2 x_2 }{a_1+a_2}[/tex]

The amplitude of this varies over time as well, according to simulation, although to much less degree. In our simulations, we've taken to using the naive version and just measuring at the beginning of each period, which gives acceptable results that agree with analytics, and no one else seems interested in this question, but now it is bugging me.

What is some equation for [tex]x_{cm}[/tex] that would not vary over time? If you stretched a line from [tex]x_1[/tex] to [tex]x_2[/tex], what point on that line would not move?

Homework Equations

We have a coupled system of ODEs. In one dimension, the stokes force and Oseen interactions reduce to

[tex]\frac{d}{dt}x_1 = \frac{F}{6\pi\eta a_1} - \frac{F}{4\pi\eta(x_2-x_1)}[/tex]
[tex]\frac{d}{dt}x_2 = \frac{F}{4\pi\eta(x_2-x_1)} - \frac{F}{6\pi \eta a_2}[/tex]

Find some point [tex]x_{cm}(x_1, x_2)[/tex] which is invariant with respect to time, representing the "position" of the dumbbell.

The Attempt at a Solution

This seems like the sort of thing which would be quite easy, but for some reason it's vexing me. I'm having difficulty framing the question; I feel like it might be easy once I frame it properly.

One way to think about it is "If the force was always contractile (presuming that [tex]x_2 > x_1[/tex]), at what point [tex]x_{cm}[/tex] would the spheres come into contact?" -- I think that's a valid idea, because that point should be the same no matter where in the sinusoidal cycle the dumbbell was. Unfortunately I haven't made much progress here. I've also thought about, as mentioned above, imagining a line stretched between the spheres and finding a point on it which doesn't move, but that hasn't led me far either.

Another way might be to imagine some function [tex]x_{cm}(x_1, x_2)[/tex] such that [tex]x_{cm}(x_1(t), x_2(t)) = x_{cm}(x_1(t+dt), x_2(t+dt) )[/tex], but since I'm trying to guess the form of [tex]x_{cm}[/tex], that hasn't helped much.

I thought comparing the derivatives, might be helpful, in that

[tex]\frac{ \frac{d}{dt}x_1 }{ \frac{d}{dt}x_2 } = - \frac{ a_2( 2( x_2 - x_1 ) - 3 a_1 ) }{ a_1( 2( x_2 - x_1 ) - 3 a_2 ) }[/tex]

That's a nice-looking useful ratio, but I can't seem to do much useful with it.

I'm also more interested in the general approach, since the off-diagonal terms will change slightly if we use the RPY tensor instead of Oseen... and the equation also must extend to 3d coordinates.

This is really rather embarassing since moving to computers from my analytic undergraduate days... I can simulate hundreds of these guys in parallel using CUDA and applied stochastic noise, but I can't even frame the above correctly. What horribly obvious thing am I missing?

 

[vputz]

My latest idea, which may work:

Define some point on the line between x1 and x2:

[tex]x_{cm} = x_1 + y( x_2 - x_1 )[/tex]

Then say that this point doesn't move, set d/dt of it to zero, solve for y; since we know d/dt of [tex]x_1[/tex] and [tex]x_2[/tex], this may work (I'm imagining that y does not change with time, since it's a proportion, so assuming that dy/dt=0).

 

[vputz]

...sadly, the above works brilliantly for simple stokes velocity when everything is nice and linear, but the introduction of the cross-terms means this doesn't work (so the idea of a "stationary point on the dumbbell" seems invalid, or that it moves over time so dy/dt above is not zero). Hmm. So I may be back to the idea of "if it contracts until the spheres are touching, that is the "position" of the dumbbell... no progress yet.

 

[vputz]

Okay, I thought about it a bit more, and there are a couple ways forward.

Thinking abstractly, we have something shaped like this:

[tex]\frac{dx_1}{dt} = f(t)\left( g(a_1)-h(x_2-x_1)\right)[/tex]
[tex]\frac{dx_2}{dt} = -f(t)\left( g(a_2)-h(x_2-x_1)\right)[/tex]

In our case, [tex]g(a)[/tex] is constant over time, while [tex]h[/tex] is difficult since it varies with [tex]x_1[/tex] and [tex]x_2[/tex].

What we want is some function of [tex]x_1[/tex] and [tex]x_2[/tex] that is constant in time which we can use as a position.

One way to proceed, to get rid of the tricky [tex]h[/tex] dependence, is to add the derivatives, which would give us something like

[tex]\frac{d}{dt}(x_1+x_2) = f(t)\left( g(a_1)-g(a_2) \right)[/tex]

Since [tex]g(a)[/tex] is constant in time, we can eventually get to

[tex]\frac{x_1+x_2}{g(a_1)-g(a_2)} - \int f(t)dt = C[/tex]

... so long as we exactly know the nature of [tex]f(t)[/tex] in the simulation. Well, we do know; it's just a sinusoidal applied force. So this could suit my needs.

BUT IT'S NOT ELEGANT. And that bothers me. As long as the force is paired as it is above (equal and opposite) it shouldn't matter what the force is. So we'd like to eliminate [tex]f(t)[/tex] in our expression, and the best way I can think of doing that is by dividing and eliminating [tex]t[/tex] altogether:

[tex]\frac{dx_1}{dx_2} = - \frac{g(a_1)-h(x_2-x_1)}{g(a_2)-h(x_2-x_1)}[/tex]

...which gets us the equation in the first post ("looks like a useful ratio") for a simple version of [tex]g[/tex] and [tex]h[/tex]. But I can't figure out how to integrate this (incidentally, if [tex]h=0[/tex] or [tex]a_1 = a_2[/tex], the above is awfully easy indeed). It looks as if it should be easy, but I wind up flailing. Worse, the next correction would make [tex]h[/tex] quadratic in [tex]x_2-x_1[/tex].

Grr. Any suggestions?