Читать книгу Computational Geomechanics - Manuel Pastor - Страница 41

To complete the numerical solution, it is necessary to integrate the ordinary differential Equations (3.23), (3.27), and (3.28) in time by one of the many available schemes. Although there are many multistep methods available (see, e.g., Wood 1990), they are inconvenient as most of them are not self‐starting and it is more difficult to incorporate restart facilities which are required frequently in practical analyses. On the other hand, the single‐step methods handle each step separately and there is no particular change in the algorithm for such restart requirements.

Two similar, but distinct, families of single‐step methods evolved separately. One is based on the finite element and weighted residual concept in the time domain and the other based on a generalization of the Newmark or finite difference approach. The former is known as the SSpj – Single Step p^th order scheme for j^th order differential equation (p ≥ j). This was introduced by Zienkiewicz et al. (1980b, 1984) and extensively investigated by Wood (1984a, 1984b, 1985a, 1985b). The SSpj scheme has been used successfully in SWANDYNE‐I (Chan, 1988). The later method, which was adopted in SWANDYNE‐II (Chan 1995) was an extension to the original work of Newmark (1959) and is called Beta‐m method by Katona (1985) and renamed the Generalized Newmark (GNpj) method by Katona and Zienkiewicz (1985). Both methods have similar or identical stability characteristics. For the SSpj, no initial condition, e.g. acceleration in dynamical problems, or higher time derivatives are required. On the other hand, however, all quantities in the GNpj method are defined at a discrete time station, thus making transfer of such quantities between the two equations easier to handle. Here we shall use the later (GNpj) method, exclusively, due to its simplicity.

In all time‐stepping schemes, we shall write a recurrence relation linking a known value ϕ_n (which can either be the displacement or the pore water pressure), and its derivatives , ,… at time station t_n with the values of Φ_n+1, , ,…, which are valid at time t_n + Δt and are the unknowns. Before treating the ordinary differential equation system (3.23), (3.27), and (3.28), we shall illustrate the time‐stepping scheme on the simple example of (3.6) by adding a forcing term:

(3.34)

From the initial conditions, we have the known values of Φ_n, . We assume that the above equation has to be satisfied at each discrete time, i.e. t_n and t_n+1. We can thus write:

(3.35)

and

(3.36)

From the first equation, the value of the acceleration at time t_n can be found and this solution is required if the initial conditions are different from zero.

The link between the successive values is provided by a truncated series expansion taken in the simplest case as GN22 as Equation (3.34) is a second‐order differential equation j and the minimum order of the scheme required is then two: as (p ≥ j)

(3.37)

Alternatively, a higher order scheme can be chosen such as GN32 and we shall have:

(3.38)

In this case, an extra set of equations is required to obtain the value of the highest time derivatives. This is provided by differentiating (3.35) and (3.36).

(3.39)

and

(3.40)

In the above equations, the only unknown is the incremental value of the highest derivative and this can be readily solved for.

Returning to the set of ordinary differential equations we are considering here, i.e. (3.23), (3.27), and (3.28) and writing (3.23) and (3.28) at the time station t_n+1, we have:

(3.41)

(3.42)

assuming that (3.27) is satisfied.

Using GN22 for the displacement parameters and GN11 for the pore pressure parameter ^w, we write:

(3.43a)

and

(3.43b)

where and are as yet undetermined quantities. The parameters β₁, β₂, and are usually chosen in the range of 0 to 1. For β₂ = 0 and = 0, we shall have an explicit form if both the mass and damping matrices are diagonal. If the damping matrix is non‐diagonal, an explicit scheme can still be achieved with β₁ = 0, thus eliminating the contribution of the damping matrix. The well‐known central difference scheme is recovered from (3.41) if β₁ = 1/2, β₂ = 0 and this form with an explicit and implicit scheme has been considered in detail by Zienkiewicz et al. (1982) and Leung (1984). However, such schemes are only conditionally stable and for unconditional stability of the recurrence scheme, we require

The optimal choice of these values is a matter of computational convenience, the discussion of which can be found in literature. In practice, if the higher order accurate “trapezoidal” scheme is chosen with β₂ = β₁ = 1/2 and = 1/2, numerical oscillation may occur if no physical damping is present. Usually, some algorithmic (numerical) damping is introduced by using such values as

or

Dewoolkar (1996), using the computer program SWANDYNE II in the modelling of a free‐standing retaining wall, reported that the first set of parameters led to excessive algorithmic damping as compared to the physical centrifuge results. Therefore, the second set was used and gave very good comparisons. However, in cases involving soil, the physical damping (viscous or hysteretic) is much more significant than the algorithmic damping introduced by the time‐stepping parameters and the use of either sets of parameters leads to similar results.

