[Nek5000-users] Mantle convection (infinite Prandtl, variable viscosity) with Nek5000

[nek5000-users at lists.mcs.anl.gov]

Hi Neks,

I'm trying to use Nek to simulate convection in planetary mantles (infinite Prandtl-number, variable viscosity...) in the way Jan has requested a few months ago:

>> Hi all,
>> I was wondering if it is possible to use Nek5000 for simulations in 
>> the infinite Prandtlnumber regime (e.g. planetary mantles). This 
>> would require to solve the NS Stokes equation without the time 
>> derivative and advection term and the heat equation for T:
>>
>> 1) Solve visc * laplace u - grad p + buoyancy = 0 and div u = 0, 
>> where bouyancy depends on a given T field and a Rayleigh number
>> 2) Update T by solving dT/dt + u grad T = diff * laplace T (+ other 
>> sources)
>> 3) go back to 1 with the new T field, repeat
>>
>> More realistic models would then use a varying viscosity.
>> I looked at the steady state example (kov_st_state), but this is 
>> missing the time depencence via the Temperature. Any hints on the 
>> possibility would be much appreciated.
>>
>> Thanks,
>> Jan

Initially I tried the whole thing with a simple 2D rectangular geometry and a fluid with temperature-dependent viscosity ([udiff=param(2)*exp(-log(param(70))*temp)] were param(70) is the viscosity contrast between temp=1 and temp=0) cooled from above (T=0) and heated from below (T=1) .

At first I tried the way Paul suggested:

>> Hi Jan,

>> I just saw this and the follow-on posts.  (Apologies - have been on the road for
>> several weeks.)
>> 
>> I've not tried this, but it appears to me that you simply will set IFNAV = F in the .rea 
>> file and can set param(1) [ rho ] to be a very small number.  Param(2) is the dynamic
>> viscosity, so you can set that independently of param(1).   This will have the affect
>> of forcing the velocity to relax very quickly (its inertia is negligible compared to the
>> other terms in the momentum equation).  So, you could conceivably set param(1)
>> to be 1.e-30, say.
>> 
>> Note that, on the same line where IFNAV is set, you would set IFADV to be "T" for
>> the second field.  It should look something like:
>> 
>> F T F F F F F F F F F IFNAV & IFADVC (convection in P.S. fields)

>> Also, you might find it of value in this scenario to try the Pn-Pn-2 formulation with
>> plan1 ...  In this case, you would indeed want to edit the source --- specifically:

>> change the call to plan3() to call plan1()  in drive2.f ...    As I say, I've not experimented
>> with this flow regime so am only making suggestions to try, rather than firm
>> recommendations based on experience.

>> Paul

but I got the problem that Nek worked well only for one condition, either for an infinite Prandtlnumber (i.e. density=1.e-30) or variable viscosity. 
As I wanted Nek to solve both, I got some strange results... for a sufficient resolution very thin boundary layers evolve (they should be much thicker with for instance Rayleigh numbers like 10^4 and a viscosity contrast of 1000) and the boundary layers had a rough structure, but they should be quite smooth.
For lower resolution I got no or completely wrong results, either the temperature was oscillating very strong (with values higher than the maximum temperature of 1) or Nek crashed at all.

Maybe someone has an idea how to fix it?!


However, after that I also tried the way lailai suggested:

> Hello, Jan,
>
> We have in fact solved a similar problem as you face now. In our case, 
> we have steady Stokes equation at each time step,
> in the absence of time derivative and convection term. However we have 
> time-dependent volume forcing and boundary
> conditions, as we solve for fluid-structure interactions in the Stokes 
> regime.
>
> You only need to change the source code for one thing. By changing the 
> flags in the .rea file, you can easily turn
> off the two terms. However, the code only runs for 1 time step even if 
> you specify for example 10 total time steps, which makes sense
> as NEK realizes this is a steady problem. NEK5000 does not see that 
> the next time step, BCs and forcing will be changed, and so
> as the flow.
>
> So I remember I only commented two lines of the source code to achieve 
> this, in the file 'connect2.f', subroutine 'rdparam',
> this is what I have done (i am using a quite old version of the NEK, 
> so...)
>
>       IF (.NOT.IFTRAN) THEN
> c         PARAM(11) = 1.0 !lailai comment for steady stokes
> c         PARAM(12) = 1.0 !lailai comment for steady stokes
>          PARAM(19) = 0.0
>       ENDIF
>
> Of course you need to set the flags right in .rea file, which you can 
> find the details on the webpage. You also need to
> solve for heat equations to get T. Then in the userchk subroutine, you 
> calculate the T-dependent buoyancy forcing at each time step.
> Solve the flow with new forcing terms. This should be all manageable.
>
> cheers and good luck,
>
> lailai

where I've included Jan's advice as well:

> Hi,
> thanks for the advice. I found that one also has to comment the line
> IF (.NOT.IFTRAN) NSTEPS=1
> in setvar in drive2.f and use a fixed timestep so that the timeloop does 
> not stop too soon.
> 
> Thanks again,
> Jan

Using this method I got more satisfying results, at least for the tested parameters (Ra=10^4, viscosity contrast of 1000), but to achieve a statistically steady state it took a very long time (non-dimensional) of t>1 (independent of the maximum stepsize dt) and therefore very much iteration steps (~10^6) and CPU-time. For comparison other mantle convection codes need a non-dimensional time of t<0.1 and less than 10000 steps (with a comparable maximum stepsize). 
Unfortunately I've no idea why it takes that much time using Nek...

I would be vary grateful for any hints or suggestions to make at least one method work (either Paul's or lailai's way).

Thanks,
Philipp

----------------------
Philipp Hellenkamp
Institut für Geophysik
Corrensstr. 24
D-48149 Münster