[aspect-devel] cylindrical rigid body rotation

[Wolfgang Bangerth]

On 02/27/2014 10:31 AM, Timo Heister wrote:
>> Citcom has the same phenomenon, and they just subtract off rigid rotations
>> at every timestep.  I do have some code that addresses this which is not in
>> the mainline.
>
> The method used in citcom is described in Zhong et al 2008 (G^3), section 2.4.1.
>
> I wonder if there is some need for eliminating the rotational modes in
> a sublayer instead of the whole geometry (as described at the end of
> that section).

If the goal is simply to remove the rotation from the solution (i.e., to 
eliminate the kernel of the operator), then it's enough to run the 
computation of the rotational mode over the entire domain. I think 
Shijie only presents the formula as a general case, but I don't read it 
as saying that it should be applied to anything other than the whole mantle.

Formulas 24 and 25 in the paper correspond to the case of constant 
density but we should implement the more general case. The rotational 
moment of the *corrected* velocity should equal zero, i.e., we need an 
omega so that

   \int  \rho(x)  x \times [u(x) - x \times omega]  =  0

i.e., an omega so that

   \int \rho(x)  x \times [x \times \omega] = \int \rho(x)  x \times u(x)

This gives

   \omega = I^{-1} \int \rho(x)  x \times u(x)

where the matrix I is defined by

   I_{im} = \int \rho(x)  eps_{ijk} x_j eps_{klm} x_l

and eps is the totally antisymmetric tensor. (I'm sure there is a 
simpler form of the integrand with the two eps, but I'm not seeing it 
right now.) Anyway, this should be simple enough to compute, and then 
you can use formula 26 to correct the velocity.

Has someone already implemented something in this form?

Cheers
  W.


-- 
------------------------------------------------------------------------
Wolfgang Bangerth               email:            bangerth at math.tamu.edu
                                 www: http://www.math.tamu.edu/~bangerth/