김경종
17 KiB
text_image
10m 20m K = 2(1 + 2v) K = 1 + v
Fig. 3.10 Groundwater flow example—Problem 3.3.
text_image
φ = 100 K = 1 - 10 dφ/dr 1000m Well point
Fig. 3.11 Water flow in a confined aquifer—Problem 3.6.
The circular region shown in Fig. 3.11 has a central well point at which water is being extracted at a rate of 200 \, m^{3}/day . Determine the steady state potential distribution for this system assuming the material permeability to be nonlinear in the manner shown.
3.7 The relationship between stress, \sigma , and strain, \epsilon , for a certain locking material is given by the relationship
\sigma = \frac {E _ {0} \epsilon}{\epsilon_ {L} (\epsilon_ {L} - \epsilon)}, \tag {3.43}
in which E_{0} is the elastic modulus and \epsilon_{L} is the limiting strain value of the material. Implement this relation in program NONLAS documented in Appendix I, Section A1.3, by modifying the strain derivative function in Section 3.11.1. Also allow the behaviour of certain elements to be linear elastic. Use this modified program to determine the force displacement/relationship of the central node in Fig. 3.12 for a total applied load of 100 units.
text_image
100 Cross-sectional area, A = 1 for both members 2 2 Nonlinear elastic E₀ = 1000 εₗ = 0.1 Linear elastic material, E = 1000
Fig. 3.12 Nonlinear elastic example—Problem 3.7.
3.8 Use program ELPLAS, for which user instructions are provided in Appendix I, Section A1.4, to solve the one-dimensional elasto-plastic problem shown in Fig. 3.13.
3.9 Develop the elastic stiffness matrix, K^{(e)} , for a two-node finite element in the form of a thin disc of thickness t which is to be subjected to axisymmetric in-plane loading. Assume a linear variation between nodes, as shown in Fig. 2.7, and note the following relationships
\epsilon_ {r} = \frac {d u}{d r} = \frac {1}{E} (\sigma_ {r} - \nu \sigma_ {\theta})
\epsilon_ {\theta} = \frac {u}{r} = \frac {1}{E} (\sigma_ {\theta} - \nu \sigma_ {r}), \tag {3.44}
other
Prescribed axial displacement. φ = 0.25
| Section | Material | E | A | σγ | H' |
|---|---|---|---|---|---|
| Material 1 | 10 | 1000 | 1.0 | 10.0 | 0.0 |
| Material 2 | 10 | 1000 | 1.0 | 15.0 | -100.0 |
Fig. 3.13 Elasto-plastic example—Problem 3.8.
in which u is the radial displacement and E and \nu are respectively the elastic modulus and Poisson's ratio of the material. Also express the stresses \sigma_{r} and \sigma_{\theta} in terms of the nodal displacements \phi_{1} and \phi_{2} .
3.10 Use the stiffness matrix evaluated in Problem 3.9 to modify program ELPLAS to allow solution of one-dimensional axisymmetric problems by the initial stiffness method. Assume a Tresca yield criterion (discussed in Chapter 7) where yielding is assumed to begin when the maximum shearing stress reaches a critical value. For the present application this implies commencement of yielding when either \sigma_{r} or \sigma_{\theta} reaches the uniaxial yield stress, \sigma_{Y} .
3.11 Employ the program developed in Problem 3.10 to determine the elasto-plastic stress distribution in a thin disc, of thickness 1 mm, subjected to internal pressure loading. Take the internal and external
text_image
Axis of symmetry σθ dr σr dr σr + dσr dr dθ r f1 φ1 t φ2 f2 r1 r2
Fig. 3.14 Axisymmetric membrane element—Problem 3.9
radii of the disc as 5 cm and 10 cm respectively, the elastic modulus E=2\times10^{5} N/mm ^{2} , Poisson's ratio \nu=0\cdot3 and the uniaxial yield stress, \sigma_{Y}=300 N/mm ^{2} . Compare your solution with the theoretical expressions given in Ref. 8.
3.14 References
- RALSTON, A., A First Course in Numerical Analysis, McGraw-Hill, 1965.
- JENNINGS, A., Matrix Computation for Engineers and Scientists, John Wiley, 1977.
- HINTON, E. and OWEN, D. R. J., Finite Element Programming, Academic Press, London, 1977.
- SALVADORI, M. G. and BARON, M. L., Numerical Methods in Engineering, Prentice-Hall, New Jersey, 1964.
- HINTON, E. and OWEN, D. R. J., An Introduction to Finite Element Computations, Pineridge Press, Swansea, U.K., 1979.
- BARRER, R. M., Measurement of diffusion and thermal conductivity 'constants' in non-homogeneous media, and in media where these constants depend respectively on concentration or temperature, Proc. Phys. Soc. 58, 321–331 (1946).
- OWEN, D. R. J., PRAKASH, A. and ZIENKIEWICZ, O. C., Finite element analysis of non-linear composite materials by use of overlay systems, Computers and Structures, 4, 1251–1267 (1974).
- HILL, R., The Mathematical Theory of Plasticity, Oxford University Press, 1950.
Chapter 4 Viscoplastic problems in one dimension
4.1 Introduction
In this chapter the basic concepts of viscoplasticity are introduced by the consideration of one-dimensional situations. This topic is then studied further in Chapter 8 where the case of a general continuum is treated.
Viscoplastic theory allows the modelling of time rate effects in the plastic deformation process. Thus after initial yielding of the material the plastic flow, and the resulting stresses and strains, are time dependent. Such effects are always present to some degree in all materials but they may or may not be significant depending on the physical situation being considered.
The basic theory of viscoplasticity in one dimension is developed and a numerical solution process is then described. All the essential features of viscoplasticity can be demonstrated with reference to one-dimensional behaviour. Finally the solution process is coded in FORTRAN to form a working program and the basic characteristics of a viscoplastic material response are illustrated by the solution of numerical examples.
4.2 Basic theory
The concept of viscoplastic material behaviour is best introduced by means of the one-dimensional rheological model illustrated in Fig. 4.1. The friction slider component develops a stress \sigma_{p} , becoming active only if \sigma > Y , where \sigma is the total applied stress and Y is some limiting yield value. The excess stress \sigma_{d} = \sigma - \sigma_{p} is carried by the viscous dashpot. Instantaneous elastic response is, of course, provided by the linear spring. The presence of the dashpot allows the stress level to instantaneously exceed the value predicted by plasticity theory, the solution tending to this equilibrium level as steady state conditions are achieved in the system.
The total strain in the model is given by the sum of the elastic and viscoplastic components as
\epsilon = \epsilon_ {e} + \epsilon_ {v p}. \tag {4.1}
The stress in the linear spring is equal to the total applied stress and is
text_image
σd = σ - σp. σp εvp εe σp Inactive if σp < Y σp
Fig. 4.1 Basic one-dimensional elastic-viscoplastic model.
related to the elastic strain by
\sigma_ {e} = \sigma = E \epsilon_ {e}, \tag {4.2}
where E is the elastic modulus of the linear spring.
The stress level in the friction slider depends on whether or not the threshold or yield stress, Y, has been reached. The onset of viscoplastic deformation is governed by a uniaxial yield stress \sigma_{Y} . The stress level for continuing viscoplastic flow depends on the strain-hardening characteristics of the material. Restricting discussion to a linear strain-hardening response as discussed in Section 2.5, the stress level for viscoplastic yielding at any stage is given by
Y = \sigma_ {Y} + H ^ {\prime} \epsilon_ {v p}, \tag {4.3}
in which H' is the slope of the strain hardening portion of the stress-strain curve after removal of the elastic strain component and \epsilon_{vp} is the current viscoplastic strain. Thus the stress in the friction slider is
\begin{array}{r l} \sigma_ {p} = \sigma & \text { if } \left\{ \begin{array}{l} \sigma_ {p} < Y \\ \sigma_ {p} \geqslant Y. \end{array} \right. \\ = Y \end{array} \tag {4.4}
The stress in the viscous dashpot, \sigma_{d} , is related to the viscoplastic strain by
\sigma_ {d} = \mu \frac {d \epsilon_ {v p}}{d t}, \tag {4.5}
where \mu is a viscosity coefficient and t denotes the time. We note that
\sigma = \sigma_ {d} + \sigma_ {p}. \tag {4.6}
Before the onset of viscoplastic yielding \epsilon_{vp}=0 , giving \sigma_{d}=0 from (4.5) and consequently \sigma_{p}=\sigma . It now remains to establish the constitutive relationship for the model under both elastic and elasto-viscoplastic conditions.
Before viscoplastic yielding, \epsilon_{vp} = 0 and from (4.1) and (4.2) we have the elastic stress–strain relation to be
\sigma = E \epsilon . \tag {4.7}
Substituting from (4.4) and (4.5) in (4.6) gives
\sigma_ {Y} + H ^ {\prime} \epsilon_ {v p} + \mu \frac {d \epsilon_ {v p}}{d t} = \sigma . \tag {4.8}
Substituting for \epsilon_{vp} from (4.1) and using (4.2) results in
H ^ {\prime} E \epsilon + \mu E \frac {d \epsilon}{d t} = H ^ {\prime} \sigma + E (\sigma - \sigma_ {Y}) + \mu \frac {d \sigma}{d t}, \tag {4.9}
which is a first order ordinary differential equation defining the time-dependent relationship between stress and strain under viscoplastic conditions. At this stage we introduce a fluidity parameter, \gamma , such that
\gamma = \frac {1}{\mu}. \tag {4.10}
Substituting in (4.9) and rearranging
\dot {\epsilon} = \frac {\dot {\sigma}}{E} + \gamma [ \sigma - (\sigma_ {Y} + H ^ {\prime} \epsilon_ {v p}) ], \tag {4.11}
in which ( \cdot ) denotes derivative with respect to time, t. Or
\dot {\epsilon} = \dot {\epsilon} _ {e} + \dot {\epsilon} _ {v p}, \tag {4.12}
where
\dot {\epsilon} _ {e} = \frac {\dot {\sigma}}{E}, \tag {4.13}
and
\dot {\epsilon} _ {v p} = \gamma [ \sigma - (\sigma_ {Y} + H ^ {\prime} \epsilon_ {v p}) ]. \tag {4.14}
Expression (4.14) defines the viscoplastic strain rate in terms of the portion of stress in excess of the steady state yield value.
It is instructive to consider the closed form solution to (4.9). Consider the case when a constant applied stress \sigma = \sigma_{A} is applied to the model. Then (4.9) reduces, (using (4.10)), to
\gamma H ^ {\prime} \epsilon + \frac {d \epsilon}{d t} = \frac {\gamma H ^ {\prime}}{E} \sigma_ {A} + \gamma (\sigma_ {A} - \sigma_ {Y}). \tag {4.15}
The solution to this first-order ordinary differential equation is elementary and is
\epsilon = \frac {\sigma_ {A}}{E} + \frac {(\sigma_ {A} - \sigma_ {Y})}{H ^ {\prime}} [ 1 - \mathrm{e} ^ {- H ^ {\prime} \gamma t} ], \tag {4.16}
line
| Time, t | Strain, ε | | ------- | --------- | | 0 | 0 | | t | σ_A/E | | t | σ_A - σ_Y/H' |(a)
line
| Time, t | Strain, ε |
|---|---|
| 0 | 0 |
| 1 | 1 |
| 2 | 2 |
| 3 | 3 |
| 4 | 4 |
| 5 | 5 |
| 6 | 6 |
| 7 | 7 |
| 8 | 8 |
| 9 | 9 |
| 10 | 10 |
(b)
Fig. 4.2 Strain response with time for the model of Fig. 4.1 due to a constant applied load. (a) Linear strain hardening material. (b) Perfectly plastic material.
provided that H' is nonzero. The form of the response is shown in Fig. 4.2(a). Following an initial elastic response, the strain in the model attains the steady state value indicated in an exponential fashion.
The case of a perfectly viscoplastic material in which H' = 0 , can be obtained by taking the limit as H' tends to zero in (4.16) and applying L'Hopital's rule. This results in
\epsilon = \frac {\sigma_ {A}}{E} + \left(\sigma_ {A} - \sigma_ {Y}\right) \gamma t. \tag {4.17}
This response is shown in Fig. 4.2(b). In this case it is seen that a steady state condition is not achieved and that viscoplastic deformation continues indefinitely at a constant strain rate. The different behaviour shown in Figs. 4.2(a) and 4.2(b) arises from the fact that for a strain hardening material the viscoplastic yield stress increases according to (4.3) until it reaches the applied stress level \sigma_{A} at which stage the viscoplastic strain rate becomes zero. On the other hand, for a perfectly viscoplastic material there is always a stress imbalance of \sigma_{A}-\sigma_{Y} in the system which does not reduce and consequently steady state conditions cannot be achieved.
We note that in (4.16) and (4.17) that the time t only enters the expressions through the term \gamma t . Therefore the solution for a material with a different fluidity parameter \gamma can be obtained by a simple adjustment of the time scale.
4.3 Numerical solution process
Viscoplasticity is a transient phenomenon and therefore the essential objective of a numerical solution process is to determine the displacement, strains and stresses throughout the time interval of interest. Consequently some time stepping or time marching scheme must be introduced in order to allow the solution to be advanced from a time t_{n} to time t_{n+1} = t_{n} + \Delta t_{n} , where subscripts n and n+1 denote successive times and \Delta t_{n} the interval between. The simplest method of incrementing quantities over a time interval is afforded by Euler's rule. In this the mean rate of change over the interval is taken as the value at the beginning of the interval and thus the predicted value of some quantity X at time t_{n+1} is extrapolated from the value at time t_{n} to be
X ^ {n + 1} = X ^ {n} + (\dot {X}) ^ {n} \Delta t _ {n}. \tag {4.18}
This scheme becomes unstable for time steps exceeding a critical value and estimation of the limiting step length is discussed in Section 4.4. The Euler method, however, remains attractive due to its simplicity.
With the viscoplastic strain rate defined by (4.14) we can define the strain increment \Delta\epsilon_{vp}^{n} occurring in a time interval \Delta t_{n}=t_{n+1}-t_{n} , using (4.18), as
\Delta \epsilon_ {v p} ^ {n} = \dot {\epsilon} _ {v p} ^ {n} \Delta t _ {n}. \tag {4.19}
We note that the time step length can, in general, be different for each time interval.
text_image
Δφ₁ f₁ Δφ₂ f₂ 1 2 L
Fig. 4.3 One-dimensional two-noded element with linear displacement variation.
With reference to Fig. 4.3, consider the behaviour of a linear displacement element, which is of length L and has a cross-sectional area, A. The change of length in this element associated with strain increment (4.19) is
\Delta \phi^ {n} = \Delta \epsilon_ {v p} ^ {n} L, \tag {4.20}
or adding the displacement change due to a change in applied loading \Delta f^{n} occurring between times t_{n} and t_{n+1} we obtain the total change in element length to be
\Delta \phi^ {n} = \Delta \epsilon_ {v p} ^ {n} L + \frac {L}{A E} \Delta f ^ {n}. \tag {4.21}
This can be rewritten in matrix form, in terms of the nodal displacements and forces as
\Delta \varphi^ {n} = [ K ] ^ {- 1} \Delta V ^ {n}, \tag {4.22}
where
\Delta \varphi^ {n} = \left[ \begin{array}{l} \Delta \phi_ {1} ^ {n} \\ \Delta \phi_ {2} ^ {n} \end{array} \right], \tag {4.23}
\Delta V ^ {n} = A E \dot {\epsilon} _ {v p} ^ {n} \Delta t _ {n} \left[ \begin{array}{c} 1 \\ - 1 \end{array} \right] + \Delta f ^ {n}, \tag {4.24}
and
\boldsymbol {K} ^ {(e)} = \frac {E A}{L} \left[ \begin{array}{c c} 1 & - 1 \\ - 1 & 1 \end{array} \right]. \tag {4.25}
In the above, \Delta V^{n} are termed the pseudo forces and \Delta\varphi^{n} and \Delta f^{n} are respectively the incremental changes in the nodal displacements and applied forces for the element.