Previous Up Next

2  Equations

2.1  Convection and diffusion of heat

2.1.1  Convection-diffusion equation



ρ C ( T + βi
T
xi
) = k (
2 T
x12
+
2 T
x22
+
2 T
x32
) − a T + f


The primary unknown is the condif_temperature T. Further notation: ρ group_condif_density; C group_condif_capacity; x space coordinate; βi group_condif_flow in i-direction; k group_condif_conductivity; a group_condif_absorption; f force_element_volume. Typical applications are heat conduction and heat conduction in a flow.

2.1.2  Convection to environment



qc = αc ( TTc )


Here qc is the condif_convection heat flux, αc is the convection coefficient and Tc is the environmental temperature for convection.

2.1.3  Radiation to environment



qr = αr ( T4Tr4 )


Here qr is the condif_radiation heat flux, αr is the radiation coefficient and Tr is the environmental temperature for radiation.

2.2  Material deformation and flow



ρ vi =
∂ σij
xj
+ (1−β T) ρ gid
vi
xi
no   sum   on   i + fi


Notations: ρ materi_density; vi materi_velocity in i-direction; σij materi_stress matrix; x space coordinate; β group_materi_expansion_volume; T (optional) condif_temperature; gi force_gravity; d group_materi_damping coefficient; fi force_element_volume. The equation is given for space coordinates following the material velocities vi.

TOCHNOG allows you to build your favorite material, by adding separate contributions to the stresses σij. In this way you can build solids or fluids or things in between. The separate contributions will be listed below. First two typical examples are given.

Nearly incompressible Navier Stokes:

...
materi_velocity
materi_stress
end_initia
...
group_type 0 -materi
group_materi_elasti_compressibility 0 1.0
group_materi_viscosity 0 1.2
Linear solid:

...
materi_velocity
materi_strain_total
materi_stress
end_initia
...
group_type 0 -materi
group_materi_elasti_young 0 1.e10
group_materi_elasti_poisson 0 0.2
group_materi_memory 0 -updated_without_rotation

2.2.1  Memory

The -updated Lagrange formulation

Deformations (i.e. the incremental deformation matrix F) refers to the previous time point. TOCHNOG decomposes the incremental deformation tensor with a polar decomposition into F = R U with F the incremental deformation matrix, R the incremental rotation matrix and U the incremental stretch matrix. The incremental stretch matrix U is used to determine the incremental strain matrix 0.5(U+UT)−I with I the identity tensor. The velocities between the two time points are the unknowns to be solved. Stresses are calculated from adding incremental stresses to the old stresses. Incremental stresses are caused by the incremental strain matrix and a rigid body rotation by the incremental rotation matrix of the old stresses.

The -updated_without_rotation Lagrange formulation

Deformations (i.e. the incremental deformation matrix F) refers to the previous time point. Any rigid body rotation between the two time points are neglected, so TOCHNOG decomposes the incremental deformation tensor with a polar decomposition into F = U with F the incremental deformation matrix, and U the incremental stretch matrix. The velocities between the two time points are the unknowns to be solved. Stresses are calculated from adding incremental stresses to the old stresses. Incremental stresses are caused by the incremental strain matrix.

The -total Lagrange formulation

Deformations (i.e. the total deformation matrix F) refers to the time 0. TOCHNOG decomposes the total deformation tensor with a polar decomposition into F = R U with F the total deformation matrix, R the total rotation matrix and U the total stretch matrix. The total stretch matrix U is used to determine the total strain matrix 0.5(U+UT)−I with I the identity tensor. The displacements between the current time and time 0 are the unknowns to be solved. Stresses are calculated from the total displacements.

The -total_piola Lagrange formulation

Deformations (i.e. the total deformation matrix F) refers to the time 0. TOCHNOG decomposes the total deformation tensor with a polar decomposition into F = R U with F the total deformation matrix, R the total rotation matrix and U the total stretch matrix. The total stretch matrix U is used to determine the total Green-Lagrange strain matrix 0.5 ( FT FI ) = 0.5 ( UT UI ) Stresses are calculated from the total displacements.

The -total_linear Lagrange formulation

Deformations (i.e. the total deformation matrix F) refers to the time 0. TOCHNOG neglects any rigid body rotations and uses linear engineering strains 0.5 ( F + FT ) − I. The displacements between the time current time and time 0 are the unknowns to be solved. Stresses are calculated from the total displacements.

