Thermodynamics - TU Berlin

acridboneMechanics

Oct 27, 2013 (3 years and 10 months ago)

156 views

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 dierent 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 (non-zero source)
which is always there if the material undergoes a motion.One possible aim is to calculate the mass density,
the velocities and the specic 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 specic (per mass) energy,e,is additive so that the specic internal energy,u,
and specic kinetic energy,e
kin.
,are separable and uncoupled.The author believes that the kinetic energy
is perturbing the material in the non-recoverable way and the change of the internal energy
2
is recoverable.
There may be dierent formulations,however,none of them is wrong because we cannot measure dierent
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 space-time.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 state-space and plot
4
the specic
internal energy u = u(T;"
ij
) in that ten-dimensional 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 three-axis 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 dierential and the 1
st
Law of thermodynamics states that the internal energy is a
total dierential,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 dened.Since the assumption that
the set fT;"
ij
g describes the state fully,the stress 
ij
and the entropy density S should be dened also in that
state-space,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 stiness 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 coecients in three-dimensional 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 stress-strain 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
specic 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 specic
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 three-dimensional case the 81 components of stiness
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 dicult 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) dierential 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 dierential perfect or total.By having a total dierential 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 dierential form is essential in analysis.Suppose that there are additional conditions
to be fullled,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 AB  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 AB  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 dened 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 state-space fT;"
ij
g so that its total dierential 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 state-space):
@
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 coecients for the heat developed in the body,
measuring the thermal pressure is sucient.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 dene the stress 
ij
or P
rj
,the heat ux q
i
or Q
i
and the specic 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 stress-strain relation for a specic temperature is linear,i.e.,the plot of stress-strain curve on each tensile
testing (remember we need many of them) is a line
7
,then the coecients 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 non-linearity is expressed so
that the coecients 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 dierential 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 Piola-Kirchhoff 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 coecient 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 dierential 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 non-isothermal 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 co-and contravariant notation is not
needed.Moreover a double contraction of a symmetric tensor with a skew-symmetric 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 coecients  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 conguration 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 co-and con-
travariant notation due to the Cartesian coordinates and using a compact notation for dierentiation:
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
+2u
(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 Green-Lagrange strain measure,which we implemented.There is also the
true strain to be calculated as the displacement over the actual length,known as the Euler-Alamnsi 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 elasto-viscous 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 coecient of thermal expansion,,is negative so that s
temperature rise induces a volume contraction instead of an expansion,however,it is contra-intuitive 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 specic 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 @tube r l i n.de )"
dat e
="20120927"
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.0E9
#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.2E6
#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.E3
#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 (xMaxxMin)
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,D1)
i nt e r i o r.s e t
a l l ( 0)
e xt e r i or = MeshFunction("ui nt",mesh,D1)
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 ( ('At','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 (u2.u0+u00 ) [ j ]/dt/dt n
+ del u [ j ].dx( l ) ( ( lambadau[ k ].dx( k)+(3. lambada+2.mu) al pha (TTref ) ) del t a [ j,l ] n
+ 2.musym( grad (u ) ) [ l,j ] + eta/dt (sym( grad (u))sym( grad ( u0 ) ) ) [ l,j ] ) n
 del u [ j ]  rho0 f [ j ] n
 delT(TTref ) ( 3. lambada+2.mu) al pha/dt (sym( grad (u))sym( grad ( u0 ) ) ) [ i,i ] n
+ delTrho0 capaci ty (TT0)/dt + delT.dx( l )kappaT.dx( l )  delTrho0r n
 delTetaeta/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)  delThds ( 1)  delThds ( 2)  delThds ( 3) n
 delThds ( 4)  delThds ( 5)  delThds ( 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 e8
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