김경종
23 KiB
parameter \kappa . Substituting from (12.4) into (12.3), and using (7.42) to express a , results in
\dot {\epsilon} = \gamma \langle \Phi (\sqrt {3} (J _ {2} ^ {\prime}) ^ {1 / 2} - \sigma_ {Y}) \rangle \sqrt {3 / 2} (J _ {2} ^ {\prime}) ^ {1 / 2} \sigma^ {\prime} = \Gamma (\sigma^ {\prime}) \sigma^ {\prime}, \tag {12.5}
in which \sigma' are the deviatoric stresses and \Gamma(\sigma') is a symmetric viscoplastic compliance matrix whose form can be explicitly determined on prescription of the appropriate flow function \Phi . Thus a relationship has been established between the total strain rate and the deviatoric stresses.
The strain rate can be expressed in terms of the displacement velocities v by taking the differential form of the standard strain/displacement relationship, to give
\dot {\epsilon} = B v. \tag {12.6}
We assume, as for the viscoplastic case of Chapter 8, that the flow velocities are sufficiently slow to neglect inertia effects and that the following standard static equilibrium equations therefore hold.
\int_ {V} \pmb {B} ^ {T} \pmb {\sigma} d V + \pmb {f} = 0, \tag {12.7}
in which f are the applied forces comprising body forces b and boundary tractions, t. Thus a complete analogy exists between the above problem and the case of an elastic material in which the relationship between stress and strain is nonlinear according to
\sigma = D (\sigma) \epsilon . \tag {12.8}
Table 12.1 Correspondence between small strain nonlinear elastic problems and viscoplastic flow situations
| Small strain nonlinear elasticity | Flow problem |
| Displacements, d | Velocities, v |
| Stresses, σ | Stresses, σ |
| Strains, ε | Strain rates, ε |
| Applied forces, f | Applied forces, f |
| Nonlinear elastic compliance matrix, [D(σ)]-1 | Viscoplastic compliance matrix, Γ(σ) |
This analogy is indicated in Table 12.1. Therefore flow problems, in which the elastic components of deformation are negligible, can be solved by use of a linear elastic computer code which includes a facility for dealing with a stress dependent D matrix. Obviously the steady state solution to the flow problem must be arrived at in an iterative manner and a similar procedure must be employed in the corresponding elastic solution. The simplest approach
is to proceed by the method of direct iteration, as described in Chapters 2 and 3, and to base the value of the compliance matrix \Gamma on the current value of \sigma . This solution procedure can be summarised as follows:
(1) From the stresses \sigma^n at iteration n evaluate the viscoplastic compliance matrix \Gamma(\sigma^n) = \Gamma^n .
(2) Compute the element stiffness matrix of each element as
\int_ {V} \boldsymbol {B} ^ {T} [ \boldsymbol {\Gamma} ^ {n} ] ^ {- 1} \boldsymbol {B} d V
and also the consistent nodal applied forces, f^{(e)} .
(3) Assemble and solve the stiffness equations to give the improved velocity estimate, v^{n+1} .
(4) Compute the strain rates, \dot{\mathbf{e}}^{n + 1} = \mathbf{B}\mathbf{v}^{n + 1} .
(5) Compute the stresses, \sigma^{n + 1} = [\Gamma^n ]^{-1}\dot{\epsilon}^{n + 1} .
(6) Return to Step 1 and repeat the process until convergence takes place (i.e. v^{n+1} \approx v^n ).
The procedure described above is most suitable when boundary and body forces produce the forcing action. For the case when the problem is defined in terms of prescribed boundary velocities the compliance matrix \Gamma must be expressed in terms of the current strain rate, \dot{\epsilon} . ^{(12)}
For metal forming problems, the situation is complicated by the fact that the geometry of the deforming solid is continually varying throughout the process. For such problems the transient form of the flow equations must be used and an incremental procedure can be adopted by which the coordinates of the finite element mesh are sequentially updated during solution. ^{(13)}
It should be noted that no volumetric strain rate exists for some viscoplastic flow laws, as generally defined by (12.3), and this is indeed the case for the Von Mises criterion employed in (12.5). Consequently the viscoplastic compliance matrix \Gamma cannot be inverted as required by Step 2 above and the same numerical difficulties that exist in incompressible elastic problems are encountered. However these can be readily overcome by the use of selective integration techniques whereby the element stiffness matrix is separated into volumetric and deviatoric components. ^{(14)} The near singularity arising in the former term as incompressible behaviour is approached is then numerically removed by employing a low order Gaussian integration rule.
An important application of the above solution process is to the flow of non-Newtonian fluids, in which the material viscosity depends nonlinearly on the shear strain rate. Practical examples of such flow can be found in Refs. 15 and 16. Deviations from Newton's law of viscosity are best illustrated by means of flow curves and some of the most important cases are shown in Fig. 12.6. The effective stress, \bar{\sigma} , and effective strain rate, \bar{\epsilon} , are defined by (7.12) and (7.22) respectively.
line
| Model | Effective Stress (σ̄) | Effective Strain Rate (ε̅) |
|---|---|---|
| Bingham plastic | γ(σ̄ - σy) | Not labeled |
| Eyring pseudoplastic | B sinh(σ̄/A) | Not labeled |
| power law pseudoplastic | γ/ε̅M⁻¹ σ̄ | Not labeled |
| ideal plasticity | γσ̅ | Not labeled |
| Newtonian ε̅ | γσ̅ | Not labeled |
Fig. 12.6 Various flow curves for non-Newtonian fluids.
The Bingham fluid is seen to be a particular form of viscoplastic relation (12.3) or (12.5). Writing in terms of the effective stress and strain rate, (12.5) can be expressed as
\bar {\sigma} = \mu \dot {\bar {\epsilon}}, \tag {12.9}
where the apparent viscosity \mu is given by
\frac {1}{\mu} = \frac {\sqrt {(3)} \gamma}{2 (J _ {2} ^ {\prime}) ^ {1 / 2}} \langle \Phi [ (\sqrt {3}) (J _ {2} ^ {\prime}) ^ {1 / 2} - \sigma_ {Y} ] \rangle . \tag {12.10}
For the Bingham plastic we can write from the expression given in Fig. 12.6 and using (12.9) that
\mu = \frac {\overline {{{\dot {\epsilon}}}} / \gamma + \sigma_ {Y}}{\overline {{{\dot {\epsilon}}}}}. \tag {12.11}
As \gamma \to \infty , ideal plasticity behaviour is approached resulting in
\mu = \frac {\sigma_ {Y}}{\dot {\bar {\epsilon}}}. \tag {12.12}
Similarly for a Power Law pseudoplastic we have from Fig. 12.6
\mu = \frac {\overline {{{\dot {\epsilon}}}} ^ {M - 1}}{\gamma}. \tag {12.13}
Thus for each case the problem again reduces to an elastic problem in which the shear modulus is dependent on the current strain rate and can be solved
by use of the analogy indicated in Table 12.1. Solution can be achieved by use of the method of direct iteration or by the Newton–Raphson process described in Chapters 2 and 3.
As an example of viscous flow analysis ^{(17)} the problem of the flow of a Bingham fluid in a cylindrical annulus is illustrated in Fig. 12.7, where the geometry and finite element mesh employed are also indicated. Steady state flow is induced parallel to the axis of the cylinder by the application of an axial pressure gradient. The finite element velocity distributions obtained by a direct iteration solution scheme are shown in Fig. 12.8 for different values of the pressure gradient. The flow velocities are in good agreement with the theoretical solution of Ref. 18.
text_image
s o d p b u g R₂ a R₁ π/8 section shown below ∂v/∂n = 0 finite element mesh v = 0 v = 0 10 parabolic elements with 2×2 gauss integration
Fig. 12.7 Flow of Bingham fluid in an annulus under an axial pressure gradient showing finite element mesh idealisation.
12.4.2 Nonlinear fracture mechanics
A class of elasto-plastic problems which require special attention is that of crack propagation in ductile materials. Figure 12.9 illustrates the types of problem which demand solution and it is seen that a geometrical singularity exists at the crack tip. The numerical techniques presented in Chapter 7 allows the elasto-plastic stress field to be determined in the vicinity of the crack tip (for Modes I and II at least) but a criterion for propagation of the crack must be established in some way.
line
| r/R₁ | v/v_AVG (f.e.m. nodal) |
|---|---|
| 1.0 | 0.0 |
| 1.1 | 0.3 |
| 1.2 | 0.6 |
| 1.3 | 0.9 |
| 1.4 | 1.1 |
| 1.5 | 1.2 |
| 1.6 | 1.2 |
| 1.7 | 1.1 |
| 1.8 | 0.9 |
| 1.9 | 0.5 |
| 2.0 | 0.0 |
Fig. 12.8 Steady state velocity profile for the problem of Fig. 12.7 for various applied pressure gradients.
For linear elastic fracture problems crack advance can be monitored by specifying a critical value of a quantity, K, termed the stress intensity factor* which characterises the stress field in the vicinity of the crack tip according to ^{(20)}
\sigma = K f (\theta) / \sqrt {(2 \pi r)} + \text { terms of order } r ^ {0}. \tag {12.14}
A separate K parameter exists for each fracture mode, designated by K_{I} , K_{II} and K_{III} respectively and they are functions only of geometry and loading conditions. A crack in any mode is then assumed to propagate when K attains a critical value K_{c} which is treated as a material parameter.
We now seek a similar criterion for elasto-plastic material behaviour. The most widely accepted principle in present use is the so-called J contour integral attributed to Rice ^{(21)} and which was originally formulated for nonlinear elastic applications. The J integral is defined to be
J = \int_ {\Gamma} \omega d y - T _ {i} \frac {d u _ {i}}{d x} d S, \tag {12.15}
for a crack aligned in the x direction. Here \Gamma is any contour from the lower crack face leading anticlockwise around the crack tip to the upper face, S is the path length around this contour and T_{i}du_{i} is the work contribution
text_image
y r θ x Γ mode I mode II mode III
Fig. 12.9 Basic modes of fracture.
of traction components T_{i} on \Gamma moving through displacements du_{i} . The term \omega is the strain energy density defined as
\omega = \int_ {0} ^ {\epsilon} \sigma_ {i j} d \epsilon_ {i j}. \tag {12.16}
The J integral is independent of the choice of path \Gamma provided that the faces of the crack are stress free.
For Mode I opening in a strain-hardening nonlinear elastic material the near tip solution for the stress, strain and displacement can be shown to be of the form ^{(22-24)}
\sigma = C \frac {1}{r ^ {1 / (N + 1)}} \sigma (\theta)
\epsilon_ {p} = C \frac {1}{r ^ {N / (N + 1)}} \epsilon (\theta)
u = C r ^ {1 / (N + 1)} u (\theta), \tag {12.17}
where
C = \left(\frac {J E}{\sigma_ {Y} {} ^ {2} I}\right) ^ {1 / (N + 1)}. \tag {12.18}
The term N is a constant which measures the strain hardening of the material, E the elastic modulus, \sigma_{Y} the stress denoting the limit of linearity and I is a tabulated constant whose value depends on N.
For loading situations, nonlinear elastic behaviour is identical to that of a material obeying the laws of 'deformation' plasticity ^{(25)} in which the current stiffness is a function only of the current state of deformation and not of the loading path by which this condition has been reached. Furthermore for monotonic loading, experience indicates that there is no significant difference between solutions obtained by use of 'deformation' theories and the incremental theory adopted in Chapter 7. By this argument it is concluded that expressions (12.17) and (12.18) are applicable to elasto-plastic solids. Consequently crack propagation in elasto-plastic materials is governed by a critical value of the J integral.
One of the difficulties of numerical fracture studies is that a reasonably accurate prediction of the stress field in the vicinity of the crack tip is required. This is a computationally expensive process for elasto-plastic problems and in some instances economies can be made by use of special crack tip elements. For example, in Mode II deformation under plastic conditions, a shear strain singularity of order 1/r develops, which has been modelled by Levy et al. ^{(26)} by coalescing two nodes of a linear quadrilateral isoparametric element and treating their displacements independently. This approach has also been employed by Rice et al. ^{(27)}
12.4.3 Coupled-field problems
The transient analysis of many engineering systems involves the formulation of the semi-discrete coupled-field equations of motion which are then solved by a time-stepping procedure. ^{(28)} Coupled-field equations involving plasticity arise in the modelling of structure–fluid interaction, soil–fluid interaction, structure–structure interaction, etc. There are two main sources of difficulty in solving such problems:
(i) The isolated fields may display quite different response characteristics which may only be analysed efficiently by different time integration algorithms and/or different time steps.
(ii) Most engineering software has been developed for the treatment of single-field problems. The term ‘partitioned transient analysis procedures’ has been used to describe methods which allow the direct time integration of the entire equations to be performed by either sequential or parallel execution of single-field analyzers.
We have discussed partitioned procedures for structural dynamic problems in Chapter 11. We described an implicit-explicit partition through which meshes that exhibit high (low) frequency response characteristics are treated by implicit (explicit) integration formulae. Park ^{(29)} has recently extended the approach described in Chapter 11.
Park et al. ^{(30)} have studied implicit–implicit partitions in certain types of fluid-structure interaction problems. The solution of these coupled-field equations was obtained by a sequential execution of fluid and structural analyzers which gave rise to the term ‘staggered solution procedures.’
Hughes ^{(31)} has summarised recent work on transient fluid-structure interaction problems. In particular he mentions work on procedures known as mixed, or arbitrary, Lagrangian–Eulerian methods.
In recent work on soil liquefaction problems, Zienkiewicz et al. ^{(32)} have devised a model which couples the soil and pore-fluid behaviour during earthquakes. Pore pressure build up and pore water migration are both accurately modelled.
Many other coupled-field problems involving elasto-plastic behaviour have been reported in the literature. It should however be emphasised that care should be taken in considering the stability of such schemes.
12.4.4 Elasto-plastic and geometrically nonlinear analyses of plates and shells
The linear and nonlinear finite element analysis of plates and shells has attracted much attention in the last decade. Two basic approaches have been adopted:
(i) The classical procedure
Here a plate or shell theory is used as a basis for the finite element formulation. Let us briefly summarise such an approach. We begin with the field equations of the three-dimensional theory and make various assumptions which lead to the plate or shell theory. In the reduction from three to two dimensions we include an analytical integration over the thickness. We then base our finite element discretisation process on the plate or shell theory. The surface geometry (in the case of shells) and the field variables are approximated using discrete nodal values and suitable interpolation functions. Integration of the various element stiffness and force terms is carried out over the reference surface. Stresses may then be obtained from the stress resultants. Examples of such an approach include the simple facet element and the many elements derived from classical thin plate theory, Mindlin/Reissner plate theory, shallow shell theory or even higher order shell theories. ^{(33,34)} There are very many examples of the application of the classical procedures in nonlinear finite element analysis of plates and shells. We include a brief sample in the list of references to this chapter. ^{(35-38)} For elasto-plastic problems many research workers express the yield function in terms of the stress resultants (cf. the non-layered approach in Chapter 9). For example, Crisfield ^{(39-44)} uses a
modified Ilyushin yield criterion expressed in terms of the bending moments [M_{x}, M_{y}, M_{xy}]^{T} and the membrane forces [N_{x}, N_{y}, N_{xy}]^{T} . To allow for the gradual spread of plasticity over the plate or shell thickness, a modified classical procedure may be adopted in which integration through the thickness is performed numerically during the finite element stiffness and force evaluation rather than analytically prior to the finite element discretisation. Gauss–Legendre, Lobatto and the mid-ordinate rules are frequently used for this purpose. To allow for geometrically nonlinear effects, total or updated Lagrangian approaches are adopted. ^{(45-55)}
(ii) Ahmad and related elements
Here isoparametric elements with independent rotational and displacement degrees of freedom are used. This concept originally introduced by Ahmad et al. ^{(56)} was later extended to allow for the linear analysis of thin as well as moderately thick shells by Zienkiewicz et al. ^{(57)} by the use of the reduced integration technique.*
Ahmad elements were originally developed because of the computational difficulties encountered in the use of the usual three-dimensional elements for the analysis of plates and shells. In the three-dimensional elements the stiffness coefficients corresponding to the transverse displacement degrees of freedom are very much larger than those corresponding to the longitudinal displacements. Erroneous strain energy corresponding to the normal stresses in the thickness direction are also introduced. Both of these difficulties are overcome in Ahmad elements. Normals to the plate or shell reference surface before deformation are assumed to remain straight but not necessarily normal to the reference surface after deformation. Furthermore, the normal stresses in the direction of the shell thickness are ignored and suitably modified constitutive equations are adopted.
Various nonlinear problems have been solved using Ahmad shell elements by Ramm ^{(67)} , Krakeland ^{(68)} , Bathe and Bolourchi ^{(69)} and others ^{(70-73)} . As in the modified classical procedures, to allow for the gradual spread of plasticity over the plate or shell thickness, numerical integration techniques are adopted. For geometrically nonlinear behaviour both total and updated
- The Mindlin plate elements described in Chapters 6 and 9 are simply plate versions of the Ahmad elements in which integration has been carried out analytically through the plate thickness. Much work on reduced and selective integration techniques ^{(58-65)} eventually led to the recognition that the use of selective integration techniques is equivalent to the use of a special type of mixed formulation. ^{(66)} Defects in the Ahmad elements have now been widely acknowledged and the use of the 9-node heterosis Mindlin plate element and the 16-node cubic Ahmad element are usually recommended. Other Ahmad/Mindlin C(0) elements should be used with caution as they are known to give overstiff solutions for thin plates and shells and to develop mechanisms (zero energy modes) or near mechanisms (artificially low energy modes) when reduced or selective integration is used.
Lagrangian schemes have been used. Special techniques have been incorporated to allow for large rotations in the total Lagrangian formulations. ^{(67-69)}
The Ahmad shell concept has been developed further by its originator Irons with the introduction of the Semiloof element. ^{(90)} Irons adopted a convenient nodal configuration involving rotational degrees of freedom at 'Loof' nodes on the curved boundaries of the element. By imposing a series of constraints to eliminate transverse shear effects (reminiscent of the discrete Kirchhoff hypothesis), a highly effective thin shell element is obtained. Various research workers ^{(74-76)} have successfully extended this work into the nonlinear range.
Both classical and Ahmad procedures may be used as a basis for the nonlinear analysis of reinforced concrete plates and shells using the layering concept described in Chapter 9. Special constitutive relationships are required to represent the concrete and steel reinforcing bars are treated as a 'smeared' layer with uni-directional elasto-plastic properties. Much work has been completed in this area. ^{(77-85)}
Elasto-viscoplastic plates and shells are easily developed using the concepts described in Chapters 8 and 9. ^{(86,87)}
12.5 Equation solving techniques
12.5.1 Standard and modified Newton method
Before considering some alternative nonlinear solution procedures which may be used in elastoplastic finite element analysis we review the techniques described earlier.
As we have already seen, most elasto-plastic finite element programs are simply extensions of elastic finite element programs with linearised load increments. Some form of iterative procedure is usually adopted to dissipate the out-of-balance nodal forces.
The standard and variety of modified Newton methods were described earlier in Part I. Recall that the standard Newton method involves iterations in which
\boldsymbol {K} ^ {(i)} \left[ \boldsymbol {d} ^ {(i + 1)} - \boldsymbol {d} ^ {(i)} \right] = \psi \left(\boldsymbol {d} ^ {(i)}\right), \tag {12.19}
where d is the vector of nodal displacements and the equations \psi(d)=0 express a force balance (internal forces = external forces; either for an increment of loading or for the whole applied load). The matrix K in the standard Newton method is the Jacobian of \psi ; which is the tangential stiffness matrix K_{T}=[\partial\psi(d^{(t)}/\partial d] evaluated at the displacements described by d^{(t)} .
The modified Newton method works with a variety of approximations to K, the most simple of which is the initial elastic stiffness matrix K_{0} evaluated at the first iteration of the first load increment.