A remark on the total Lagrange models. Normally stresses are calculated from the total displacements (and thus the total strains). The old stresses are not used. This means that any initial stresses are neglected. This type of stress calculation for the total Lagrange models is used whenever materi_strain_total is not initialised. However, in case materi_strain_total is initialised, the difference between the total displacements (and thus strains) between two time points is used to determine incremental stresses, which are added to the stresses at the previous time point. And thus, in case materi_strain_total is initialised, the old stresses are used and any initial stresses are not neglected.

See also group_materi_memory.

2.2.2  Elasticity

The elastic stress rate is

Cijkl єklelas


where Cijkl is the elastic modulus tensor (which is a doubly symmetric tensor: Cijkl=Cjikl, Cijkl=Cijlk and Cijkl=Cjilk), and єklelas is the elastic strain rate. See the plasticity section for a definition of the elastic strain rate.

For an isotropic material

C0000 = C1111 = C2222 =
E(1−ν)
(1+ν)(1−2ν)
C0011 = C0022 = C1122 =
Eν
(1+ν)(1−2ν
C0101 = C0202 = C1212 =
E
1+ν


with E group_materi_elasti_young modulus and ν group_materi_elasti_poisson ratio (the remaining non-zero moduli follow from the double symmetry conditions).

For a transverse isotropic material the material has one unique direction (think of an material with fibers in one direction). Here we take 'a' as the unique direction; 'b' and 'c' are the transverse directions. The material is fully defined by Caaaa, Cbbbb, Caabb, Cabab and Cbcbc and the unique direction in space (see group_materi_elasti_transverse_isotropy). The other non-zero moduli follow from Ccccc=Cbbbb, Cacac=Cabab, Cbbcc=Cbbbb−2Cbcbc and from the double symmetry conditions.

The nonlinear elasticity polynomials is a strain dependent model. In this model, the 'young's stiffness' modulus is made dependend of the size of the strains via a series of polynomials

E = E0 + E1 є1 + E2 є2 + …     (1)


where

є =
(
єij єij )     (2)


with єij the components of the strain matrix. The parameters E0 etc. need to be specified in the group_materi_elasti_young_polynomial record.

The power law nonlinear elasticity is a stress dependent model which typically is used to model the elastic behavior of granular materials. It can be combined with plastic models, by example with the di Prisco plasticity model for soils, and with a poisson ratio.

In this model, the 'young's stiffness' modulus is made a function of the average stress state:

E = E0 (p/p0)α    (3)


where p is the pressure. Furthermore, E0 is the reference stiffness at reference pressure p0, and α is a soil dependent power coefficient. The parameters E0, p0, and α need to be specified in the group_materi_elasti_young_power record.

The Lade nonlinear elasticity is a stress dependent model which typically is used to model the elastic behavior of granular materials. It can be combined with plastic models, by example with the di Prisco plasticity model for soils.

The stress rates are linked to the strain rates by the equation:

єij =
W2
∂ σij ∂ σhk
  σhk     (4)


where the function W is

W =
X1−λ
2 B (1−λ)


where

X = p2 + R* abs( sij sij )


with pressure p=(σ112233)/3 and deviatoric stresses sij = σijp δij.

The model contains three user specified constants B, R, λ which need to be specified in the group_materi_elasti_lade record. B and λ are defined by means of an isotropic unloading test, and R by means of an unloading-standard-triaxial-compression test. For example for a loose sand B=1028, R=0.25, λ=0.28. See [7] for the details.

The model cannot be used in combination with a poisson ratio.

2.2.3  Elasto-Plasticity

Plastic strain


In plastic analysis, the materi_strain_elasti rate follows by subtracting from the materi_strain_total rate the materi_strain_plasti rate

єijelas = єij − єijplas


where the materi_strain_total rate is

єij = 0.5 (
vi
xj
+
vj
xi
)


The materi_strain_plasti rate follows from the condition that the stress cannot exceed the yield surface. This condition is specified by a yield function fyieldij)=0. The direction of the plastic strain rate is specified by the stress gradient of a flow function ∂ fflow/∂ σij. If the yield function and flow function are chosen to be the same, the plasticity is called associative, else it is non-associative.

Von-Mises is typically used for metal plasticity. Mohr-Coulomb and Drucker-Prager are typically used for soils and other frictional materials. The plasticity models can freely be combined; the combination of the plasticity surfaces defines the total plasticity surface.

First some stress quantities which are used in most of the plasticity models are listed.

Equivalent Von-Mises stress:
σ =
sijsij
2
Mean stress:
σm =
σ11 + σ22 + σ33
3
Deviatoric stress:
sij = σij − σm δij


Modified CamClay plasticity model


