Category WATER IN ROAD STRUCTURES

The phenomena considered in this section are much more complex as they associate multiphase fluid flow and hydro-mechanical coupling (cf. relevant sub-sections of Section 11.3.2) as well as temperature effects. All the features described in the pre­ceding section are to be considered here, associated with some new points.

Heat diffusion has to be modelled. Temperature variation affects fluid flow, by a modification of the fluid specific weight or viscosity. Moreover, if the two fluids concerned are a liquid and a gas (e. g. water and air), then equilibrium between the phases has to be modelled: dry air – vapour equilibrium.

Heat transfer is governed not only by conduction but also by advection by the liquid and gas movements. Similarly, transfers of vapour and dry air in the gas phase are not only governed by diffusion and gradient of species density, but also by advection of the global gas movements.

In such an analysis, the total number of d. o.f. per node is 5 for a 2D problem: 2 displacements, 2 fluid pore pressures and the temperature.

The coupled flow of two, different, fluids in a partly saturated rigid media is now considered. Unsaturated soils provide a common example, where the two fluids are water and air. Often, the air phase is considered to be at constant pressure, which is generally a relevant approximation as air pressure doesn’t highly affect the water flow. In such an arrangement only one d. o.f. per node is sufficient, and the classical diffusion equations (see Eq. 6.10) are relevant, with parameters depending on the suction or saturation level.

Mixing between different fluids is sometimes possible. Then two or more d. o.f. per node are to be considered. The permeability and storage equations of each phase depend on the suction or saturation level, and so the problem may be highly non­linear. However, it is not difficult to numerically develop coupling, as the formula­tion is similar for each phase.

11.3.2.2 Diffusion and Transport Coupling

Heat and single fluid flow in a rigid porous media are considered here. The fluid specific weight and viscosity are dependant on the temperature or salt concentration, and the heat or salt transport by advection-diffusion process is dependant on the fluid flow. Then a diffusion process and an advection-diffusion process have to be solved simultaneously. In this case there are 2 d. o.f. per node: fluid pore pressure and salt concentration or temperature.

A large number of different phenomena may be coupled. It is impossible to discuss here all potential terms of coupling, and we will restrict ourselves to some basic cases often implied in environmental geomaterial mechanics. In the following para­graphs, some fundamental aspects of potential coupling are briefly described.

11.3.2.1 Hydro-Mechanical Coupling

In the case of hydro-mechanical coupling, the number of d. o.f. per node will be 3 (2 displacements +1 pore pressure) for 2D analysis and 4 (3 displacements +1 pore pressure) for 3D analysis.

Coupling mechanical deformation of soils or rock mass and water flow in pores is a frequent problem in geomechanics. The first coupling terms are related to the influence of pore pressure on mechanical equilibrium through Terzaghi’s postulate (or through any other effective stress concept or net stress use, cf. Eq. 9.20):

a = a’ + uI

with the effective stress tensor a’ related to the strain rate tensor thanks to the con­stitutive Eq. 11.3, and the identity tensor I.

The second type of coupling concerns the influence of the solid mechanics be­haviour on the flow process, which comes first through the storage term. Storage of water in saturated media is mainly due to pores strains, i. e. to volumetric changes in the soil/rock matrix:

Another effect, which may be considered, is the permeability change related to the pore volume change, which may, for example, be modelled by the Kozeny – Carman law as a function of the porosity K = K (n).

Biot proposed an alternative formulation for rocks where contacts between grains are much more important than in soils. Following Biot, the coupling between flow and solid mechanics are much more important (Detournay & Cheng, 1991; Thimus etal., 1998).

The time dimension may cause some problems. First, an implicit scheme is used for the solid mechanics equilibrium and various solutions are possible for the pore pressure diffusion process. Consistency would imply use of fully explicit schemes for the two problems. Moreover, it has been shown that time oscillations of the pore pressure may occur for other time schemes. Associated to Terzaghi’s postulate, oscillations could appear also on the stress tensor, which can degrade the numerical convergence rate for elasto-plastic constitutive laws.

When using isoparametric finite elements, the shape functions for geometry and for pore pressure are identical. Let us consider, for example, a second or­der finite element. As the displacement field is of second order, the strain rate field is linear. For an elastic material, the effective stress tensor rate is then also linear. However, the pore pressure field is quadratic. Then Terzaghi’s postulate mixes the linear and quadratic fields, which is not very consistent. Some authors have then proposed to mix in one element quadratic shape functions for the ge­ometry and linear shape functions for pore pressure. But then problems arrive with the choice of spatial integration points (e. g. should there be 1 or 4 Gauss points?).

