[aspect-devel] Thermodynamic consistency of Aspect's temperature and momentum equations
Rene Gassmoeller
rengas at gfz-potsdam.de
Wed Feb 13 04:27:32 PST 2013
Hi all,
Juliane and myself had a quite lengthy discussion about this issue (it
started yesterday and it is actually not finished yet ;-)) but I just
wanted you to know what would be our opinion up to now (to save double
effort in thinking).
We think Ian had an interesting point here, although Thomas is right
about the correctness of the currently used term in a way as well. Ian's
term is derived from the momentum equation and exactly balances out the
viscous dissipation (since it is derived that way). However, there is no
guarantee that this term actually is(!) the adiabatic heating in all
material models / approximations, except that it seems to be the right
term for the boussinesq / anelastic liquid approximation with a simple
material model.
The term that is currently used in Aspect on the other hand is derived
from the actual thermodynamic equations (dS/dt = ... the root of our
equation of conservation of energy as stated in "Mantle convection in
the Earth and Planets") and is therefore consistent, but it uses at
least one approximation, namely that the pressure is only lithostatic
and there is no dynamic pressure. The right term (in a thermodynamic
sense) would not be:
Q_a = ( velocity * gravity ) * alpha * density *
temperature (1)
but rather:
Q_a = ( velocity * pressure_gradient ) * alpha *
temperature (2)
The gravity and density in the currently used term (1) just came out of
the assumption that p = rho * g * h, and by the way, this justifies the
use of the real rho instead of the reference_rho in the equation (1)
although it does not seem to fit with the viscous heating. The use of
the above mentioned formula (2) reduces the error in a test case
(attached) from around 3% to 0.2% ... that is not perfect considering
the fit of Ian's term of 1e-13 or so. Still I would currently rather use
term (2) since it also makes us independent of the choice of a right
reference density and it is also consistent in the compressible case
(because it makes no assumptions on the material model at all, it just
needs the real alpha und the local solution variables). Additionally
with this term we also capture the heating of material that moves
laterally along a pressure gradient and therefore also cools/heats up
adiabatically.
Nevertheless, the misfit of 0.2% seems too large to be satisfied. I
first thought that the inconsistency between the definition of alpha = -
(1 / rho) * drho / dT = const and the definition of rho = rho_ref (1 -
alpha * (T - T_ref)) leads to the misfit and tried a definition of rho =
rho_ref * exp (- alpha * (T - T_ref)), which is the solution of the
definition of alpha above. But this does not seem to make a large
difference in the error. So currently we are thinking in two ways:
1. Find the reason for the misfit between thermodynamically derived
adiabatic heating term and viscous dissipation that is computed out of
the stokes velocity field.
2. If we can derive that Ian's term actually is the adiabatic heating
term in a general sense (perhaps with an additional or without an
currently used approximation in term (2)), then of course that term
would be more precise, although then we end up asking how to use this in
a compressible model.
Currently, I am a bit sceptical about the possibility of point 2.
therefore I am thinking of point 1.. The next part is speculation but I
find it a bit intriguing that Leng and Zhong (2008) also end up with a
difference of the order of 0.2 % although they do not consider the
formula (2) but rather Ian's term extended by an compressible
correction. Could it be that the remaining 0.2% originate from an
approximation that is not bound to the formulation of the adiabatic
heating term but rather something influencing the velocity field (and
therefore the viscous dissipation)?
By the way during this discussion we came across the formulation of
compressible stokes flow in aspect and think there is an inconsistency
as well. In the right-hand side of the stokes equation we have the term
rho*g as well as an additional compressible term. Considering a
compressible term in the compressible case seems reasonable if the first
term just covers the temperature dependency of rho, but since in Aspect
the material model takes care of all the influences on the density, this
term should be moved to any material model that wants to be capable of
doing compressible convection. For example a material model could take
the density values directly from a thermodynamic calculation rather than
a simplified equation of state. In that case the current formulation
would double-count the pressure effect (without any chance for the user
to know this, since he thinks the material model need to take care of
all the influences).
We are continuing working on this, but please let us know, if you have
an opinion or suggestion on this (and also in case something in this is
wrong ;-)),
Cheers,
Rene
--
Rene Gassmoeller
2.5/Geodynamic Modelling
Tel.: +49 (0)331/288-28744
Fax: +49 (0)331/288-1938
Email: rengas at gfz-potsdam.de
___________________________________
Helmholtz Centre Potsdam
German Research Centre for Geosciences GFZ
Telegrafenberg, 14473 Potsdam
On 02/13/2013 01:21 AM, Magali Billen wrote:
> Hello All,
>
> Each of the different approximations of the equations have specific
> sets of terms that drop in or out together,
> so you need to be careful adding back in just one term or another
> without taking into account which approximation
> has leads to that term being assumed small. When people state
> Boussinesq approximation in mantle convection
> calculations, this has historically meant an approximation that does
> not include adiabatic heating. Also, I agree with
> Thomas that it doesn't really make sense to have such a term in an
> incompressible convection case, since without compression there is no
> physical cause for an adiabatic gradient.
>
> Speaking for myself, I don't have the specific approximations
> memorized, however, I've found Chapter 6 of
> "Mantle Convection in the Earth and Planets" by Schubert, Turcotte and
> Olson, very helpful when trying to understand
> the origin and loss of terms for boussinesq, extended-boussinesq and
> TALA approximations - it goes through
> each of these approximation in detail and explains the specific
> assumptions for each approximation and which
> terms drop out. I worked through all these once about a year ago and I
> think it might help with this particular question.
>
> Magali
>
> On Feb 12, 2013, at 4:09 PM, Ian Rose wrote:
>
>> Hmm, I am not sure I agree. Di is frequently assumed to be zero in
>> mantle convection problems, but that is not a result of the
>> Boussinesq approximation. That is to say, the "work done by gravity"
>> term in the kinetic energy equation arises just fine with Boussinesq
>> (there are just more terms that come from the div velocity terms in
>> the compressible case).
>>
>> Even though this term (along with viscous dissipation) are likely to
>> be smallish, I see no reason not to allow them to be turned on and
>> off with flags as they are now. But if it is turned on, it should be
>> consistent with what you get from integrating the momentum equation.
>>
>> Cheers,
>> Ian
>>
>>
>>
>> On Tue, Feb 12, 2013 at 3:28 AM, Thomas Geenen <geenen at gmail.com
>> <mailto:geenen at gmail.com>> wrote:
>>
>> he Timo,
>>
>> there is no such thing as adiabatic heating in the incompressible
>> Boussinesq case Di (alpha*g/cp) is assumed zero .
>> for extended Boussinesq there should also be no problem since
>> there is no density in the net adiabatic heating term.
>>
>> setting thermal diffusion, viscous dissipation and internal
>> heating to zero (dS/dt=0) we end up with
>> rhocp(dT/dt) - alphaTdP/dt=0
>> or
>> rho*cp*(dT/dt) - alpha*rho*g*u_r*T=0
>>
>> this will give for an adiabatic temperature profile
>> T(r) = T_0*exp(alpha*g*r/cp)
>>
>> iow the density does not play a role since its devided out of the
>> equation.
>>
>> this also holds for the compressible case i would say.
>>
>> cheers
>> Thomas
>>
>>
>>
>> On Tue, Feb 12, 2013 at 5:57 AM, Timo Heister
>> <heister at math.tamu.edu <mailto:heister at math.tamu.edu>> wrote:
>>
>> Hey everyone,
>>
>> Ian approached me about this and I asked him to write it down
>> here.
>> Does anyone have any feedback about this, especially
>> (assuming this is
>> correct), what to do in the compressible case?
>>
>> On Wed, Feb 6, 2013 at 6:33 PM, Ian Rose
>> <ian.rose at berkeley.edu <mailto:ian.rose at berkeley.edu>> wrote:
>> > 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)
>> >
>> >
>> > _______________________________________________
>> > Aspect-devel mailing list
>> > Aspect-devel at geodynamics.org
>> <mailto:Aspect-devel at geodynamics.org>
>> > http://geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>>
>>
>>
>> --
>> Timo Heister
>> http://www.math.tamu.edu/~heister/
>> <http://www.math.tamu.edu/%7Eheister/>
>> _______________________________________________
>> Aspect-devel mailing list
>> Aspect-devel at geodynamics.org
>> <mailto:Aspect-devel at geodynamics.org>
>> http://geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>>
>>
>>
>> _______________________________________________
>> Aspect-devel mailing list
>> Aspect-devel at geodynamics.org <mailto:Aspect-devel at geodynamics.org>
>> http://geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>>
>>
>> _______________________________________________
>> Aspect-devel mailing list
>> Aspect-devel at geodynamics.org <mailto:Aspect-devel at geodynamics.org>
>> http://geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>
> --------------------------------------------------
> Associate Professor & Chancellor Fellow, U.C. Davis
> Chair, Geology Graduate Program
> Department of Geology & KeckCAVES
> 2129 Earth & Physical Sciences Bldg
> Davis, CA 95616
>
> E-mail: mibillen at ucdavis.edu <mailto:mibillen at ucdavis.edu>
> Phone: no extension, please e-mail
> --------------------------------------------------
>
>
>
>
>
>
>
>
>
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://geodynamics.org/pipermail/aspect-devel/attachments/20130213/170d1add/attachment-0001.htm
-------------- next part --------------
# At the top, we define the number of space dimensions we would like to
# work in:
set Dimension = 2
# There are several global variables that have to do with what
# time system we want to work in and what the end time is. We
# also designate an output directory.
set Use years in output instead of seconds = false
set End time = 1e17
set Output directory = output
# Then there are variables that describe the tolerance of
# the linear solver as well as how the pressure should
# be normalized. Here, we choose a zero average pressure
# at the surface of the domain (for the current geometry, the
# surface is defined as the top boundary).
set Linear solver tolerance = 1e-15
set Temperature solver tolerance = 1e-15
set Pressure normalization = surface
set Surface pressure = 0
set Adiabatic surface temperature = 0
# Then come a number of sections that deal with the setup
# of the problem to solve. The first one deals with the
# geometry of the domain within which we want to solve.
# The sections that follow all have the same basic setup
# where we select the name of a particular model (here,
# the box geometry) and then, in a further subsection,
# set the parameters that are specific to this particular
# model.
subsection Geometry model
set Model name = box
subsection Box
set X extent = 2890000
set Y extent = 2890000
end
end
# The next section deals with the initial conditions for the
# temperature (there are no initial conditions for the
# velocity variable since the velocity is assumed to always
# be in a static equilibrium with the temperature field).
# There are a number of models with the 'function' model
# a generic one that allows us to enter the actual initial
# conditions in the form of a formula that can contain
# constants. We choose a linear temperature profile that
# matches the boundary conditions defined below plus
# a small perturbation:
subsection Initial conditions
set Model name = function
subsection Function
set Variable names = x,z
set Function constants = p=100, L=2890000, pi=3.1415926536, k=1
set Function expression = 1613.0 + p*cos(k*pi*x/L)*sin(pi*z/L)
end
end
# Then follows a section that describes the boundary conditions
# for the temperature. The model we choose is called 'box' and
# allows to set a constant temperature on each of the four sides
# of the box geometry. In our case, we choose something that is
# heated from below and cooled from above. (As will be seen
# in the next section, the actual temperature prescribed here
# at the left and right does not matter.)
subsection Boundary temperature model
set Model name = box
subsection Box
set Bottom temperature = 2600
set Left temperature = 0
set Right temperature = 0
set Top temperature = 273
end
end
# We then also have to prescribe several other parts of the model
# such as which boundaries actually carry a prescribed boundary
# temperature (as described in the documentation of the `box'
# geometry, boundaries 2 and 3 are the bottom and top boundaries)
# whereas all other parts of the boundary are insulated (i.e.,
# no heat flux through these boundaries; this is also often used
# to specify symmetry boundaries).
subsection Model settings
set Fixed temperature boundary indicators = 2,3
# The next parameters then describe on which parts of the
# boundary we prescribe a zero or nonzero velocity and
# on which parts the flow is allowed to be tangential.
# Here, all four sides of the box allow tangential
# unrestricted flow but with a zero normal component:
set Zero velocity boundary indicators =
set Prescribed velocity boundary indicators =
set Tangential velocity boundary indicators = 0,1,2,3
# The final part of this section describes whether we
# want to include adiabatic heating (from a small
# compressibility of the medium) or from shear friction,
# as well as the rate of internal heating. We do not
# want to use any of these options here:
set Include adiabatic heating = false
set Include shear heating = false
set Radiogenic heating rate = 0
end
# The following two sections describe first the
# direction (vertical) and magnitude of gravity and the
# material model (i.e., density, viscosity, etc). We have
# discussed the settings used here in the introduction to
# this cookbook in the manual already.
subsection Gravity model
set Model name = vertical
subsection Vertical
set Magnitude = 10 # = Ra / Thermal expansion coefficient
end
end
subsection Material model
set Model name = simple # default:
subsection Simple model
set Reference density = 3300
set Reference specific heat = 1200
set Reference temperature = 1613
set Thermal conductivity = 4
set Thermal expansion coefficient = 3e-5
set Viscosity = 1e22
end
end
# The settings above all pertain to the description of the
# continuous partial differential equations we want to solve.
# The following section deals with the discretization of
# this problem, namely the kind of mesh we want to compute
# on. We here use a globally refined mesh without
# adaptive mesh refinement.
subsection Mesh refinement
set Initial global refinement = 4
set Initial adaptive refinement = 0
set Time steps between mesh refinement = 0
end
# The final part is to specify what ASPECT should do with the
# solution once computed at the end of every time step. The
# process of evaluating the solution is called `postprocessing'
# and we choose to compute velocity and temperature statistics,
# statistics about the heat flux through the boundaries of the
# domain, and to generate graphical output files for later
# visualization. These output files are created every time
# a time step crosses time points separated by 0.01. Given
# our start time (zero) and final time (0.5) this means that
# we will obtain 50 output files.
subsection Postprocess
set List of postprocessors = velocity statistics, temperature statistics, heat flux statistics, visualization, energetics
subsection Visualization
set Time between graphical output = 0.01
set List of output variables = nonadiabatic pressure, density
end
end
# things to play with:
# - discretization
# - different Ra number
#subsection Discretization
# set Stokes velocity polynomial degree = 2
# set Temperature polynomial degree = 2
#end
-------------- next part --------------
A non-text attachment was scrubbed...
Name: energetics.h
Type: text/x-chdr
Size: 1487 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/aspect-devel/attachments/20130213/170d1add/attachment-0001.h
-------------- next part --------------
A non-text attachment was scrubbed...
Name: energetics.cc
Type: text/x-c++src
Size: 8341 bytes
Desc: not available
Url : http://geodynamics.org/pipermail/aspect-devel/attachments/20130213/170d1add/attachment-0001.cc
More information about the Aspect-devel
mailing list