Here we give the equations of the Modified Cam Clay model (Roscoe and Burland, 1968, summarized e.g. by Wood, 1990, see [12]). All stresses are effective (geotechnical) stresses, i.e.compression is positive! Definitions of variables:

p = (σ112233)/3


q = {
1
2
[ (σ11−σ22)2 + (σ22−σ33)2 + (σ33−σ11)2 ] + 3 ( σ122 + σ232 + σ312 ) }1/2


The CamClay yield rule, which is also the flow rule, reads:

f = g = q2M2 [ p (p0p) ] = 0


M is a soil constant and p0 is a history (hidden) variable which corresponds to the preconsolidation mean pressure. The hardening function, evolution, of p0 reads:

d p0 =
p0 (1+e) dεvp
λ−κ


in which

dεvp = dε11p+dε22p+dε33p


and λ and κ are user specified soil constants. Further e is the void ratio with the evolution equation:

de = −dεv (1+e)


in which

dεv = dε11+dε22+dε33


The poisson ratio ν reads:

ν =
3K − 2G
2G+6K


in which the elastic bulk modulus K is given by:

K = (1+e) p / κ


and the Young's modulus E:

E = 2.*G*(1+ν)


in which G is a user specified soil constant, By using this ν and E the classical isotropic stress-strain law is used to calculate the stresses.

The soil constants M, κ, λ need to be specified in group_materi_plasti_camclay. The soil constant G, need to be specified in group_materi_elasti_camclay_g. For an alternative see group_materi_elasti_camclay_poisson. The history variables e, p0 need to be initialized by materi_history_variables 2 record (and given initial values in node_dof records).

Remark 1: An additional parameter N can be often found in textbooks on the Cam Clay model. We don't include it since it is linked to other model parameters via:
1+e = N − λ lnp0 + κ ln(p0/p)


Remark 2: If you apply a geometrical linear analysis, see section 8.4, then also the calculation of de void ratio development is linearized, and so will contain some error as compared to the exact void ratio change. Hence for very large deformations, say above 10 percent or so, don't use such geometrical linear analysis.

Cap plasticity model


This model accounts for permanent plastic deformations under high pressures for granular materials. It is intended to be used in combination with shear plasticity models like Drucker-Prager, Mohr-Coulomb, etc.

First a deviatoric stress measure t and hydrostatic stress measure p are defined
t =
3
σ
p = − σm


See above for σ and σm. The yield rule for the group_materi_plasti_cap model reads:

f =
(ppa)2 +





