Computational Reality X
Thermodynamics
B.Emek Abali @ LKM  TU Berlin
Abstract
After solving the coupled thermoelastic equations in Comp.Real.IX,it is time to give a real explanation
for the internal energy density.This is actually a branch of thermodynamics,though this question never
gets its full answer in a thermodynamics'class.Here we will show one possibility of deriving the energy
formulation which is based on the phenomenological thermodynamics.There are many dierent and even
really sophisticated methods,however,froma pragmatic point of view the phenomenological way is the most
suitable one for applied mechanics.
1 Motivation
We have seen three balance equations so far,local forms of the balance of mass,balance of linear momentum
and balance of internal energy,in the actual frame expressed in an oblique coordinate system x
i
:
+
@v
i
@x
i
= 0;(v
i
)
@
ji
@x
j
f
i
= 0;u
+
@q
i
@x
i
r =
ij
@v
j
@x
i
:
(1)
The rst two has zero sources (right sides) and the internal energy has a production term (nonzero source)
which is always there if the material undergoes a motion.One possible aim is to calculate the mass density,
the velocities and the specic internal energy (;v
i
;u) by satisfying above equations.In Comp.Real.VI & VII
instead of the pressure p has been computed.Similarly in Comp.Real.IX instead of v
i
the displacements,u
i
,
and instead of u the temperature,T,has been obtained.Therefore,the rst declaration in a thermodynamic
process is the space
1
and time,and then the second declaration is the set of primitive variables which are func
tions in space x
i
and time t.The goal is to determine the set of primitive variables.Here,we want to calculate
in material frame X
i
,expressed in Cartesian coordinate system,the displacements u
i
(or strains"
ij
) and tem
perature since it is known that the mass density does not change from
0
in the material frame.We can solve the
balance equations if for
ij
,q
i
and u some constitutive equations are given depending on the primitive variables.
We try to explain the free energy in crystallographic materials and how to measure it,in principle.The rst
assumption is that the total specic (per mass) energy,e,is additive so that the specic internal energy,u,
and specic kinetic energy,e
kin.
,are separable and uncoupled.The author believes that the kinetic energy
is perturbing the material in the nonrecoverable way and the change of the internal energy
2
is recoverable.
There may be dierent formulations,however,none of them is wrong because we cannot measure dierent
parts of the energy separately,they are just models to formulate the behavior of the material.We introduce
temperature T(X
i
;t) and displacement u
i
(X
i
;t) as primitive functions in spacetime.The aim is to calculate
them with appropriate constitutive equations and we suppose that they be independent functions.Due to the
objectivity
3
we derive the gradients of displacements,i.e.,the strains,"
ij
.The measure of the strains is not
important,however,it has to be independent of the temperature,e.g.,the linear symmetric strains can be used:
"
ij
= @u
(i
=@X
j)
.The set of fT;"
ij
g is supposed to determine the thermodynamic state.Hence the internal
energy and all the other derived quantities depend on this set,we may span a statespace and plot
4
the specic
internal energy u = u(T;"
ij
) in that tendimensional space on a perpendicular axis.
1
Material frame for solid bodies,i.e.,Lagrangean frame,or the actual frame for uids,xed or moving Eulerian frame.
2
Under the assumption that no phase changes occur.
3
A displacement is evaluated in a coordinate system,it has other values in another coordinate system,thus it is not objective,
i.e.,it depends on the choice of the coordinate system.However,response of the material has to be objective,hence strain is a
better measure to that,in any coordinate system neighboring particles are neighboring particles!
4
The word plot should not confuse,we (still!) don't know how to plot in so many dimensions,but we may plot in threeaxis one
for temperature,one for an equivalent strain and one for the value of the energy.
1
2 Internal energy
The recoverable internal energy density,i.e.,called the stored energy density,can be changed (variated) by a
variation of the temperature and strains:
dU =
ij
d"
ij
+T dS;(2)
which is fully reversible in the mechanical part
ij
d"
ij
as well as in the thermal part T dS.The symbol d
means the perfect or total dierential and the 1
st
Law of thermodynamics states that the internal energy is a
total dierential,on the next section its necessity will be clear.The newly introduced Cauchy stress tensor
ij
and the entropy density S are the material relations and need to be dened.Since the assumption that
the set fT;"
ij
g describes the state fully,the stress
ij
and the entropy density S should be dened also in that
statespace,thus the change of them read:
d
ij
=
@
ij
@"
kl
T
d"
kl
+
@
ij
@T
"
dT;
dS =
@S
@"
ij
T
d"
ij
+
@S
@T
"
dT;(3)
where we want to explain the meaning of the above notation in detail before we proceed.The partial derivatives
are taken by holding the other arguments xed,thus it is super uous but often it is written for clarity.The
small change d() is simply how we measure it.The stiness tensor C
ijkl
=
@
ij
@"
kl
T
is to be measured on a
constant temperature by variating the strains d"
ij
and recording the stress changes d
ij
.Think about a tensile
loading of a bar where we measure two components (Young's modulus and Poisson's ratio),say,starting from
zero strains up to 0:001 = 0:1%.The tensor of rank four C
ijkl
has 81 coecients in threedimensional space,
so we need couple of more tests than solely a tensile test in one direction.Let us assume that we measure at
T = 300 K and nd all the C
ijkl
components for that temperature,from 0 to 0:1%.The same procedure can
be done at T = 301 K and so on up to the 400 K.Then the stressstrain relation has been determined in the
temperature range of [300 K;400 K] and in the strain range of [0;0:1%].This is of importance,the notation
includes that C
ijkl
may depend on the temperature,it surely does,but each measurement has been done by
variating the strains and holding the temperature xed.This means that the temperature is not variated in the
measurement,i.e.,dT = 0 and the thermal pressure p
ij
=
@
ij
@T
"
drops o when the C
ijkl
are obtained for one
specic temperature.The thermal pressure p
ij
is the pressure when all the strain components of a body is xed,
i.e.,all of its boundaries are clamped,and then the temperature is variated.The material tries to expand or
shrink with the variation of the temperature,thus it applies a pressure due to the thermal change for one specic
strain value.It again implies that the thermal pressure may depend on the strains and we have to perform the
measurement for every strain components"
11
;"
12
;"
13
;"
21
;:::in the range [0;0:1%].We cannot measure the
entropy directly but the heat (the thermal energy) dQ = T dS in a calorimeter can be measured.Hence on
constant strains a variation on the temperature,dT,reads to a change on the heat,dQ,according to the heat
capacity:dQ = C
"
dT which implies C
"
= T
@S
@T
"
.In threedimensional case the 81 components of stiness
tensor,the 9 components of pressure tensor as well as the heat capacity for a range of strains and temperature
can be determined regarding thousands of measurement.At least in principle it is\doable".However,it is not
possible to measure the
@S
@"
ij
T
over the heat energy,since the heat energy is measured always indirectly via
temperature change itself.Therefore the dicult question is,how to obtain the heat developed in the body due
to the recoverable strains.
3 Lagrange multipliers
Before we move on,the method of Lagrange multipliers needs to be outlined.Suppose that there is a function
with a total (or perfect) dierential in a given set of quantities q
= fq
1
;q
2
;:::g:
U =
U(q
):(4)
2
When we search for the minimum(or extremal) value,then it is the value where U gets stationary (or invariant)
under the change (or variation) of these quantities:
U = dU =
@U
@q
q
= 0:(5)
The variation of the quantities dq
is arbitrary,thus the gradient of U in q
has to vanish.The condition as
above dU = 0 is needed to call the dierential perfect or total.By having a total dierential the solution can be
obtained immediately using a rst integral:U =
R
dU +const.up to an integration constant.Hence bringing
the problem into a total dierential form is essential in analysis.Suppose that there are additional conditions
to be fullled,which we can write as residuals:
A =
A(q
) = 0;B =
B(q
) = 0;C =
C(q
) = 0;(6)
thus their variations simply read:A = 0;B = 0;C = 0.The method of Lagrange is to combine them all
by utilizing some constants (constant in q
) as:
U AB C = 0:(7)
The minus between is just a fact of a generalization by Pfaffian form which is not important at the moment.
Though the problem is correct,we have lost the perfect form,however,by rewriting the latter equation:
@U
@q
@A
@q
@B
@q
@C
@q
q
= 0 (8)
and rst by recalling that ;; = const.j
q
and second introducing a new function:
= U AB C;
@
@q
=
@U
@q
@A
@q
@B
@q
@C
@q
;(9)
we get again a perfect or Pfaffian form:
@
@q
q
= 0 = d :(10)
Since the q
is arbitrary,the gradient of should vanish.The newly introduced potential has again a total
diferential and is dened in the same space q
as U.Hence the goal of minimizing U in q
is exactly the same
as trying to minimize the in the same space.
4 Free energy
Now we want to apply the same procedure to the internal energy density U.The notation above is not by
chance.Suppose that there exists a potential,the Helmholtz free energy density:
= U TS;(11)
which is again described fully in the same thermodynamic statespace fT;"
ij
g so that its total dierential reads:
d =
@
@"
kl
T
d"
kl
+
@
@T
"
dT;
d = dU dTS T dS =
ij
d"
ij
+T dS dTS T dS =
ij
d"
ij
S dT;(12)
@
@"
kl
T
=
ij
;
@
@T
"
= S:
3
Since the free energy is a scalar quantity (invariant in statespace):
@
2
@T@"
kl
=
@
2
@"
kl
@T
;
@
@T
kl
=
@
@"
kl
S;(13)
p
kl
=
@S
@"
kl
;
the heat developed in the body is also the same as the (negative) thermal pressure.Thus the thermal pressure
p
ij
may not depend neither the temperature nor the strains and the problem about the feasibility of the
measurements has been solved.We do not need to measure the coecients for the heat developed in the body,
measuring the thermal pressure is sucient.The Eq.(13)
2
is called the Maxwell relation.
5 Constitutive relations
The balance equations can be transformed onto the material frame:
J =
0
;
0
(v
j
)
@P
rj
@X
r
0
f
j
= 0;
0
u
+
@Q
r
@X
r
0
r = P
rj
@v
j
@X
r
;
(14)
where
F
i
r
=
@x
i
@X
r
;J = det F
i
r
;P
rj
= (F
i
r
)
1
J
ij
;Q
r
= (F
i
r
)
1
Jq
i
;
(15)
thus we need to dene the stress
ij
or P
rj
,the heat ux q
i
or Q
i
and the specic internal energy u regarding
primitive variables in actual or material frame.Let us use again the Fourier type heat equation and a
viscoelastic stress in material frame such that:
Q
i
=
ij
@T
@X
j
;P
ij
=
ij
+
ij
viscous
;
(16)
holds generally for anisotropic materials
5
.Especially the stress should be emphasized since a decomposition
into an elastic part
ij
which is fully recoverable and a viscous part
ij
viscous
which is irreversible is just a model.
Now start with the most generic case as explained in x2 how to measure
6
it:
ij
= C
ijkl
"
kl
+p
ij
T:(17)
If the stressstrain relation for a specic temperature is linear,i.e.,the plot of stressstrain curve on each tensile
testing (remember we need many of them) is a line
7
,then the coecients of the C
ijkl
is just some numbers (in
the unit of stress,i.e.,in MPa).If the plot reads a curve rather than a line this nonlinearity is expressed so
that the coecients are functions of strains or their invariants.These hyperelastic materials do not need more
experiments,they just need to be t on a function
8
depending on the invariants.From a thermodynamic point
of view the hyperelastic materials are not admissible.The dierential d
ij
can bring the constitutive relation
above if by the integration of it C
ijkl
is constant in"
kl
and p
ij
is constant T.The thermal pressure is a constant
due to the maxwell relations,however,the hyperelastic materials depend on the strains.We leave this aside
and employ here linear elastic materials.
For practical purposes the thermal pressure,p
ij
is not possible to measure at T = 0 K therefore the measure
ment takes place regarding a reference temperature T
ref.
.As an example let us choose an isotropic material,
like steel,aluminium,copper,etc.,and apply Hooke's law (linear model) for the elastic part and Duhamel
Neumann extension (again linear model) for the thermal part:
ij
= C
ijkl
"
kl
+p
ij
(T T
ref.
) = (
ij
kl
+
ik
jl
+
il
jk
)"
kl
+p
ij
(T T
ref.
) =
=
"
k
k
+p(T T
ref.
)
ij
+2"
ij
;"
kl
=
@u
(k
@X
l)
;
(18)
5
Try to distinguish between the rst PiolaKirchhoff stress P
ij
and the thermal stress p
ij
6
Of course we have to start from a (reference) state fT;"
ij
g where we know the stress,
ij
.The basic (but not practical)
assumption is that on zero strains and on zero absolute temperature stress vanishes.
7
For many engineering materials like AISI steels,copper,titanium,...this is the case up to the 0:35 of the melting temperature.
8
Cf.Comp.Real.XV
4
and use a Newtonian uid for the viscous part:
ij
viscous
= ("
ij
)
=
@v
(i
@X
j)
:(19)
The 1
st
Law does not include any statement for the latter viscous part.Therefore we are free to choose anything
at the moment.In short,however,we will see that the 2
nd
Law brings some restrictions to that,actually what
we have chosen here is not just any model at all!
We have explained the internal energy density:
dU =
0
du =
ij
d"
ij
+T
0
ds;
0
u
=
ij
@v
(i
@X
j)
+T
0
s
;
(20)
and by using the free energy made the entropy measurable:
dS =
0
ds = p
ij
d"
ij
+
C
"
T
dT;
0
s
= p
ij
@v
(i
@X
j)
+
0
c
"
T
T
:
(21)
For the isotropic material this reads:
p
ij
=
@
ij
@T
=
@
ij
@"
kl
@"
kl
@T
= C
ijkl
kl
= (
ij
kl
+
ik
jl
+
il
jk
)
kl
= (3 +2)
ij
;
(22)
where is the coecient of thermal expansion.Finally we get the sought after equation:
0
u
=
ij
@v
(i
@X
j)
T(3 +2)
ij
@v
(i
@X
j)
+
0
c
"
T
:
(23)
Now the balance of the internal energy,i.e.,the Eq.(14)
3
can be rewritten with the Eq.(23) such that:
ij
@v
(i
@X
j)
T(3 +2)
ij
@v
(i
@X
j)
+
0
c
"
T
+
@Q
r
@X
r
0
r = P
rj
@v
j
@X
r
;
(24)
by using P
rj
=
ij
+
ij
viscous
,since
ij
=
ji
and the Eq.(16)
1
,where
ij
=
ij
due to the isotropic material
response,the dierential equation to be solved consists of only primitive variables:
T(3 +2)
ij
@v
(i
@X
j)
+
0
c
"
T
@
@X
r
@T
@X
r
0
r =
@v
(r
@X
j)
@v
j
@X
r
;v
i
= u
i
:(25)
The 2
nd
Law of Thermodynamics asserts that any process should produce a positive entropy source.Hence
we need to bring the entropy in a form of a balance equation.We take the Eq.(14)
3
and employ the Eq.(20)
2
into:
ij
@v
(i
@X
j)
+T
0
s
+
@Q
r
@X
r
0
r = P
rj
@v
j
@X
r
;
0
s
+
1
T
@Q
r
@X
r
1
T
0
r =
rj
viscous
@v
j
@X
r
;
0
s
+
@
@X
r
Q
r
T
1
T
0
r =
rj
viscous
@v
j
@X
r
+Q
r
@
@X
r
1
T
;
0
s
+
@
@X
r
Q
r
T
=
1
T
0
r +
rj
viscous
@v
j
@X
r
Q
r
T
2
@T
@X
r
;
(26)
and get a balance for the entropy.The right side is the entropy source and according to the 2
nd
Law it is greater
than zero for irreversible and is equal to zero for reversible processes.If the body has everywhere the same
temperature and no viscous part in the stress tensor (elasticity!) then the process supposed to be reversible
if no supply term exists.Here we see that the Comp.Real.IX has a small misinterpretation by assuming that
the entropy does not change in elasticity,since the laser supply as well as the nonisothermal process should
produce entropy!
5
The most important restriction from the 2
nd
Law can be seen if we utilize the Eqs.(19) and (16),i.e.,the
constitutive equations of an isotropic body for viscous yielding and heat exchange to the entropy source:
=
1
T
0
r +
rj
viscous
@v
j
@X
r
Q
r
T
2
@T
@X
r
=
1
T
0
r +
@v
(j
@X
r)
@v
j
@X
r
+
T
2
@T
@X
r
@T
@X
r
:(27)
Since the material framework is expressed in Cartesian coordinates,the coand contravariant notation is not
needed.Moreover a double contraction of a symmetric tensor with a skewsymmetric tensor vanishes so that
the latter:
=
1
T
0
r +
@v
(j
@X
r)
@v
(j
@X
r)
+
T
2
@T
@X
r
@T
@X
r
;(28)
is always positive,i.e., 0 if the material coecients and are positive valued.Temperature is in Kelvin,
thus always positive,there exists also not minus supply term for energy so that a suction machine could bring
a body (volumetrically) into a lower temperature.The negative heat exchange performs over the boundary but
not volumetrically.
6 Application
We want to summarize and write the variational form for a generic case in material conguration expressed in
Cartesian coordinates.The balance Eqs.(14)
2
and (25) contracted with appropriate test functions should be
computed for a continuum B
0
over time:
Z
B
0
u
j
0
(v
j
)
@P
rj
@X
r
0
f
j
+T
T(3 +2)
ij
@v
(i
@X
j)
+
0
c
"
T
@
@X
r
@T
@X
r
0
r
@v
(r
@X
j)
@v
j
@X
r
!
dV
0
= 0:
(29)
Now by using integration by parts twice the latter reads:
Z
B
0
u
j
0
@v
j
@t
+
@u
j
@X
r
P
rj
u
j
0
f
j
TT(3 +2)
ij
@v
(i
@X
j)
+T
0
c
"
T
+
+
@T
@X
r
@T
@X
r
T
0
r T
@v
(r
@X
j)
@v
j
@X
r
!
dV
0
+
I
@B
0
u
j
P
rj
T
@T
@X
r
!
N
r
dA
0
= 0;
(30)
which is then in a discretized form for the displacements u
i
and temperatures T by omitting the coand con
travariant notation due to the Cartesian coordinates and using a compact notation for dierentiation:
X
ele.
Z
B
ele.
0
u
j
0
u
j
2u
0
j
+u
00
j
dt dt
+u
j;r
u
k;k
+(3 +2)(T T
ref.
)
jr
+2u
(r;j)
+
dt
(u
(r;j)
u
0
(r;j)
)
u
j
0
f
j
T(T T
ref.
)(3 +2)
dt
u
(i;i)
u
0
(i;i)
+T
0
c
"
T T
0
dt
+T
;r
T
;r
T
0
r
T
dt dt
u
(r;j)
u
0
(r;j)
u
j;r
u
0
j;r
!
dV
0
+
I
@B
0
u
j
^
t
j
T
^
h
!
dA
0
= 0;
(31)
with given boundary conditions
^
t
j
= P
rj
N
r
and
^
h = T
;r
N
r
.This can be solved for a tensile test as in Fig.1.
Suppose that there is a beam of 100 20 20 mm expressed in a Cartesian coordinate system,where
x 2 [0;100] on its axial and y;z 2 [10;+10] on its normal surface are positioned.By holding it on its left,
x = 0,clamped,i.e.,u
i
= 0 and loading on its right,x = 100,by a traction linear in time,
^
t = 800 t MPa,
the beam will be longer.The traction (force per area) is the controlled parameter,i.e.,the tensile loading
machine is steered by the force and the displacement at the right tip is measured,i.e.,observed parameter.It
is convenient to plot stress vs.strain,where the stress (on the right tip) is the (axial) traction and the (normal
axial) strain,i.e.,"
xx
is the displacement over initial length (100 mm).This is the engineering strain,it is easy
6
Figure 1:Tensile loading without viscosity.
to calculate and is known as the GreenLagrange strain measure,which we implemented.There is also the
true strain to be calculated as the displacement over the actual length,known as the EulerAlamnsi strain
measure (cf.Comp.Real.II for their derivation).
The time,t,is starting fromzero and it is of importance for viscous but not for the elastic case.If the material
is purely elastic, = 0,then the material responds to the strain and the same curve as in Fig.1 comes out for
any rate of loading.This is not the same in the viscous case, 6= 0,where the relation is curved (cf.Fig.2)
and the curvature depends on the rate of loading.Theoretically the temperature is coupled to the strains and
strain rates.When we apply a Neumann condition so that the initially from the room temperature (293:15
K) started beam is interchanging energy (heat ux) over all of its boundaries to obtain the outer temperature
(293:15 K),then this realistic conditions show us the temperature change in the beam.In Fig.1 and Fig.2,
respectively for the purely elastic and elastoviscous case,it is obvious that the deformation induces a tem
perature (the observed parameter) drop.However,it is so small that this change can easily be neglected for
materials like steel,aluminum,copper,titanium type of crystallographic engineering materials.This happens
due to the elastic dilation,"
ii
,i.e.,the volume contraction in small deformations,which reduces the entropy.
The mechanical part of the entropy change,p
ij
"
ij
,for the isotropic materials,p
ij
= p
ij
,results a reduction
in entropy,p
ij
"
ij
= p"
ii
< 0,since the thermal pressure,p = (3+2) > 0,is a positive quantity for many
materials.There are some fancy materials where the coecient of thermal expansion,,is negative so that s
temperature rise induces a volume contraction instead of an expansion,however,it is contraintuitive because
it is happening so rarely that only a small amount of material scientists are\enjoying"them in labs.
For practical purposes,a deformation which reads a volume enlargement (Poisson's ratio = 0:5 reads an
isochoric,volume preserving,and < 0:5 reads an enlargement),lowers the entropy which induces a temperature
drop,depending on the specic heat capacity c
"
.Though this is so small that it be neglected and this is good
news because the elasticity components,; for isotropic materials (81 in the generic case,C
ijkl
components),
need to be measured in constant temperature (cf.Eq.(3)).Hence we can use a tensile test to determine the
material parameters.
7
Figure 2:Tensile loading with viscosity.
7 Code
"""Co mput a t i o na l Re a l i t y X,t he r mo dy na mi c a l a n a l y z e o f a t e n s i l e l o a d i n g"""
aut hor
="B.
Emek
Abal i
( abal i @tube r l i n.de )"
dat e
="20120927"
c opyr i ght
="Copyri ght
(C)
2012,
B.
Emek
Abal i"
l i c e n s e
="GNU
GPL
Versi on
3
or
any
l a t e r
ver s i on"
#Co p y r i g h t ( C) 2012 B.Emek Ab a l i
#
#Thi s c o de i s f r e e s o f t wa r e:you can r e d i s t r i b u t e i t and/o r mo di f y
#i t unde r t he t e r ms o f t he GNU L e s s e r Ge n e r a l Pu b l i c L i c e n s e as p u b l i s h e d by
#t he Fr e e S o f t wa r e Foundat i on,e i t h e r v e r s i o n 3 o f t he Li c e n s e,o r
#( a t your o p t i o n ) any l a t e r v e r s i o n.
#
#Thi s c o de i s d i s t r i b u t e d i n t he hope t h a t i t w i l l be u s e f u l,
#but WITHOUT ANY WARRANTY;wi t h o u t e ve n t he i mp l i e d wa r r a nt y o f
#MERCHANTABILITY o r FITNESS FOR A PARTICULAR PURPOSE.Se e t he
#GNU L e s s e r Ge n e r a l Pu b l i c L i c e n s e f o r more d e t a i l s.
#
#For t he GNU L e s s e r Ge n e r a l Pu b l i c L i c e n s e s e e <ht t p://www.gnu.o r g/l i c e n s e s/>.
#
#
from dol f i n import
import numpy
import pyl ab
parameters ["f orm
compi l er"] ["cpp
opti mi ze"] = True
parameters ["f orm
compi l er"] ["opti mi ze"] = True
8
parameters ["f orm
compi l er"] ["opti mi ze"] = True
parameters ["f orm
compi l er"] ["r epr es ent at i on"] ="quadrature"
parameters ["f orm
compi l er"] ["quadrature
degree"] = 2
parameters ['l i near
al gebr a
backend']='uBLAS'#'PETSc','uBLAS','Epe t r a','STL'
#a b s t r a c t i n t e r f a c e f o r a Newton s o l v e r,
#i t e r a t e i s an o b j e c t and has i t s own a t t r i b u t e s s e l f.a t t r i b u t e
cl ass i t e r a t e ( Nonl i nearProbl em):
def
i n i t
( s e l f,a,L,bc,exter
B,i nt er
B ):
Nonl i nearProbl em.
i n i t
( s e l f )
s e l f.L = L
s e l f.a = a
s e l f.bc = bc
s e l f.ext er = exter
B
s e l f.i nt e r = i nt er
B
def F( s e l f,b,x ):
assembl e ( s e l f.L,t ens or=b,e xt e r i or
f ac e t
domai ns=s e l f.ext er )
for condi t i on in s e l f.bc:condi t i on.appl y (b,x)
def J( s e l f,A,x ):
assembl e ( s e l f.a,t ens or=A,e xt e r i or
f ac e t
domai ns=s e l f.ext er )
for condi t i on in s e l f.bc:condi t i on.appl y (A)
cl ass i ni t
va l ue s ( Expressi on ):
def eval ( s e l f,vec,x ):
vec [ 0 ] = 293.15#K
vec [ 1 ] = 0.0#m
vec [ 2 ] = 0.0
vec [ 3 ] = 0.0
def val ue
shape ( s e l f ):
return ( 4,)
#u n i t s:mm,1000 kg=ton,s,MPa,mJ,K
del t a = I dent i t y ( 3)
f = Constant ( ( 0.0,0.0,9810.))
r = 0.
Tref = 293.15
#Ma t e r i a l da t a o f P265GH,
#a.k.a.St 4 5.8 [ VDI Wae r me at l as ],a t 2 9 3.1 5 K
#De n s i t y [ kg/mm^ 3 ]
rho0 = 7850.0E9
#Ther mal c o n d u c t i v i t y kappa [ mJ/( s mm K) ]
kappa = 57.0
#S p e c i f i c he a t c a p a c i t y c [ mJ/( t on K) ]
capaci ty = 430.0E6
#t h e r ma l e x p a n s i o n [ 1/K ] at 3 7 3.1 5 K
al pha = 12.2E6
#Young's modul us [ MPa ]
E = 211.E+3
#Po i s s o n's r a t i o [ ]
nu = 0.28
#t h e r ma l f l u x o v e r boundar y [ mJ/( s m^ 2 ) ]
h = 10.E3
#v i s c o s i t y [ MPa/s ]
eta = 3.E+6#v i s c o e l a s t i c c a s e
#e t a = 0.#e l a s t i c c a s e
tMax = 1.0
dt = 0.1
t = 0.0
9
#t he g e o me t r y [ mm]
xMin = 0.0
xMax = 100.0
xElements = 10
yMin = 10.
yMax = +10.
yElements = 10
zMin = +10.
zMax = 10.
zEl ements = 10
mesh = Box(xMin,yMin,zMin,xMax,yMax,zMax,xElements,yElements,zEl ements )
n = FacetNormal ( mesh)
l ength = abs (xMaxxMin)
D = mesh.topol ogy ( ).dim( )
u
Space = VectorFuncti onSpace ( mesh,'CG',1)
T
Space = Functi onSpace ( mesh,'CG',1)
Space = MixedFunctionSpace ( [ T
Space,u
Space ] )
l e f t = compi l e
subdomai ns ('x [ 0 ]
==
l')
l e f t.l = xMin
r i ght = compi l e
subdomai ns ('x [ 0 ]
==
l')
r i ght.l = xMax
back = compi l e
subdomai ns ('x [ 1 ]
==
l')
back.l = yMin
f r ont = compi l e
subdomai ns ('x [ 1 ]
==
l')
f r ont.l = yMax
bottom = compi l e
subdomai ns ('x [ 2 ]
==
l')
bottom.l = zMin
top = compi l e
subdomai ns ('x [ 2 ]
==
l')
top.l = zMax
al l
boundar i e s = compi l e
subdomai ns ('on
boundary')
i nt e r i o r = MeshFunction("ui nt",mesh,D1)
i nt e r i o r.s e t
a l l ( 0)
e xt e r i or = MeshFunction("ui nt",mesh,D1)
e xt e r i or.s e t
a l l ( 0)
l e f t.mark( ext er i or,1)
r i ght.mark( ext er i or,2)
bottom.mark( ext er i or,3)
top.mark( ext er i or,4)
f r ont.mark( ext er i or,5)
back.mark( ext er i or,6)
tHat = Expressi on ( ('At','0.','0.'),A=8.E2,t =0.)
bc = [ Di ri chl etBC( Space.sub ( 1),Constant ( ( 0.0,0.0,0.0) ),ext er i or,1),n
Di ri chl etBC( Space.sub ( 1 ).sub ( 1),Constant ( 0.0),ext er i or,2),n
Di ri chl etBC( Space.sub ( 1 ).sub ( 2),Constant ( 0.0),ext er i or,2) ]
dunkn = Tri al Functi on ( Space )
t e s t = TestFuncti on ( Space )
delT,del u = s p l i t ( t e s t )
unkn = Functi on ( Space )
unkn0 = Functi on ( Space )
unkn00 = Functi on ( Space )
T,u = s p l i t ( unkn)
T0,u0 = s p l i t ( unkn0)
T00,u00 = s p l i t ( unkn0)
10
unkn
i ni t = i ni t
va l ue s ( )
unkn00.i nt e r pol at e ( unkn
i ni t )
unkn0.as s i gn ( unkn00)
unkn.as s i gn ( unkn0)
i,j,k,l = i ndi c e s ( 4)
lambada = E nu/( 1.+ nu)/( 1. 2. nu)
mu = 0.5 E/( 1.+ nu)
Form = ( del u [ j ] rho0 (u2.u0+u00 ) [ j ]/dt/dt n
+ del u [ j ].dx( l ) ( ( lambadau[ k ].dx( k)+(3. lambada+2.mu) al pha (TTref ) ) del t a [ j,l ] n
+ 2.musym( grad (u ) ) [ l,j ] + eta/dt (sym( grad (u))sym( grad ( u0 ) ) ) [ l,j ] ) n
del u [ j ] rho0 f [ j ] n
delT(TTref ) ( 3. lambada+2.mu) al pha/dt (sym( grad (u))sym( grad ( u0 ) ) ) [ i,i ] n
+ delTrho0 capaci ty (TT0)/dt + delT.dx( l )kappaT.dx( l ) delTrho0r n
delTetaeta/dt/dt (sym( grad (u))sym( grad (u ) ) ) [ l,j ] ( grad (u)grad (u ) ) [ j,l ] )dx n
del u [ j ] tHat [ j ] ds ( 2) delThds ( 1) delThds ( 2) delThds ( 3) n
delThds ( 4) delThds ( 5) delThds ( 6)
Gain = de r i vat i ve (Form,unkn,dunkn)
pwd='/bea/c al c ul/CompReal/'
f i l e
u = Fi l e (pwd+'CR10
di spl.pvd')
f i l e
T = Fi l e (pwd+'CR10
temp.pvd')
s t r ai n,s t r e s s,temp =[ ],[ ],[ ]
f i g = pyl ab.f i gur e (1,f i g s i z e =(14,10))
ax1 = f i g.add
subpl ot (111)
ax1.gr i d ( True,axi s='x')
ax1.s e t
xl a be l ('engi neer i ng
s t r ai n
$u
1/l
0$')
ax1.s e t
yl a be l ('engi neer i ng
s t r e s s
$F/A$',c ol or='r')
ax1.ti ck
params ( axi s='y',c ol or s='r')
ax1.gr i d ( True,axi s='y',c ol or='r')
ax2 = ax1.twinx ( )
ax2.s e t
yl a be l ('temperature
$T$',c ol or='b')
ax2.ti ck
params ( axi s='y',c ol or s='b')
ax2.gr i d ( True,axi s='y',c ol or='b')
#ax2.s e t
y l i m ( ( 0.,0.1 ) )
while t <= tMax:
print'time:
',t
tHat.t = t
problem = i t e r a t e ( Gain,Form,bc,ext er i or,i nt e r i o r )
s ol ve r = NewtonSolver ('l u')
s ol ve r.parameters ['c onve r ge nc e
c r i t e r i on'] ='i ncremental'
s ol ve r.parameters ['r e l a t i ve
t o l e r a nc e'] = 1.0 e8
s ol ve r.parameters ['maxi mum
i terati ons'] = 25
s ol ve r.s ol ve ( problem,unkn.vector ( ) )
f i l e
T << ( unkn.s pl i t ( ) [ 0 ],t )
f i l e
u << ( unkn.s pl i t ( ) [ 1 ],t )
s t r ai n.append( unkn.s p l i t ( ) [ 1 ] ( xMax,0.,0.) [ 0 ]/l ength )
s t r e s s.append( tHat (xMax,0.,0.) [ 0 ] )
temp.append( unkn.s pl i t ( ) [ 0 ] ( xMax/2.,0.,0.) )
ax1.pl ot ( s t r ai n,s t r e s s,'o',c ol or='r')
ax2.pl ot ( s t r ai n,temp,'d',c ol or='b')
#f i g.s a v e f i g ('C R 1 0
t e n s i l e t e s t
e l a s t i c.png')
f i g.s ave f i g ('CR10
t e ns i l e t e s t
vi s c ous.png')
11
unkn00.as s i gn ( unkn0)
unkn0.as s i gn ( unkn)
t += dt
#e o f
12
Enter the password to open this PDF file:
File name:

File size:

Title:

Author:

Subject:

Keywords:

Creation Date:

Modification Date:

Creator:

PDF Producer:

PDF Version:

Page Count:

Preparing document for printing…
0%
Comments 0
Log in to post a comment