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 (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 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 non-recoverable 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 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 specic

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 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

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 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 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

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 three-dimensional 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 state-space 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 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 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 stress-strain relation for a specic 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 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 non-linearity 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 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 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 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 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 co-and 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 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 coecient 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 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

## Comments 0

Log in to post a comment