R t
(1+α−
α
cosϕ






2






 
R ( c + pa tanϕ )


Here c is the cohesion and ϕ is the friction angle which should be taken equal to the values in the shear flow rule which you use. The parameter pa follows from

pa =
pbRc
1 + R   tanϕ


where the hydrostatic compression yield stress pb is to be defined with an table of volumetric plastic strains єvp versus pb with єvp = є11p + є22p + є33p. As always, positive strain denote extension whereas negative strains denote compression.

Associative flow is used, so the flow rule is taken equal to the yield rule.

Summarizing the group_materi_plasti_cap model needs as input the cohesion c, the friction angle ϕ, the parameter α (typically 1.  10−2 up to 5.   10−2), and a table єvp versus pb.

Compression limiting plasticity model


This group_materi_plasti_compression model uses a special definition for the equivalent stress

σ =
σ12 + σ22 + σ32


where σ1, σ2 and σ3 are the first, second and third principal stress respectively. Each of these is only incorporated if it is a compression stress. The model now reads

σ − σy = 0


This plasticity surface limits the allowed compression stresses.

Tension limiting plasticity model


This group_materi_plasti_tension model uses a special definition for the equivalent stress

σ =
σ12 + σ22 + σ32


where σ1, σ2 and σ3 are the first, second and third principal stress respectively. Each of these is only incorporated if it is a tension stress. The model now reads

σ − σy = 0


This plasticity surface limits the allowed tension stresses.

A simple model for concrete can be obtained as follows. Use group_materi_plasti_tension to limit the tension strength ft. Use group_materi_plasti_compression to limit the compressive strength fc. The tension strength could be softened to zero over an effective plastic strain κ of, say, 1 percent. The compressive strength could be softened to zero over an effective plastic strain κ of, say, 10 percent.

Another possibility for concrete is to combine group_materi_plasti_tension to limit the tension strength ft, and use the group_materi_plasti_vonmis model to limit the compressive strength fc.

di Prisco plasticity model


The di Prisco model is an non-associative plastic model for soils, which can be typically combined with the 'Lade elastic model'. This di Prisco model is a rather advanced soil model, which is explained in more detail in [4] and [5]. The yield rule reads:

f = 3 βf (γ − 3) ln


r
rc



− γ J3 η* +
9
4
( γ − 1 ) J2 η*


and the flow rule yields:

g = 9 ( γ − 3 ) ln


r
rg



− γ J3 η* +
9
4
( γ − 1 ) J2 η*


This is an anisotropic model in which the first and second invariant of the stress rate η* are defined relative to the rotation axes χ.

r = σij χij
J* = ηij* ηjk* ηki*
J* = ηij* ηij*
ηhk* =
3
shk*
r


where s* follows from

shk* = σhk*r χhk


Further rg=1.

The history variables are χij ( rotation axes, 9 values), β (yield surface form factor), and rc (preconsolidation mean pressure). The evolution laws for these history variables can be found in the papers listed above. The history variables χij (9 values), β, rc need to be initialized by the materi_history_variables 11 record (and should be given initial values in node_dof records). In a normally consolidated sand with isotropic initial conditions χij = δij / √3 , β=0.0001 and rc equals √3 times the means pressure.

The total model, yield rule and flow rule and evolution laws for history variables, contains a set of soil specific constants. In group_materi_plasti_diprisco you need to specify these constants. These constants are explained in more detail in the papers mentioned above, but here we give a short explanation. The constants θc, θe, ξc and ξe are linked to the dilitancy and the stress state during failure (standard triaxial compression and extension test in drained conditions). The constants γ, cp, βf and βf0 are defined by means of the experimental curves ( qaxial, єvolaxial) obtained by performing a standard compression test in drained conditions. Moreover, βf, βf0 and tp can also be determined by means of the effective-stress path obtained by performing a standard triaxial compression test in undrained conditions.

A cohesion C can be also be introduced if required.

Finally bp can determined from an isotropic compression test. For a loose sand θc=0.253, θe=0.0398, ξc=−0.2585, ξe=−0.0394, γ=3.7, cp=18., βf=0.5, βf0=1.1, tp=10., and bp=0.0049.

Drucker-Prager plasticity model


The group_materi_plasti_druckprag model reads

3 α σm + σ − K = 0
α =
2 sin( ϕ )
3
( 3 − sin(ϕ) )
K =
6 c cos( ϕ )
3
( 3 − sin(ϕ) )


Here c is the cohesion, which needs to be specified both for the yield function and the flow rule; by choosing different values non-associative plasticity is obtained.

Gurson plasticity model


The group_materi_plasti_gurson model reads
3 σ2
σy2
+ 2 q1 f* cosh( q2
3 σm
2 σy
) − (1 + ( q3 f* ) 2 ) = 0


Here f* is the volume fraction of voids. The rate equation
f* = ( 1 − f*) f* єkkplas
defines the evolution of f* if the start value for f* is specified. Furthermore, q1, q2 and q3 are model parameters.

HLC plasticity model


The group_materi_plasti_hlc Hau-Liu-Chang model is nearly similar to the Gurson's equation, and it is reads

3 σ2
σy2
+ f* ( +
1
m20
) m1 exp(
3 σm
2 σy
) − 1 = 0


Here f* is the volume fraction of voids.

The rate equation
f* = ( 1 − f*) f* єkkplas
defines the evolution of f* if the start value for f* is specified. The variable m20 is a function of the porosity:

m20 = ( m21m22 f*) exp( m21
σm
σ
)


Furthermore, σy, m21, m22, m23 and m1 are model parameters. The parameters can be calculated numerically by simple numerical simulations.

See http://www.tam.nwu.edu/wkl/paper/suhao-paper2.html for further info about the model.

Von-Mises plasticity model


The group_materi_plasti_vonmis model reads

3
  σ − σy = 0


where without hardening the yield value is fixed σy = σy0 .

If however the group_materi_plasti_vonmis_nadai hardening law for Von-Mises plasticity is specified then

σy = σy0 + C ( κ_0 + κ ) n


where C, κ_0 and n are parameters for the hardening law, and κ is the isotropic hardening parameter (see later). The parameter σy0 is specified by group_materi_plasti_vonmis.

Modified Matsuoka-Nakai model plasticity model


The Matsuoka-Nakai model [8] is a perfectly plastic model thus the fixed yield surface represents the failure surface as well. The model is based on experimental results with soils and can be formulated in terms of three stress invariants
f = I3 +
cos2 ϕ
9−sin2 ϕ
  I1   I2 = 0


where

I1 = tr(σij) = σ112233 = σ1 + σ2 + σ3 = 3 σm
I2 =
1
2

tr ( σikσkj ) − I12
= −σ1 σ2 − σ2 σ3 − σ3 σ1
I3 = det(σij) = σ1 σ2 σ3


σ1, σ2 and σ3 are the principal stresses (all stresses are effective; compressive stresses are negative). The parameter ϕ is equal to the angle of internal friction in axisymmetric (triaxial) compression [11].

For axisymmetric stress states the Matsuoka-Nakai model corresponds to the Mohr-Coulomb model. Nevertheless, the Matsuoka-Nakai model is described by a smooth surface in the stress space and thus it is more suitable from the computational aspect.

When the cohesion c is considered in the model, the yield condition is formulated for a modified stress [9]
σij = σij−σ0δij


with

σ0 = c cotϕ  .


Notice that straightforward usage of the function f above also as flow rule g would lead to a loss of the information on ϕ in the derivative of g, since the ϕ only appears in a constant.

For that reason we apply as flow rule g the Drucker-Prager function (see elsewehere in this theoretical part). A separate ϕflow can be specified which enters this Drucker-Prager g. Since we use the Drucker-Prager function as flow rule g, the present plasticity model has been named Modified Matsuoka-Nakai.

Mohr-Coulomb plasticity model


The group_materi_plasti_mohrcoul model reads

0.5 ( σ1 − σ3 ) + 0.5 ( σ1 + σ3 ) sin( ϕ ) − c   cos( ϕ ) = 0


Here c is the cohesion, σ1 is the maximal principal stress and σ3 is the minimal principal stress. The angle ϕ needs to be specified for both the yield condition and the flow rule; by choosing different values, non-associative plasticity is obtained.

For a numerically more stable solution, consider using Matsuoka-Nakai plasticity in stead of Mohr-Coulomb.

Mohr-Coulomb softening plasticity model


The group_materi_plasti_mohrcoul_softening model is the same as the standard Mohr-Coulomb model. Now, however, the parameters c and ϕ (both for the yield rule and for the flow rule) are softened on the the effective plastic strain κ.

By example, for the cohesion a linear variation is taken between the initial value c0 at κ=0, up to c1 at a specified critical value of κ, and constant c1 for larger values of κ. The same is done for ϕ for the yield rule and for the flow rule.

Isotropic Hardening


The size of the plastic strains rate is measured by the materi_plasti_kappa parameter
κ =
0.5 єijplas єijplas


This parameter can be used for isotropic hardening. Use the dependency_diagram for this.

Kinematic Hardening


The materi_plasti_rho matrix ρij, governs the kinematic hardening in the plasticity models. It is used in the yield rule and flow rule to get a new origin by using the argument σij − ρij:

fyield = fyieldij − ρij)
fflow = fflowij − ρij)
where the rate of the matrix ρij is taken to be

