[aspect-devel] Thermodynamic consistency of Aspect's temperature and momentum equations

[Ian Rose]

Hi Aspect folks,

I was working through some tests with Aspect and came across what I believe
is an inconsistency in the governing equations.

For incompressible Boussinesq flow, the global viscous dissipation should
exactly cancel the global adiabatic heating.  This can be seen by
multiplying the momentum equation by velocity and integrating over the
domain.

As it stands in assembly.cc, the formula used for calculating adiabatic
heating is different from that you would get by integrating the momentum
equation.  I wrote a simple postprocessor that compares the two integrated
quantities which I am attaching.  The difference is quite a lot for the
current formula.

Put another way, this is the formula that is currently used:

   Q_a = ( velocity * gravity ) * alpha * density * temperature

The density at this point however, has already been adjusted for
temperature, so we are in effect double counting the thermal expansion.
 Instead, I believe it should be

  Q_a = ( velocity * gravity ) * ( density - reference_density )


The compressible case, too, should require some thought, though I have not
gone through the paces there.

Thoughts?

Best,
Ian

PS, for some details on the derivations, I refer you to Leng and Zhong
(2008)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://geodynamics.org/pipermail/aspect-devel/attachments/20130206/9bdb7116/attachment.htm 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: energetics.cc
Type: application/octet-stream
Size: 7949 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/aspect-devel/attachments/20130206/9bdb7116/attachment.obj 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: energetics.h
Type: text/x-chdr
Size: 1486 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/aspect-devel/attachments/20130206/9bdb7116/attachment.h