Inserting the relationships (3.43) into Equations (3.41) and (3.42) yields a general nonlinear equation set in which only and remain as unknowns.

This set can be written as

(3.44a)

(3.44b)

where and can be evaluated explicitly from the information available at time t_n and

(3.45)

In this, must be evaluated by integrating (3.27) as the solution proceeds. The values of _n+1 and _n+1 at the time t_n+1 are evaluated by Equation (3.43).

The equation will generally need to be solved by a convergent, iterative process using some form of Newton–Raphson procedure typically written as

(3.46a)

where 1 is the iteration number and

The Jacobian matrix can be written as

(3.47)

where

which are the well‐known expressions for tangent stiffness matrix. The underlined term corresponds to the “initial stress” matrix evaluated in the current configuration as a result of stress rotation defined in (2.5).

Two points should be made here:

1 that in the linear case, a single “iteration” solves the problem exactly

2 that the matrix can be made symmetric by a simple scalar multiplication of the second row (provided KT is itself symmetric).

In practice, it is found that the use of various approximations of the matrix J is advantageous such as, for instance, the use of “secant” updates (see, for instance, Crisfield (1979), Matthies, and Strang (1979) and Zienkiewicz et al (2013).

A particularly economical computation form is given by choosing β₂ = 0 and representing matrix M in a diagonal form. This explicit procedure was first used by Leung (1984) and Zienkiewicz et al. (1980a). It is, however, only conditionally stable and is efficient only for phenomena of short duration.

The process of the time‐domain solution of (3.44) can be amended to that of successive separate solutions of the time equations for variables and , respectively, using an approximation for the remaining variable. Such staggered procedures, if stable, can be extremely economical as shown by Park and Felippa (1983) but the particular system of equations presented here needs stabilization. This was first achieved by Park (1983) and, later, a more effective form was introduced by Zienkiewicz et al. (1988).

Special cases of solution are incorporated in the general solution scheme presented here without any modification and indeed without loss of computational efficiency.

Thus, for static or quasi‐static, problems, it is merely necessary to put M = 0 and immediately the transient consolidation equation is available. Here time is still real and we have omitted only the inertia effects (although with implicit schemes, this a priori assumption is not necessary and inertia effects will simply appear as negligible without any substantial increase of computation). In pure statics, the time variable is still retained but is then purely an artificial variable allowing load incrementation.

In static or dynamic undrained analysis, the permeability (and compressibility) matrices are set to zero, i.e. H · f ⁽²⁾ = 0, and usually S = 0 resulting in a zero‐matrix diagonal term in the Jacobian matrix of Equation (3.47).

The matrix to be solved in such a limiting case is identical to that used frequently in the solution of problems of incompressible elasticity or fluid mechanics and, in such studies, places limitations on the approximating functions N ^u and N ^p used in (3.19) if the Babuska–Brezzi (Babuska 1971, 1973, Brezzi 1974) convergence conditions or their equivalent (Zienkiewicz et al. 1986b) are to be satisfied. Until now, we have not referred to any particular element form, and, indeed, a wide choice is available to the user if the limiting (undrained) condition is never imposed. Due to the presence of first derivatives in space in all the equations, it is necessary to use “C₀‐continuous” interpolation functions and Figure 3.2 shows some elements incorporated in the formulation. The form of most of the elements used satisfies the necessary convergence criteria of the undrained limit (Zienkiewicz 1984). Though the bi‐linear u and p quadrilateral does not, it is, however, useful when the permeability is sufficiently large.

Figure 3.2 Elements used for coupled analysis, displacement (u) and pressure (p) formulation: (a) (i) quadratic for u; (ii) linear for p; (b) (i) biquadratic for u; (ii) bilinear for p: (c) (i) linear for u; (ii) linear for p: (d) (i) linear (with cubic bubble) for u; (ii) linear for Element (c) is not fully acceptable at incompressible–undrained limits.

We shall return to this problem in Chapter 5 where a modification is introduced allowing the same interpolations to be used for both u and p.

In that chapter, we shall discuss a possible amendment to the code permitting the use of identical u ‐p interpolation even in incompressible cases.

We note that all computations start from known values of and ^w possibly obtained as the result of static computations by the same program in a manner which will be explained in the next Section 3.2.4. The incremental computation allows the various parameters, dependent on the solution history, to be updated.

Thus, for known Δ increment, the Δσ ^″ are evaluated by using an appropriate tangent matrix D and an appropriate stress integration scheme.

Further, we note that if p_w ≥ 0 (full saturation, described in Section 2.2 of the previous chapter), then we have

and the permeability remains at its saturated value

However, when negative pressures are reached, i.e. when p_w < 0, the values of S_w, χ_w, and k have to be determined from appropriate formulae or graphs.