ρij = a    єijplas
where a is a user specified factor (see group_materi_plasti_kinematic_hardening).

Plastic heat generation


The plastic energy loss can be partially turned into heat rate per unit volume q:

q = η σij єijplas


where η is a user specified parameter (between 0 and 1) specifying which part of the plastic energy loss is turned into heat (see group_materi_plasti_heatgeneration).

2.2.4  Hypo-Plasticity

In hypoplasticity a direct relation is used between strain rates and stress rates. Specifically:

σij = Lijkl єij + Nij
єkl єkl


Here the part with Lijkl gives a linear relation between strain rates and stress rates and the part with Nij gives a nonlinear relation. The constitutive tensors Lijkl and Nij are functions of the effective stress tensor σij and void ratio e. The effective stress tensor σij follows from the total stress tensor σij minus any pore pressures (see groundflow). Rigid body rotations (objectivity) are treated elsewhere (see the section on memory).

Basic law, Wolffersdorff

The law proposed by Wolffersdorff [11] is used.

Lijkl = fs
1
σmnσmn

F2 Iijkl + a2   σijσkl
 
Nij = fs   fd
F a
σklσkl

 σij  +   σij*  
 
 
with    σij = σij / (σmnδmn)    ,    σij* = σij
1
3
 δij    ,    Iijkl = δik δjl  ,
 
 
a =
3
(3−sinφc)
2
2
sinφc
 
 
F  = 
1
8
tan2ψ+
2−tan2ψ
2+
2
tanψcos3θ
1
2
2
tanψ  ,
 
 
   tanψ=
3
σij*σij*
  ,    cos3θ=−
6
 
σij*σjk*σki*
[ σmn*σmn* ] 3/2
   .
For σij*=0 is F=1.

The scalar factors fs and fd take into account the influence of mean pressure and density:
fs =
hs
n
 


ei
e



β



 
1+ei
ei
 