Numerical locking problems may also appear for isoparametric finite elements when the two phase material (water plus soil) is quite incompressible, i. e. for very short time steps with respect to the fluid diffusion time scale. Specific elements have to be developed for such problems.

11.3.1 Finite Element Modelling: Monolithical Approach

Modelling the coupling between different phenomena should imply the need to model each of them and, simultaneously, all the interactions between them. A first approach consists in developing new finite element and constitutive laws especially dedicated to the physical coupled problem to be modelled. This approach allows taking accurately all the coupling terms into account. However there are some draw­backs that will be discussed in Section 11.3.4. Constitutive equations for coupled phenomena will be discussed in the following sections.

The number of basic unknowns and, consequently, the number of degrees of freedom – d. o.f. – per node are increased. This has a direct effect on the computer time used for solving the equation system (up to the third power of the total d. o.f. number). Coupled problems are highly time consuming.

Isoparametric finite elements will often be considered. However, some specific difficulties may be encountered for specific problems. Nodal forces or fluxes are computed in the same way as for decoupled problems. However, the stiffness matrix evaluation is much more complex, as interactions between the different phenomena are to be taken into account. Remember that the stiffness or iteration matrix, Eq. 11.25, is the derivative of internal nodal forces/fluxes with respect to the nodal unknowns (displacements/pressures/etc…). The complexity is illustrated by the following scheme of the stiffness matrix, restricted to the coupling between two problems.

The part of the stiffness matrix in cells 1-1 and 2-2 are similar, or simpler, than the ones involved in uncoupled problems. The two other cells, 1-2 and 2-1, are new and may be of a greater complexity. Remember also that the derivative considers internal nodal forces/fluxes as obtained numerically, i. e. taking into account all numerical integration/derivation procedures. On the other hand, the large difference of orders of magnitude between different terms may cause troubles in solving the problem and so needs to be checked.

Numerical convergence of the Newton-Raphson process has to be evaluated care­fully. It is generally based on some norms of the out-of-balance forces/fluxes. How­ever, coupling often implies the mixing of different kinds of d. o.f., which may not be compared without precaution. Convergence has to be obtained for each basic problem modelled, not only for one, which would then predominate in the computed indicator.

Let us first consider a purely advective process. In this case, the transport is gov­erned by the advection Eq. 11.11 and by the balance Eq. 11.6. Associating these two equations, one obtains:

(VTC). ffluld + C = 0 (11.45)

diff v ‘

which is a hyperbolic differential equation. It cannot be solved by the finite element or finite difference problem, but by characteristic methods. The idea is to follow the movement of a pollutant particle by simply integrating step by step the fluid velocity field. This integration has to be accurate enough, as errors are cumulated from one step to the next.

On the other hand, if advection is very small compared to diffusion, then the finite element and finite difference methods are really efficient.

For most practical cases, an intermediate situation holds. It can be checked by Peclet’s number, Eq. 11.13, which is high for mainly advective processes and low for mainly diffusive ones. As diffusion has to be taken into account, the numeri­cal solution must be based on the finite element method (the finite difference one may also be used but will not be discussed here). However, numerical experiments show that the classical Galerkin’s formulation gives very poor results with high spatial oscillations and artificial dispersion. Thus, new solutions have been proposed (Zienkiewicz & Taylor, 1989, Charlier & Radu, 2001). A first solution is based on the use in the weighted residual method of a weighting function that differs from the shape one by an upwind term, i. e. a term depending in amplitude and direction on the fluid velocity field. The main advantage of this method is the maintenance of the finite element code formalism. However, it is never possible to obtain a highly accurate procedure. Numerical dispersion will always occur.

Other solutions are based on the association of the characteristic method for the advection part of the process and of the finite element method for the diffusive part (Li et al., 1997). The characteristic method may be embedded in the finite element code, which has a strong influence on the finite element code structure. It is also possible to manage the two methods in separated codes, as in a staggered procedure (cf. Section 11.3.4).

The theoretical analysis of a time integration scheme accuracy and stability is generally based on a simplified problem (Zienkiewicz et al., 1988). Let us con­sider diffusion phenomena restricted to the linear case. Introducing the discre – tised field (Eq. 11.15) into the constitutive equations gives Darcy law (Eq. 11.7) (neglecting here the gravity term for the sake of simplicity and using the more gen­eral pressure p in place of u) in the following form:

fi = – K dtp = – K (dtNL )pL (11.32)

p p

Similarly the storage law (linear case) can be re-written:

S = rp = rNLpL (11.33)

where r is a storage parameter (cf. Eq. 11.8). Neglecting source terms, the weak form of the balance equation, Eq. 11.21, then produces:

j [SSp – ft dt (8p)]dV

V

к

– Nl pl diNK 8 pKdV = 0

p

VV

Considering that nodal values are not affected by the integration, this becomes:

(11.35)

which is valid for any arbitrary perturbation 8p. Thus:

RklP l + KklPl = 0 (11.36)

which is a simple system of linear equations with a time derivative, a storage matrix R (of which Rkl is an element) and a permeability matrix K (of which KKL is an element). One can extract the eigenvalues of this system and so arrive at a series of scalar independent equations of similar form:

pL + a2LpL = 0 (no summation) (11.37)

where L represents now the number of the eigenmode with the eigenvalue aL al­though it will be omitted in the following. The exact solution for Eq. 11.35 is a decreasing exponential function of time, t:

p(t) = p(t0)e-2t (11.38)

This problem represents the damping of a perturbation for a given eigenmode. Numerically, the modelling is approximated and numerical errors always appear. If Eq. 11.36 is well modelled, any numerical error will be rapidly damped, if the error source is not maintained. Following this analysis, the whole accuracy and stability discussion may be based on these last scalar Eqs. 11.37 and 11.38.

Introducing the time discretisation (Eqs. 11.26 and 11.27) into Eq. 11.37 gives:

PB-PA + a2 [(1 – r)Pa + трв] = 0 (11.39)

which allows evaluation of the end-of-step pressure as a function of the pressure at the beginning of the step:

Pb = UpA

with the amplification factor, U:

1 – (1 – t) a2 At
1 + та2 At

To ensure the damping process of the numerical algorithm, which is the sta­bility condition, it is strictly necessary that the amplification factor remains lower than unity:

-1 < U < 1

This condition is always verified if т > У2, and conditionally satisfied otherwise:

This last equation is not easy to verify, as it depends on the eigenvalues, which are generally not computed. Therefore, for classical diffusion processes considered in geomaterials, the condition т > /2 is generally used.

It should be noted that the amplification factor becomes negative for large time steps, except for the fully implicit scheme. Then, the perturbed pressure de­creases monotonically in amplitude but with changes of sign. This may be ques­tionable for some coupled phenomena, as it could induce oscillation of the coupled problem.

Let us now consider the accuracy of the numerical schemes. Developing in Tay­lor’s series the exact and numerical solution allows a comparison:

1 2 1 3

Uexact = 1 – X + X – X + …

2 6

^numerical = 1 – X + вX2 – в2X3 + …X = a2At (11.44)

It appears that only the Crank-Nicholson scheme т = /2 has second order accu­racy properties. However this conclusion is limited to infinitesimal time steps. For larger time steps, as in most numerical models, the Galerkin’s scheme т = 2 /3 gives the optimal compromise and should generally be used.

The whole discussion related to the stability and accuracy of the proposed time numerical schemes was based on eigenmodes of a linear problem. Can we extrapo­late them to general problems? The eigenvalue passage is only a mathematical tool allowing consideration of scalar problems, and has no influence on our conclusions. Conversely, the non-linear aspects could sometimes modify our conclusions. How­ever, it is impossible to develop the analysis for a general non-linear problem, and the preceding conclusions should be adopted as guidelines, as they appear to be fruitful in most cases.

For solid mechanics problems, the constitutive law form (Eqs. 11.3 and 11.4) is an incremental one and differs from the ones for diffusion problems (Eq. 11.7). The knowledge of the stress tensor at any time implies that a time-integrated constitutive law is required. The stress tensor is a state variable that is stored and transmitted from step to step based on its final/initial value, and this value plays a key role in the numerical algorithm.

Then, in nearly all finite element codes devoted to modelling, equilibrium is ex­pressed at the end of the time steps, following a fully implicit scheme (t = 1), and using the end of step stress tensor value.

However, integrating the stress history with enough accuracy is crucial for the numerical process stability and global accuracy. Integration of the first order differ­ential equation (Eq. 11.4):

tB

ctb = cta + J Eeps dt (11.31)

