김경종
28 KiB
and J = \det \mathbf{X} = \det \mathbf{X}^E = {}^0\rho /{}^\prime \rho . If we assume, for the moment, that \mathbf{X}^p is known, we have
\mathbf {X} ^ {E} = \mathbf {X} \left(\mathbf {X} ^ {P}\right) ^ {- 1} \tag {6.287}
Also, since the velocity gradient \mathbf{L} = \dot{\mathbf{X}}\mathbf{X}^{-1} [see (6.40)], we may write
\mathbf {L} = \mathbf {L} ^ {E} + \mathbf {L} ^ {P} \tag {6.288}
where by substituting from (6.285), the elastic and plastic parts of the velocity gradient are, respectively,
\mathbf {L} ^ {E} = \dot {\mathbf {X}} ^ {E} (\mathbf {X} ^ {E}) ^ {- 1}; \quad \mathbf {L} ^ {P} = \mathbf {X} ^ {E} \dot {\mathbf {X}} ^ {P} (\mathbf {X} ^ {P}) ^ {- 1} (\mathbf {X} ^ {E}) ^ {- 1} \tag {6.289}
The variables that characterize the large strain elastoplastic response are therefore X, X^{P} , \tau , and \sigma_{y} , where \tau are the Cauchy stresses and \sigma_{y} is the current yield stress (including the effect of hardening). The evaluation of \tau must be based on evolution equations for X^{P} and \sigma_{y} (where in comparison to small strain plasticity X^{P} is used instead of e^{P} ). Since the Cauchy stresses are calculated throughout the solution, the constitutive description is used efficiently in an updated Lagrangian formulation.
Since \det X^{P} = 1 , we have J = \det X^{E} > 0 and we can calculate the polar decomposition
\mathbf {X} ^ {E} = \mathbf {R} ^ {E} \mathbf {U} ^ {E} \tag {6.290}
The logarithmic strain is used in the large strain one-dimensional response characterization (see Fig. 6.16), hence it is natural to use the elastic Hencky or logarithmic strain,
\mathbf {E} ^ {E} = \ln \mathbf {U} ^ {E} \tag {6.291}
in the multidimensional characterization of the Cauchy stress. We note that the evaluation of \mathbf{E}^E requires a spectral decomposition of \mathbf{U}^E [see (6.32) and (6.52)].
Since E^{E} is associated with the (conceptual) intermediate configuration \tau , we next define a stress measure, \overline{\tau} , corresponding to that configuration,
\overline {{{\boldsymbol {\tau}}}} = J \left(\mathbf {R} ^ {E}\right) ^ {T} \boldsymbol {\tau} \mathbf {R} ^ {E} \tag {6.292}
With these stress and strain measures we have the elastic work conjugacy (see S. N. Atluri [A] and Exercise 6.86),
\overline {{{\boldsymbol {\tau}}}} \cdot \dot {\mathbf {E}} ^ {E} = J \boldsymbol {\tau} \cdot \mathbf {D} ^ {E} \tag {6.293}
where D^{E} = \text{sym} (L^{E}) . The appropriate plastic velocity gradient to use is then ^{16}
\begin{array}{l} \overline {{{\mathbf {L}}}} ^ {P} = (\mathbf {X} ^ {E}) ^ {- 1} \mathbf {L} ^ {P} \mathbf {X} ^ {E} \\ = \dot {\mathbf {X}} ^ {p} \left(\mathbf {X} ^ {p}\right) ^ {- 1} \tag {6.294} \\ \end{array}
^{16} We can prove that with the definitions of \bar{\tau} in (6.292) and \bar{L}^{p} in (6.294), we have
J _ {\boldsymbol {\tau}} \cdot \mathbf {D} = \overline {{{\boldsymbol {\tau}}}} \cdot \dot {\mathbf {E}} ^ {E} + \overline {{{\boldsymbol {\tau}}}} \cdot \overline {{{\mathbf {D}}}} ^ {P}
where D is the (total) velocity strain [see (6.41)] and \overline{D}^{p} is given in (6.298) (see Exercise 6.87).
The stress \bar{\tau} is given by the usual stress-strain relationship of isotropic elasticity. Let \bar{S} be the deviatoric stress components and \bar{\sigma}_{m} be the mean stress, then
\overline {{{\mathbf {S}}}} = 2 \mu \mathbf {E} ^ {E ^ {\prime}}; \quad \overline {{{\sigma}}} _ {m} = 3 \kappa E _ {m} ^ {E} \tag {6.295}
where the elastic deviatoric strain components are given in E^{E'} and E_{m}^{E} is the elastic mean strain component. This choice of (hyperelastic) stress-strain law uses the total elastic strain and has the advantage of providing an excellent description of the stress even when the elastic strains are of moderate size (see L. Anand [A]).
The yield condition is as in (6.222), but using the deviatoric Cauchy stresses S, hence we have ^{17}
\overline {{{\sigma}}} = \sqrt {\frac {3}{2} \mathbf {S} \cdot \mathbf {S}} \tag {6.296}
= J ^ {- 1} \sqrt {\frac {3}{2} \bar {\mathbf {S}} \cdot \bar {\mathbf {S}}}
where \bar{S} is the deviatoric stress corresponding to the relaxed configuration, \bar{S} = J(\mathbf{R}^{E})^{T}\mathbf{S}\mathbf{R}^{E} . We also note that the unit normal to the yield surface in the relaxed hypothetical configuration is
\overline {{\mathbf {n}}} = \sqrt {\frac {3}{2}} \frac {\overline {{\mathbf {S}}}}{J \overline {{\sigma}}} \tag {6.297}
To obtain the evolution equation of the plastic deformation gradient, we use (6.293) and the plastic velocity strain tensor,
\widetilde {\mathbf {D}} ^ {P} = \operatorname{sym} \left(\widetilde {\mathbf {L}} ^ {P}\right) \tag {6.298}
This strain tensor is obtained from the flow rule, in analogy to (6.225), from ^{18}
\overline {{{\mathbf {D}}}} ^ {P} = \sqrt {\frac {3}{2}} \dot {\overline {{{e}}}} ^ {P} \overline {{{\mathbf {n}}}} \tag {6.299}
where \dot{\overline{e}}^P = \sqrt{\frac{2}{3}\overline{\mathbf{D}}^P\cdot\overline{\mathbf{D}}^P} (6.300)
Substituting from (6.297) into (6.299) corresponds to (6.225) in the small strain case. Of course, the relation \dot{\overline{\sigma}} = f(\overline{e}^P) is obtained from the uniaxial stress-strain relationship in Fig. 6.16. Consistent with the other assumptions used, the modified plastic spin tensor \overline{\mathbf{W}}^P = skw(\overline{\mathbf{L}}^P) is assumed to be zero.
Hence, we notice that the basic equations used for the solution of infinitesimal strain problems have been generalized to large strains by use of the Cauchy stress–logarithmic strain relationship in Fig. 6.16, the corresponding elastic stress-strain law for the Cauchy stress and logarithmic (Hencky) strain, and a plastic deformation gradient that represents the inelastic response and conceptually a deformation from which the elastic response is measured. Since the large strain formulation is a direct extension of the formulation used for small strains (see also Exercise 6.84), the computational procedures discussed earlier are also directly applicable and Table 6.10 summarizes the sequence of solution steps (see also A. L. Eterovic and K. J. Bathe [A] and F. J. Montáns and K. J. Bathe [A]).
In the preceding discussion we assumed a general three-dimensional response, but the equations can also be used directly for two-dimensional analyses, that is, plane stress, plane strain, and axisymmetric situations, by imposing the fundamental conditions of the specific
TABLE 6.10 Large strain elastoplastic updated Lagrangian Hencky formulation
The trial elastic state:
Obtain the trial elastic deformation gradient
\mathbf {X} _ {*} ^ {E} = ^ {I + \Delta} _ {0} \mathbf {X} (_ {0} ^ {I} \mathbf {X} ^ {P}) ^ {- 1}
Perform the polar decomposition
\mathbf {X} _ {*} ^ {E} = \mathbf {R} _ {*} ^ {E} \mathbf {U} _ {*} ^ {E}
Obtain the trial elastic strain tensor
\mathbf {E} _ {*} ^ {E} = \ln \mathbf {U} _ {*} ^ {E}
Obtain the trial elastic stress tensor \overline{\tau}_{*} using the equation for the mean stress and the deviatoric stress
\operatorname{tr} \left(\overline {{{\tau}}} _ {*}\right) = 3 \kappa \operatorname{tr} \left(\mathbf {E} _ {*} ^ {E}\right)
\bar {\mathbf {S}} _ {*} = 2 \mu \mathbf {E} _ {*} ^ {E}
Obtain the trial equivalent tensile stress
\overline {{{\sigma}}} _ {*} = J ^ {- 1} \sqrt {\frac {3}{2} \overline {{{\bf S}}} _ {*} \cdot \overline {{{\bf S}}} _ {*}}
Check whether the solution step corresponds to elastic conditions:
If \overline{\sigma}_{*} < ^{\prime}\sigma_{y} , then the solution step was elastic and we conclude the solution step by setting
{ } ^ { t + \Delta t } \overline { { { \sigma } } } = \overline { { { \sigma } } } _ { * } ; \qquad { } ^ { t + \Delta t } \sigma _ { y } = { } ^ { t } \sigma _ { y } ; \qquad { } ^ { t + \Delta t } \mathbf { E } ^ { E } = \mathbf { E } _ { * } ^ { E } ; \qquad { } ^ { t + \Delta t } \overline { { { \tau } } } = \overline { { { \tau } } } _ { * } ; \qquad { } ^ { t + \Delta t } \overline { { { \tau } } } _ { * } ( \mathbf { R } _ { * } ^ { E } ) ^ { T }
Otherwise we continue as follows.
Plastic solution step:
Use the effective-stress-function algorithm to calculate t^{+\Delta t}\overline{\sigma} and t^{+\Delta t}\overline{e}^P (see 6.233)
Calculate the stress deviator of t+\Delta t\overline{\tau} ,
\lambda = \frac {3}{2} \frac {^ {t + \Delta t} \overline {{e}} ^ {P} - ^ {t} \overline {{e}} ^ {P}}{^ {t + \Delta t} J ^ {t + \Delta t} \overline {{\sigma}}}
^ {t + \Delta t} \overline {{\mathbf {S}}} = \frac {\overline {{\mathbf {S}}} _ {*}}{1 + 2 \mu \lambda}
Calculate ^{t+\Delta t}\overline{\tau} ,
^ {t + \Delta t} \overline {{{\boldsymbol {\tau}}}} = ^ {t + \Delta t} \overline {{{\boldsymbol {S}}}} + \frac {1}{3} \operatorname{tr} (\overline {{{\boldsymbol {\tau}}}} _ {*}) \mathbf {I}
Calculate the Cauchy stresses,
{ } ^ { t + \Delta t } \mathbf { \tau } = ( { } ^ { t + \Delta t } J ) ^ { - 1 } \mathbf { R } _ { * } ^ { E } { } ^ { t + \Delta t } \overline { { { \mathbf { \tau } } } } ( \mathbf { R } _ { * } ^ { E } ) ^ { T }
Update the plastic deformation gradient by integration of (6.294),
{ } _ { 0 } ^ { t + \Delta t } \mathbf { X } ^ { P } = \exp \left( \lambda ^ { t + \Delta t } \overline { { \mathbf { S } } } \right) { } _ { 0 } ^ { t } \mathbf { X } ^ { P }
two-dimensional response being considered (see Fig. 4.5). In beam and shell analyses, these equations can also be employed, but the kinematic relations in Sections 6.5 must then include the effects of the change in thickness of the elements.
The large strain formulation has been presented above for elastoplastic response, but we note that the same approach can also be employed for creep and viscoplastic response—because of the analogies mentioned in Section 6.6.3—if the appropriate substitutions of variables and experimentally obtained formulas are made (see Exercise 6.90).
Finally, we should note that because of the incompressibility constraint on the inelastic strains (6.286), it is important to employ the displacement/pressure formulations discussed in Section 6.4 when analyzing two-dimensional plane strain, axisymmetric, or fully
three-dimensional response situations. This observation holds in elastoplasticity for small strain but particularly large strain conditions. Even in small strain analysis, the plastic strains are usually much larger than the elastic strains, and this is certainly the case in large strain conditions.
6.6.5 Exercises
6.61. Consider the eight-node brick element in Fig. 6.8. Plot the force-displacement responses assuming plane stress conditions in the y and z directions.
6.62. Consider the eight-node brick element in Fig. 6.8. Show explicitly that using (6.188) to transform \tilde{E} (\equiv_{i}^{t}C_{ijrs}) in the total Lagrangian formulation (i) for the figure, the force-displacement response is as calculated in (ii) in the figure. Also, show that using (6.187) to transform \tilde{E} (\equiv_{0}^{t}C_{ijrs}) in the updated Lagrangian formulation (ii) in the figure, the force-displacement response is as calculated in (i) in the figure.
6.63. A four-node element spins without deformation about its center at a constant angular velocity \omega , as shown. Use that the Jaumann stress rate is zero to calculate the Cauchy stresses (corresponding to the axes x_{1} , x_{2} ) at any time t.
text_image
x2 100 x1 200 ω 100 200 Σij = 0 = tij + tip twpj + tjp twpi
6.64. Consider the Mooney-Rivlin material description in (6.199). Show that this formula results in a pressure, the value of which depends on \delta I_{1} and \delta I_{2} . Then show that in the description (6.203) only the last term with the bulk modulus results in a pressure.
6.65. Consider the three-term Ogden material description in (6.205). Show that this formula results in a pressure as a function of the stretches, whereas in expression (6.207) the terms under the summation sign do not affect the pressure.
6.66. Show that for the two-term Mooney-Rivlin model (6.199), in small strains, the Young's modulus is given by 6(C_1 + C_2) and the shear modulus is given by 2(C_1 + C_2) . Also, show that these moduli are given in the Ogden model (6.205) by the values \frac{3}{2} \Sigma_{n=1}^{3} \alpha_n \mu_n and \frac{1}{2} \Sigma_{n=1}^{3} \alpha_n \mu_n , respectively.
6.67. Consider the four-node element in plane strain conditions shown. Calculate the force-displacement response for case 1 and case 2. Assume the bulk modulus \kappa to be very large.
text_image
2 mm x2 x1 2 mm P/2 P/2
Plane strain condition
Case 1:
Mooney-Rivlin material
C _ {1} = 0. 3 \mathrm{MPa}
C _ {2} = 0. 2 \mathrm{MPa}
Case 2:
Ogden material
\mu_ {1} = 1 \mathrm{MPa} \quad \alpha_ {1} = 2
\mu_ {2} = - 0. 5 \mathrm{MPa} \alpha_ {2} = - 1
\mu_ {3} = 0. 1 \mathrm{MPa} \quad \alpha_ {3} = 4
6.68. Consider the deformation of the four-node element shown. Calculate the force-displacement response. Assume the bulk modulus \kappa to be very large.
text_image
2 mm x2 x1 2 mm
Plane strain condition
Mooney-Rivlin material
C _ {1} = 0. 3 \mathrm{MPa}
C _ {2} = 0. 2 \mathrm{MPa}
6.69. Assume that instead of (6.199) the following higher-order Mooney-Rivlin material model is used:
_ {0} \tilde {W} = C _ {1} (_ {0} I _ {1} - 3) + C _ {2} (_ {0} I _ {2} - 3) + C _ {3} (_ {0} I _ {1} - 3) ^ {2} + C _ {4} (_ {0} I _ {1} - 3) (_ {0} I _ {2} - 3) + C _ {5} (_ {0} I _ {2} - 3) ^ {2}
with
C _ {1} = 7 5; \quad C _ {2} = 2 5; \quad C _ {3} = 1 0; \quad C _ {4} = 1 0; \quad C _ {5} = 1 0
Plot the force-displacement relationship for the one-dimensional bar problem in Fig. E6.25.
6.70. In plane stress solutions of incompressible response, we can use the pure displacement method of finite element analysis and adjust the element thickness to always fulfill the incompressibility constraint. Use this approach and derive for the Mooney-Rivlin law in (6.199) the stress-strain relationship \delta S_{ij} = \delta C_{ijrs} \delta \epsilon_{rs} and the tangent constitutive tensor _{0}C_{ijrs} .
6.71. Use a computer program to analyze the thick cylinder shown. The constitutive behavior is given by the Mooney-Rivlin law (6.203) with C_{1} = 0.6 MPa, C_{2} = 0.3 MPa, and \kappa = 2000 MPa.
The internal pressure increases uniformly until a maximum displacement of 10 mm is reached.
Use a sufficiently fine mesh to obtain an accurate solution. (Hint: The 9/3 element is a much more effective element than the 4/1 element.)
text_image
10 mm 10 mm 10 mm Pressure p Initial configuration
6.72. Use a computer program to solve for the response of the elastic circular thick plate due to a concentrated load at its center. Increase the load P until the deflection under the load is 2 cm. (Hint: Here axisymmetric elements of the u/p formulation are effective.)
text_image
P 2.0 cm Ogden material law Radius of plate R = 10 cm
\mu_ {1} = 0. 7 \mathrm{MPa}
\mu_ {2} = - 0. 3 \mathrm{MPa}
\mu_ {3} = 0. 0 1 \mathrm{MPa}
\alpha_ {1} = 1. 8; \alpha_ {2} = - 1. 6; \alpha_ {3} = 7. 5
\kappa = 1 0 0 0 \mathrm{MPa}
6.73. Show that in the uniaxial stress condition the effective stress \overline{\sigma} and effective plastic strain \overline{e}^{P} reduce to the uniaxial (nonzero) stress and corresponding plastic strain. Then assume that in a uniaxial stress experiment the stress varies as shown by the points 1 to 6. Plot the stress–effective plastic strain relation for the stress path.
line
| Point | e11 | σ11 | |-------|--------|--------| | 1 | 0.006 | 0.006 | | 2 | 0.018 | 0.018 | | 3 | 0.012 | 0.012 | | 4 | 0.006 | 0.006 | | 5 | 0.012 | 0.012 | | 6 | 0.018 | 0.018 |6.74. Prove that for a bilinear elastoplastic material described by the von Mises yield condition and flow rule, the effective stress function method gives the solution (6.236). State how the stress state at time t + \Delta t can therefore be computed if the complete state at time t is known and the strains at time t + \Delta t have been computed.
6.75. Show that the effective stress function in (6.233) can also be written in terms of the equivalent plastic strain as
\tilde {f} (e _ {*} ^ {P}) = 3 \mu (e _ {*} ^ {P} - ^ {t} e ^ {P}) + \sigma_ {y} (e _ {*} ^ {P}) - \overline {{\sigma}} _ {*}
where the solution \tilde{f}(e_{*}^{P}) = 0 corresponds to e_{*}^{P} = t + \Delta t e^{P} . Here \sigma_{y}(.) denotes the yield stress function and \overline{\sigma}_{*} is the effective stress corresponding to the elastic stress prediction [see (6.239)].
6.76. Consider isothermal von Mises plasticity with kinematic hardening in three-dimensional stress conditions in which case the yield condition is
{ } ^ { t + \Delta t } f _ { y } ^ { v M } = \frac { 1 } { 2 } { } ^ { t + \Delta t } \tilde { \mathbf { S } } \cdot { } ^ { t + \Delta t } \tilde { \mathbf { S } } - \frac { 1 } { 3 } ( \sigma _ { y v } ) ^ { 2 } = 0
where ^{t+\Delta t}\tilde{S} is the shifted deviatoric stress due to the back stress ^{t+\Delta t}\alpha
{ } ^ { t + \Delta t } \tilde { \mathbf { S } } = { } ^ { t + \Delta t } \mathbf { S } - { } ^ { t + \Delta t } \alpha
and \sigma_{yv} is the constant (virgin) yield stress. Assume small strain conditions and derive the effective-stress-function algorithm for this case.
6.77. Assume that the bilinear elastoplastic von Mises material in Fig. 6.10 with isotropic hardening is considered. Derive the tangent stress-strain matrix consistent with the effective-stress-function algorithm as indicated in (6.245) to (6.252).
6.78. Consider the creep law (6.263) with a_{0} = 6.4 \times 10^{-18} , a_{1} = 4.4 , a_{2} = 2.0 , and a_{3} = 0.0 . The stress states are \sigma = 100 \, MPa for 0 \leq t < 4 \, hr and \sigma = 200 \, MPa for t \geq 4 \, hr . For 0 \leq t \leq 10 \, hr draw the creep strain response as calculated by the time hardening and the strain hardening methods.
6.79. Derive the relation (6.268).
6.80. Derive the effective stress function given in (6.273) and show that near the solution this function has the curvature shown in Fig. 6.14.
6.81. A one-dimensional viscoplastic response is defined by \dot{\epsilon}^{VP} = \hat{\gamma}[\sigma - (\sigma_{yv} + E_{VP}\epsilon^{VP})] . The elastic strain response is given as usual by e^{E} = \sigma / E .
Calculate analytically the total strain response when the applied stress \sigma_{\mathrm{applied}} > \sigma_{\mathrm{yo}} . Consider the cases of E_{VP} > 0 and E_{VP} = 0 ; \hat{\gamma} = \text{constant} .
6.82. The constitutive behavior of a four-node plane stress element is given by the viscoplastic material model in (6.275) with the constants \beta = 10^{-4} (1/sec), N = 1. Also, E = 20,000 MPa, \nu = 0.3 , and the material effective stress–viscoplastic strain curve is given in the figure.
Assume a one-dimensional stress situation, that the stress at time 0 is zero, and that a constant load is applied as shown. Calculate the viscoplastic strain and displacement response of the element.
line
| Effective material stress | Value | | ------------------------- | ----- | | σyv | 20 | | EVP | 2 |
text_image
Load P = 40 Time
text_image
2 mm 2 cm Thickness 1 mm P P
6.83. Consider the general large strain elastoplastic formulation in Section 6.6.4. Derive from these general equations all equations for the one-dimensional response of the bar in Fig. 6.16.
6.84. Show explicitly for each equation in Table 6.10 that the large strain elastoplasticity formulation reduces to the formulation for materially-nonlinear-only analysis if the displacements and strains are small.
6.85. Consider the 8-node brick element in Fig. 6.8. Assume elastic conditions with E = 10^{7} N/cm ^{2} and \nu = 0.30 and plot the force-displacement response using the updated Lagrangian Hencky formulation of Section 6.6.4. Also plot this response for plane stress conditions in the y and z directions.
6.86. Consider an elastic material and show that (6.293) holds with the stress and strain quantities defined in Section 6.6.4. (Refer to the general continuum mechanics relations given in Section 6.2.2.)
6.87. Show that using the definitions in Section 6.6.4, we have J\pmb{\tau} \cdot \mathbf{D} = \overline{\pmb{\tau}} \cdot \dot{\mathbf{E}}^E + \overline{\pmb{\tau}} \cdot \overline{\mathbf{D}}^P , where \mathbf{D} = \frac{1}{2} (\mathbf{L} + \mathbf{L}^T) [see (6.41)].
6.88. Use the large strain elastoplasticity options in a computer program to calculate for a four-node plane stress element the stress response corresponding to the following strain path: In (a) the element is pulled out, in (b) the element is sheared over, in (c) the sheared element is pushed down, in (d) the element is brought to its original configuration. Explain whether your results make physical sense.
\text { Assume } \Delta = \frac {1}{1 0 0}, E = 2 0 0, 0 0 0, \nu = 0. 3, \sigma_ {y v} = 4 0 0 0.
text_image
1 1 Δ
(a)
natural_image
Simple line drawing of a structural frame with supports and a triangular load, no text or symbols present
(b)
natural_image
Simple line drawing of a rectangular frame with supports at both ends (no text or symbols)
(c)
text_image
Δ Δ Thickness = 1.0
(d)
6.89. Use a computer program to calculate the elastoplastic large strain response of the plane strain cantilever shown.
text_image
1 cm 20 cm Prescribed displacement δ
Let E = 200,000 MPa, \nu = 0.3 , E_T = 200 MPa, \sigma_y = 200 MPa and increase \delta to the value \frac{1}{2} L . Plot the force-displacement relationship and the stresses at \delta = L/1000 , \delta = L/10 , \delta = L/2 . Refine the finite element model to obtain a reasonably accurate stress prediction. (Hint: The 9/3 element performs well in large strain analysis.)
6.90. Assume that the large strain creep response of a material can be described by the theory in Section 6.6.4 with X^{P} replaced by the inelastic deformation gradient X^{IN} given by the creep deformations. Modify the entries in Table 6.10 to correspond to the creep solution.
Then use a computer program to calculate the large strain creep response of the thick cylinder shown. Obtain an accurate stress prediction for various values of internal pressure p.
text_image
p 1 cm 1 cm 1 cm E = 200,000 MPa v = 0.3 e^c = 10^-16 σ^4 t
text_image
Pressure p Time
6.7 CONTACT CONDITIONS
A particularly difficult nonlinear behavior to analyze is contact between two or more bodies. Contact problems range from frictionless contact in small displacements to contact with friction in general large strain inelastic conditions. Although the formulation of the contact conditions is the same in all these cases, the solution of the nonlinear problems can in some analyses be much more difficult than in other cases. The nonlinearity of the analysis problem is now decided not only by the geometric and material nonlinearities considered so far but also by the contact conditions.
The objective in this section is to briefly state the contact conditions in the context of a finite element analysis and present a general approach for solution.
6.7.1 Continuum Mechanics Equations
Let us consider N bodies that are in contact at time t . Let ^t S_c be the complete area of contact for each body L, L = 1, \ldots, N ; then the principle of virtual work for the N bodies at time t gives
\sum_ {L = 1} ^ {N} \left\{\int_ {t _ {V}} ^ {t} \tau_ {i j} \delta_ {t} e _ {i j} d ^ {t} V \right\} = \sum_ {L = 1} ^ {N} \left\{\int_ {t _ {V}} \delta u _ {i} ^ {t} f _ {i} ^ {B} d ^ {t} V + \int_ {t _ {S _ {f}}} \delta u _ {i} ^ {S} {} ^ {t} f _ {i} ^ {S} d ^ {t} S \right\} + \sum_ {L = 1} ^ {N} \int_ {t _ {S _ {c}}} \delta u _ {i} ^ {c} {} ^ {t} f _ {i} ^ {c} d ^ {t} S \tag {6.301}
where the part given in the braces corresponds to the usual terms (see Section 6.2.3) and the last summation sign gives the contribution of the contact forces. We note that the contact force effect is included as a contribution in the externally applied tractions. The components of the contact tractions are denoted as f_i^c and act over the areas S_c , and the components of the known externally applied tractions are denoted as f_i^S and act over the areas S_f . We might assume that the areas S_f are not part of the areas S_c , although such assumption is not necessary.
Figure 6.17 illustrates schematically the case of two bodies, which we now consider in more detail. The concepts given below can be directly generalized to multiple-body contact.
text_image
Body I S^IJ Time t Time 0 tS_c of bodies I and J is part of S^IJ and S^JJ Body J S^JJ Time t Time 0 x3 x2 x1
text_image
Body I tSc of body I S^IJ (contactor surface) tfIJ tfJI Body J tSc of body J S^JI (target surface) Bodies separated (conceptually) to show contact actions
Figure 6.17 Bodies in contact at time t
In Fig. 6.17 we denote the two bodies as body I and body J. Note that each body is supported such that without contact no rigid body motion is possible. Let f^{IJ} be the vector of contact surface tractions on body I due to contact with body J, then f^{IJ} = -f^{IJ} . Hence, the virtual work due to the contact tractions in (6.301) can be written as
\int_ {S ^ {I J}} \delta u _ {i} ^ {I} {} ^ {\prime} f _ {i} ^ {I J} d S ^ {I J} + \int_ {S ^ {J I}} \delta u _ {i} ^ {J} {} ^ {\prime} f _ {i} ^ {J I} d S ^ {J I} = \int_ {S ^ {I J}} \delta u _ {i} ^ {I J} {} ^ {\prime} f _ {i} ^ {I J} d S ^ {I J} \tag {6.302}
where \delta u_{i}^{I} and \delta u_{i}^{J} are the components of the virtual displacements on the contact surfaces of bodies I and J, respectively, and
\delta u _ {i} ^ {I J} = \delta u _ {i} ^ {I} - \delta u _ {i} ^ {J} \tag {6.303}