σijδij
hs



1−n



 
 





3+a2a
3



ei0ed0
ec0ed0
 


α



 






−1






 
   ,
fd =
 


eed
eced
 


α



 
   .


Three characteristic void ratios – ei (during isotropic compression at the minimum density), ec (critical void ratio) and ed (maximum density) – decrease with mean stress:
ei
ei0
=
ec
ec0
=
ed
ed0
= exp








σijδij
hs



n



 








The range of admissible void ratios is limited by ei and ed. The model parameters can be found in Tab. 1. They correspond to Hochstetten sand from the vicinity of Karlsruhe, Germany [11].

φ [] hs [MPa] n ec0 ed0 ei0 α β
33 1000 0.25 0.95 0.55 1.05 0.25 1.0


Table 1: Basic hypoplastic parameters of Hochstetten sand.



The basic law parameters should be specified in group_materi_plasti_hypo_wolffersdorff.

Cohesion

A simplistic appraoch to include cohesion is used here. Instead of feeding the real effective stress state σij into the hypoplastic law, an alternative effective stress state σijc is used. Cohesion is modelled by subtracting in each of the normal stress components a value c representing cohesion: σ11c = σ11c, σ22c = σ22c and σ33c = σ33c. The shear stresses are not altered: σ12c = σ12, etc.

The cohesion value should be specified in group_materi_plasti_hypo_cohesion.

Intergranular strains

In order to take into account the recent deformation history, an additional tensorial state variable Sij1 is introduced.

Denoting the normalized magnitude of Sij
ρ =
SijSij
R
(R is a material parameter) and the direction of Sij
Sij =
Sij
SklSkl
(Sij=0 for Sij=0), the evolution equation for the intergranular strain tensor reads:

Sij =

( Iijkl−ρβrSij Sklkl for Sijєij > 0
єij for Sijєij ≤ 0
   ,


where Sij is the objective rate of intergranular strain. Rigid body rotations are treated elsewhere (see the section on memory). From the evolution equation (2.2.4) it follows that ρ must remain between 0 and 1.

The general stress-strain relation is now written as
σij = Mijklєkl    .


The fourth order tensor Mijkl represents the incremental stiffness and is calculated from the hypoplastic tensors Lijkl and Nij which may be modified by scalar multipliers mT and mR, depending on ρ and on the product Sijєij:

Mijkl = [ ρχ mT + (1−ρχ)mR ] Lijkl +
  +


ρχ (1−mT) Lijmn Smn Skl + ρχ Nij Skl for Sijєij > 0
ρχ (mRmT) Lijmn Smn Skl for Sijєij ≤ 0


χ is an additional material parameter.

An example intergranular parameters can be found in Tab. 2.

R mR mT βr χ
1⋅ 10−4 5.0 2.0 0.50 6.0


Table 2: Example of Intergranular hypoplastic parameters.



The intergranular parameters should be specified in group_materi_plasti_hypo_intergranularstrain. Also you need to include materi_strain_intergranular in the initialisation part.

Pressure dependent initial void ratio

You can correct the initial void ratio e0, as specified in the initial value for the history variable in the node_dof records, for the initial pressure to obtain a corrected initial void ratio e.
e
e0
= exp








σijδij
hs



n



 








See the basic law description for the parameters hs and n. The σij denotes the effective stress tensor (total stresses minus any groundflow pressure). This pressure dependent initial void ratio correction can be activated by group_materi_plasti_hypo_pressuredependentvoidratio. After the initial void ratio has been established, the development of the void ratio is governed by volumetric compression or extension of the granular skeleton.

2.2.5  Damage




In the presence of materi_damage d, the materi_stress follows:

σijdamaged = (1−d) σijundamaged





For the damage, the group_materi_damage_mazars model is available:
d = dt   αβ   +   dc   (1 − α )β


where
dt = 1. − (1−at)  
є0
єeq
at   e bteq − є0)


and
dc = 1. − (1−ac)  
є0
єeq
ac   e bteq − є0)


Here єeq contains the positive principal strains. The parameter α is given by the ratio єeq/є, where є contains the total strains (both negative and positive). The parameter є0 is the strain threshold for damage; other material parameters are β  ,  at  ,  bt  ,  ac  ,  bc. Typically for concrete:
1.e−4 < є0 < 3.e−4   ;   β = 1.    ;    1 < at < 1.5    ;    500 < bt < 2000    ;    0.7 < ac < 1.2    ;    e4 < bc < 5 e4


You can combine damage freely with plasticity models or other material behavior.

2.2.6  Average stress (hydrostatic compressibility)