tA

can be based on similar concepts as the one described in the preceding paragraph (the superscripts of ct here indicating the time at which ct is evaluated). Various time schemes based on different т values may be used for which similar reflections on stability and accuracy can be made.

When performing large time steps, obtaining enough accuracy can require the use of sub-stepping: within each global time step (as regulated by the global numer­ical convergence and accuracy problem) the stress integration is performed at each finite element integration point after division of the step into a number a sub-steps allowing higher accuracy and stability.

The time dimension appears in the form of a first order time derivative in the consti­tutive mechanical model (Eqs. 11.3, 11.4) and in the diffusion problems though the storage term (Eq. 11.6). We will here discuss the time integration procedure and the accuracy and stability problems that are involved.

11.2.6.1 Time Integration – Diffusion Problems

The period to be considered is divided into time steps. Linear development of the basic variable with respect to the time is generally considered within a time step:

t — tA tB — t

p = – PB + – B PA (11.26)

tB — tA tB — tA

where the subscripts A, B denote, respectively, the beginning and the end of a time step. Then the pressure rate is:

dp _ Pb – Pa _ A p dt tB — tA At

This time discretisation is equivalent to a finite difference scheme. It allows the evaluation of any variable at any time within a time step.

The balance equation should ideally be satisfied at any time during any time step. Of course this is not possible for a discretised problem. Only a mean assessment of the balance equation can be obtained. Weighted residual formulations have been proposed in a similar way as for finite elements (Zienkiewicz et al., 1988). However, the implementation complexity is too high with respect to the accuracy. Then the easiest solution is to assess only the balance equation at a given time, denoted tT, inside the time step tA to tB, such that a time variable, т, is defined:

tT tA

T =

tB — tA

All variables have then to be evaluated at the reference time, tT. Different classi­cal schemes have been discussed for some decades:

• Fully explicit scheme – t = 0: all variables and the balance are expressed at the time step beginning, where everything is known (from the solution of the preceding time step). The solution is, therefore, very easily obtained.

• Crank-Nicholson scheme or mid-point scheme – t = 1/2

• Galerkin’s scheme – t = 2/3

• Fully implicit scheme – t = 1

The last three schemes are functions of the pore pressure/temperature/concentration at the end of the time step, and may need to be solved iteratively if non-linear prob­lems are considered.

For some problems, phase changes, or similar large variations of properties, may occur abruptly. For example, icing or vaporising of water is associated with latent heat consumption and abrupt change of specific heat and thermal conductivity. Such rapid change is not easy to model. The change in specific heat may be smoothed using an enthalpy formulation, because enthalpy, H, is an integral of the specific heat, c. Then the finite difference of the enthalpy evaluated over the whole time step gives a mean value, c, and so allows an accurate balance equation:

From Eq. 11.24, it appears that the stiffness matrix is a derivative of the internal forces:

Two contributions will be obtained (Fig. 11.2). On the one hand, one has to derive the stress state with respect to the strain field, itself depending on the displacement field. On the other hand, the integral is performed on the volume, and the B ma­trix depends on the geometry. If we are concerned with large strains and if we are using the Cauchy’s stresses, geometry is defined in the current configuration, which is changing from step to step, and even from one iteration to the other. These two contributions, the material one, issued from the constitutive model, and the geometric one, have to be accurately computed in order to guarantee the quadratic convergence rate.

A similar discussion may be given for diffusive problems. However, the geom­etry is not modified for pure diffuse problems, so only the material term is to be considered.

Let us now concentrate on the finite element method. The fundamental equation to be solved is the equilibrium Eq. 11.1 (or the balance Eq. 11.6 for diffusion phenomena). As the numerical methods give an approximate solution, the equilib – rium/balance equation has to be solved with the best compromise. This is obtained by a global weak form of the local equation. Using weighted residuals, for solid mechanics, one obtains:

J [v, j Sel}]dV = J P, Sl, dV + j Р8Ш (11.17)

V V A

And for diffusion phenomena:

j [,SSp – f, d, (5p)]dV = j QSpdV + j qSpdA (11.18)

V V A

where p and q are surface terms of imposed loads/fluxes. The weighting functions are denoted 81 and 8p, and $£ represents a derivative of the weighting function based on the Cauchy’s strain derivate operator. An equivalent equation could be obtained based on the virtual power principle. The 81 and 8p terms would then be interpreted as virtual arbitrary displacements and pressures. Within the finite ele­ment method, the global equilibrium/balance equation will be verified for a number of fundamental cases equivalent to the degrees of freedom (d. o.f.) of the problem, i. e. the number of nodes times the number of degrees of freedom per node, minus the number of imposed values. The corresponding weighting functions will have simple forms based on the element shape functions.[26]