An extra average stress contribution on each of σ11, σ22 and σ33 is

1
co
vi
xi


where co is the group_materi_elasti_compressibility, which should not be 0. This pressure term can e.g. be used to model nearly incompressible fluids. The compressibility contribution should be combined with a contribution for the deviatoric stresses (e.g. group_materi_viscosity).

2.2.7  Thermal stresses

Temperature rates cause fictitious thermal strain rates

− α T δij     where     δij=1     if     i=j     else     δij=0


where α is the group_materi_expansion_linear coefficient and T is the condif_temperature. These fictitious thermal strain rates in turn lead to stress rates.

2.2.8  Hyper elasticity

Hyper elasticity is used to model rubbers. It should be combined with a total lagrange formulation for the memory of the material (so use -total or -total_piola for group_materi_memory).

The stresses follow from a strain energy function (with Cij components of the matrix C, and where F is the deformation tensor and U is the stretch tensor following from the polar decomposition of the deformation tensor)

2
W
Cij


C = FT F = UT U


Deviatoric contributions

To obtain a purely deviatoric function, the following strain measures are used (with I1, I2 and I3 the first, second and third invariant of the elastic strain matrix C respectively)

J1 = I1 I3
−1
3
 
        J2 = I2 I3
−2
3
 


The group_materi_hyper_besseling function reads ( with K1, K2 and α user defined constants)

W = K1 ( J1 − 3 ) α+ K2 ( J2 − 3 )


The group_materi_hyper_mooney_blatz_ko function reads (with G and β user defined constants)

W = G * 0.5 * ( I1 − 3.0 + (2.0/β) ( J−β − 1. ) );


This Blatz-Ko hyperelastic material hardens in compression, and softens slightly in tension; it models a foamlike rubber.

The group_materi_hyper_mooney_rivlin function reads (with K1 and K2 user defined constants)

W = K1 ( J1 − 3 ) + K2 ( J2 − 3 )


The group_materi_hyper_neohookean function reads (with K1 a user defined constant)

W = K1 ( J1 − 3 )


The group_materi_hyper_reducedpolynomial function reads (with Ki user defined constants)

W = Ki ( J1 − 3 )i


where a summation over i = 1,2,… is applied.

Volumetric contributions

We define J = √ I3 . Now a volumetric part can be added to the strain energy.

The group_materi_hyper_volumetric_linear contribution reads:

W =
K
2
(J−1)2


The group_materi_hyper_volumetric_murnaghan contribution reads:

W =
K
β
(
1
β−1
J− β + 1 ) J


The group_materi_hyper_volumetric_polynomial contribution reads:

W =
Ki
2
(J−1)2i


for i=0,1,….

The group_materi_hyper_volumetric_simotaylor contribution reads:

W =
K
2
( (J−1)2 + (ln J)2 )


The group_materi_hyper_volumetric_ogden contribution reads:

W =
K
β
(
1
β
(J− β −1) + ln J )


As an example, you can combine the group_materi_hyper_mooney_rivlin energy function with the group_materi_hyper_volumetric_linear so that the total strain energy function becomes:

W = K1 ( J1 − 3 ) + K2 ( J2 − 3 ) +
K
2
(J−1)2


Here the initial shear modulus and bulk modulus are included as:

initial   shear   modulus = 2 ( K1 + K2 )


and

initial   bulk   modulus = K


respectively.

2.2.9  Viscoelasticity

Viscoelasticity is modeled with n parallel group_materi_maxwell_chains. Each of the chains contains a spring with stiffness Em in line with a dash pot with relaxation time tm (m indicates the m-th maxwell chain). The viscoelastic stress rate is given by (with Cijklm is the elastic tensor modulus of the m-th maxwell chain (depending on Em and the poisson ratio))

m=n−1
Σ
m=0
( Cijklm єklelas
σijm
tm
)


2.2.10  Viscoplasticity

Viscoplasticity is a model for rate-dependent plasticity. Rate dependent plasticity is important for (high-speed) transient plasticity calculations. It should be used in combination with a plasticity law. Viscoplasticity influences the stresses via the plastic strains.

The group_materi_plasti_visco_exponential model reads

єkl plas = γ   p   eα f   
fflow
∂ σkl


where γ and α are material fluidity constants and p is the pressure. This model was first developed for visco-plastic soil behaviour.

The group_materi_plasti_visco_power model reads

єkl plas = η (
f
fref
) p   
fflow
∂ σkl


where η (fluidity constant), fref (reference value for yield function, e.g. the yield stress for Von-Mises plasticity) and p (power) are user specified parameters.