Giving a field of stress or of flux, using the weighting functions, one will obtain a value for each d. o.f., which is equivalent to a nodal expression of the equilib – rium/balance equation.

More precisely, for solid mechanics problems, one will obtain internal forces equivalent to stresses at each node, L:

Fu = f (JijBLjdV (11.19)

V

where BLj is a member of the matrix, B, of derivatives of the shape functions, N. If equilibrium is maintained from the discretised point of view, these internal forces are equal to external forces (if external forces are distributed, a weighting is necessary):

Fun = FLX (11.20)

Similarly, for diffusion phenomena the nodal internal fluxes are equivalent to the local fluxes:

Fnt = f [SNl – f diNu]dV (11.21)

V

If the balance equation is respected from the discretised point of view, these internal fluxes are equal to external ones:

FLnt = F[xt (11.22)

However, as we are considering non linear-problems, equilibrium/balance can­not be obtained immediately, but requires iteration. This means that the equations (Eqs. 11.20 and 11.22) are not fulfilled until the last iteration of each step.

Non-linear problems have been solved for some decades, and different methods have been used. From the present point of view, the Newton-Raphson method is the
reference method and probably the best one for a large number of problems. Let us describe the method. In Eq. 11.20 the internal forces Fg/ are dependant on the basic unknown of the problem, i. e. the displacement field. Similarly in Eq. 11.22 the internal fluxes are dependant on the pressure (temperature, concentration…) field.

If the external forces/fluxes don’t equilibrate, the question to be treated can be formulated in the following manner. Following the Newton-Raphson method, one develops the internal force as a first order Taylor’s series around the last approxima­tion of the displacement field:

– Fint

Fu – FL (h)) + + O2 = FLx, (11.23)

-/Kj

where the subscript (i) indicates the iteration number and O2 represents second order, infinitely small terms. This is a linearization of the non-linear equilibrium equation. It allows one to obtain a correction of the displacement field:

/А Fint 1

-g FLfQa)) – FUX) – ELikj {F^Qi)) – FUT) (11.24)

Here, the matrix, E, represented here by its member term ELi, Kj, is the so-called stiffness matrix. With the corrected displacement field, one may evaluate new strain rates, new stress rates, and new improved internal forces. Equilibrium should then be improved.

The same meaning may be developed for diffusion problems using Taylor’s development of the internal fluxes with respect to the pressure/temperature/concen – tration nodal unknowns.

The iterative process may be summarised as shown in Fig. 11.1 for a one-d. o.f. solid mechanics problem. Starting from a first approximation of the displacement field /(i) the internal forces Fint (1) (point A(1) in the figure) are computed to be lower then the imposed external forces Fext. Equilibrium is then not achieved and a new approximation of the displacement field is sought. The tangent stiffness matrix is

evaluated and an improved displacement is obtained l(2) (point B(i)) (the target being as in Eq. 11.22. One computes again the internal forces Fint(2) (point A(2)) that are again lower then the external forces Fext. As equilibrium is not yet fulfilled, a new approximation of the displacement field is sought, l(3) (point B(2)). The procedure has to be repeated until the equilibrium/balance equation is fulfilled with a given accuracy (numerical convergence norm). The process has a quadratic convergence, which is generally considered as the optimum numerical solution.

However the Newton-Raphson method has an important drawback: it needs a large amount of work to be performed as well as to be run on a computer. The stiffness matrix, E, is especially time-consuming for analytical development and for numerical inversion. Therefore other methods have been proposed:

• An approximate stiffness matrix, in which some non-linear terms are neglected.

• Successive use of the same stiffness matrix avoiding new computation and inver­sion at each iteration.

It should be noted that each alternative is reducing the numerical convergence rate. For some highly non-linear problems, the convergence may be lost, and then no numerical solution will be obtained.

Some other authors, considering the properties and the efficiency of explicit time schemes for rapid dynamic problems (e. g. for shock modelling) add an artificial mass to the problem in order to solve it as a quick dynamic one. It should be clear that such a technique might degrade the accuracy of the solution, as artificial inertial effects are added and the static equilibrium Eq. 11.1 is not checked.