2.2.11  Viscosity

The viscous contribution to the total stress is

2 ν Dij


where

Dij = 0.5 (
vi
xj
+
vj
xi
)


and divergence is neglected since we only model slightly compressible flows.

Viscous heat generation


The viscous energy loss is turned into heat rate per unit volume q:

q = 2 ν Dij Dij


See group_materi_viscosity_heatgeneration.

2.3  Contact analysis

2.3.1  Penalty formulation

In contact analysis, normal forces Fn follow from the condition that bodies cannot penetrate each other. Since we use a penalty formulation, the normal force is given by

Fn = λ un


where u is the penetration and λ is called the contact_penalty_velocity because its generates forces on the velocity unknowns. You can also impose groundflow_pressure and condif_temperature contact conditions by respectively the penalty factors contact_penalty_pressure and contact_penalty_temperature.

2.3.2  Friction and frictional heat generation

This normal force leads to a friction force Ff which equals

Ff = ν Fn


where ν is the friction coefficient (see contact_friction. The friction force causes heat generation rate Q:

Q = η Ff vf


where vf is the slip velocity, and the factor η is a user specified factor which determines which part of the frictional energy loss is transformed into heat (η is between 0 and 1; see contact_heatgeneration).

2.4  Ground water flow

2.4.1  Storage equation

The hydraulic pressure head h follows from the storage equation:

C   h = ( k1p
2 h
x12
+ k2p
2 h
x22
+ k3p
2 h
x32
) +
vi
xi
+ f


Primary unknown is the hydraulic pressure head groundflow_pressure. Further notation: C group_groundflow_capacity; kip group_groundflow_permeability in i-direction; xi space coordinate; vi material velocity (if present); f force_element_volume (hydraulic pressure source). The equation is given for space coordinates following material velocities vi (if present).

Groundflow velocities

The groundflow velocities, after initializing groundflow_velocity, follow from:
vig = kip
h
xi


Total groundwater pressure

The total groundwater pressure is by example needed to calculate the total stresses in soils (total soil stress = effective soil stress + total groundwater pressure). The total groundwater pressure follows from:

ptotal = h − ρ g z


where g is the gravitational acceleration, and ρ is the groundflow_density (Please notice that z typically is a negative number).

Static groundwater pressure

The static pressure due to gravity is:

pstatic = ρ g Δ z
where the Δ z is the distance to the groundwater level, the phreatic level. The phreatic level needs to be specified with the groundflow_phreaticlevel record. If that groundflow_phreaticlevel record is not specied, the static pressure part is not used, so that the static pressure becomes zero.

Dynamic groundwater pressure

The dynamic groundwater pressure follows from
pdynamic = ptotalpstatic


Boundary conditions

If the groundwater velocity is 0 normal to an edge (say at the interface with a rock layer it is zero), then you should prescribe nothing on that edge (Tochnog will then take care of that boundary condition for you).

At the phreatic level where the groundflow meets free air the hydraulic pressure head should become ρ g z. You can either set this yourself by using bounda_unknown combined with bounda_time or else demand that Tochnog automatically does it for you by activating the option groundflow_phreaticlevel_bounda.

At edges where you have some other hydraulic head you need to specify that head yourself with bounda_unknown and bounda_time records.

If gravity is not of importance, e.g. in biomechanics where the storage equation is used to model fluid transport in soft tissues, the hydraulic pressure head h is equal to the total pressure, and thus is zero at edges where the water meets the free air. In this case, set h to zero by using bounda_unknown combined with bounda_time.

Postprocessing

For all printing, plotting etc. you normally get the hydraulic pressure head h since it is the primary unknown solved in the storage equation. The total pressure, static pressure and dynamic pressure you can get by the post_calcul option.

Naming conventions

To connect better to conventional naming, we rematk that the capacity depends on the porosity n and water compressibility β:

C = n   β


and for the permeability:

kip =
ki
ρ   g


where ki is the permeability in i-direction.

Consolidation analysis

Look in the 'Consolidation' section of the 'Interaction analyzes and advanced analyzes' chapter in the end of this manual on how to perform consolidation analyzes (combined groundwater flow with soil stress analyzes).

2.5  Wave equation



s
t
= c2 (
2 s
x12
+
2 s
x22
+
2 s
x32
)
s
t
= s


The primary unknowns are the wave_scalar s and its first time derivative wave_fscalar s ( as TOCHNOG only solves first order in time equations, the first time derivative of s also becomes primary unknown in order to turn this second order in time equation into a set of first order in time equations). Further notation: x space coordinate, t time and c speed of sound.


Previous Up Next