김경종
28 KiB
line
| tan⁻¹(K∥/K⊥) | KII0 | MERR | MTS |
|---|---|---|---|
| 0 | 0 | 0 | 0 |
| 30 | 40 | 40 | 35 |
| 60 | 60 | 60 | 55 |
| 90 | 75 | 75 | 70 |
For the crack segment we also have the relation
G ^ {k} = \frac {1}{\bar {E}} (K _ {I} ^ {k ^ {2}} + {K _ {I I} ^ {k}} ^ {2}).
The maximum energy release rate criterion postulates that the parent crack initially propagates in the direction that maximizes G ^ { k } .
K _ {I I} = 0 \text { criterion }
This criterion simply postulates that a crack will initially propagate in the direction that makes K _ { I I } ^ { k } = 0 .
It can be seen from Figure 2.16.4-1 that the maximum energy release rate criterion and the K _ { I I } = 0 criterion predict nearly coincident crack propagation angles. By comparison, the maximum tangential stress criterion predicts smaller crack propagation angles.
2.17 Design sensitivity analysis
2.17.1 Design sensitivity analysis
ABAQUS/Design supports design sensitivity analysis (DSA) for nonperturbation, static stress problems that may include geometric nonlinearities and small-sliding, frictionless contact. DSA provides derivatives of certain response quantities with respect to specified input quantities. These derivatives are known as sensitivities. The responses available for DSA are a subset of the list of ABAQUS output variables and are known as design responses; the specified input quantities are
known as design parameters. Quantities that are functions of design parameters are referred to as being design dependent.
In the current capability only solid elements with elastic or hyperelastic properties may be made design dependent. This limits the permissible design parameters to those that affect nodal coordinates and elastic and/or hyperelastic material constants associated with solid elements. The permissible design responses are the set of "simple" output variables relevant to solid elements and static stress analysis. "Simple" responses are those that are not derived from other output variables. However, the discussion below is given in a general form without regard to the aforementioned limitations. The DSA theory is first presented from the perspective of computing the required derivatives analytically. In the final section an alternative numerical approach based on this theory is discussed.
Total displacement DSA formulation for nonlinear equilibrium problems
Let R and P be the numbers of design responses and design parameters, respectively. Let each response \phi _ { r } , r = 1 , \ldots , R . , be a function of design parameters h _ { p } , p = 1 , \ldots , P and depend on them both explicitly and via the displacement field represented here by the nodal displacement vector u ^ { N } (see the definition of finite element interpolation in ``Procedures: overview and basic equations, '' Section 2.1.1),
\phi_ {r} = \phi_ {r} \big (u ^ {N} (h _ {p}), h _ {p} \big).
The dependence u ^ { N } ( h _ { p } ) is only implicit; i.e., it is implied only by the design dependence of coefficients in the equilibrium equation system whose solution is u ^ { N } .
Assume that we have solved an equilibrium problem defined by Equation 2.1.1-2 at the end of an increment and that we have the converged solution u ^ { N } as well as values of all responses. Sensitivity of a response \phi _ { r } with respect to design parameter h _ { p } is defined as
Equation 2.17.1-1
\frac {d \phi_ {r}}{d h _ {p}} = \frac {\partial \phi_ {r}}{\partial h _ {p}} + \frac {\partial \phi_ {r}}{\partial u ^ {N}} \frac {d u ^ {N}}{d h _ {p}}.
All but one quantity in the above equation can be determined explicitly given the equilibrium solution. The only unknown is d u ^ { N } / d h _ { p } ; ; to compute it, an additional system of equations has to be solved.
Rewrite Equation 2.1.1-2 in the form
Equation 2.17.1-2
F ^ {N} (u ^ {M}) = 0,
where
F ^ {N} = - \int_ {V ^ {0}} \pmb {\beta} _ {N}: \pmb {\tau} ^ {c} d V ^ {0} + \int_ {S} \mathbf {N} _ {N} ^ {T} \cdot \mathbf {t} d S + \int_ {V} \mathbf {N} _ {N} ^ {T} \cdot \mathbf {f} d V.
All the quantities in the above equation are assumed to depend on design parameters h _ { p } explicitly or via displacement field u ^ { N } . Differentiation of the above two equations with respect to design parameters leads to the following equation:
Equation 2.17.1-3
K _ {M N} \frac {d u ^ {N}}{d h _ {p}} = \frac {\partial F ^ {M}}{\partial h _ {p}},
in which
K _ {M N} = \frac {\partial F ^ {M}}{\partial u ^ {N}}
is the tangent stiffness (Jacobian) matrix defined in Equation 2.1.1-4 and { \partial F ^ { M } } / { \partial h _ { p } } is an explicitly determinable quantity. Substituting Equation 2.17.1-3 into Equation 2.17.1-1, we obtain
Equation 2.17.1-4
\frac {d \phi_ {r}}{d h _ {p}} = \frac {\partial \phi_ {r}}{\partial h _ {p}} + \frac {\partial \phi_ {r}}{\partial u ^ {N}} K _ {N M} ^ {- 1} \frac {\partial F ^ {M}}{\partial h _ {p}},
which is the solution of the total displacement DSA problem.
The DSA algorithm used in ABAQUS is known as the direct differentiation method (DDM) and consists of the following operations. After the converged equilibrium solution is obtained, the three arrays \partial \phi _ { r } / \partial h _ { p } , \partial \phi _ { r } / \partial u ^ { M } , and { \partial F ^ { M } } / { \partial h _ { p } } have to be computed in an element-by-element manner. { \partial F ^ { M } } / { \partial h _ { p } } is often called the pseudoload since it becomes the right-hand side of the DSA problem. The final DSA solution is obtained by solving the system of Equation 2.17.1-3 for each p = 1 , \ldots , P with respect to the unknown vectors of nodal displacement sensitivity d u ^ { N } / d h _ { p } . The displacement sensitivities are then substituted into Equation 2.17.1-1 to compute d \phi _ { r } / d h _ { p } .
The coefficient matrix K _ { M N } used in the DSA computations is simply the last tangent stiffness matrix used in the equilibrium iterative algorithm. At the stage of the DSA computations this matrix is still available in the decomposed form and can be retrieved easily to perform the back substitutions for the DSA right-hand-side vectors. This makes the DSA module a very efficient add-on to the equilibrium analysis enabling sensitivity computations at a relatively low cost.
Incremental displacement DSA formulation for history-dependent equilibrium problems
The formulation of DSA presented above provides a brief introduction to the way DSA is implemented in ABAQUS; however, due to some simplifications, the discussion is not relevant to a large number of nonlinear mechanical problems, especially those involving history-dependent behavior of the structure modeled. The main difficulty in such problems is that many quantities necessary to compute the residual F ^ { N } in Equation 2.17.1-2 or to define design responses do not lend themselves to be expressed as functions of total displacement u ^ { N } . Rather, at each time increment, they are functions of certain state variables at the beginning of the increment (referred to as the time instant t) and of the
Procedures
incremental displacements, \Delta u ^ { N } :
F ^ {N} = F ^ {N} \left(\pmb {\alpha} ^ {t} (h _ {p}), \Delta u ^ {N} (h _ {p}), h _ {p}\right),
\phi_ {r} = \phi_ {r} \big (\pmb {\alpha} ^ {t} (h _ {p}), \Delta u ^ {N} (h _ {p}), h _ {p} \big);
see, for example, Kleiber et al. (1997). The notation \pmb { \alpha } ^ { t } stands for a set of state variables ® that may include tensors (stress, back stress, etc.) as well as scalar quantities (equivalent plastic strain, etc.) defined for a particular material point at time t. Some responses may also depend directly on the displacement u ^ { N } , and the beginning-of-the-increment value of u ^ { N } will, generally, also enter into the set \pmb { \alpha } ^ { t } .
In such a case Equation 2.17.1-4 takes the following form:
Equation 2.17.1-5
\frac {d \phi_ {r}}{d h _ {p}} = \frac {D \phi_ {r}}{D h _ {p}} + \frac {\partial \phi_ {r}}{\partial \Delta u ^ {M}} K _ {M N} ^ {- 1} \frac {D F ^ {N}}{D h _ {p}},
where
\frac {D (\cdot)}{D h _ {p}} \stackrel {\mathrm{def}} {=} \frac {\partial (\cdot)}{\partial h _ {p}} + \frac {\partial (\cdot)}{\partial \pmb {\alpha} ^ {t}} \frac {d \pmb {\alpha} ^ {t}}{d h _ {p}}
denotes the explicit design derivative of a quantity (¢).
The fundamental difference, from the point of view of the DSA solution algorithm, between the total and incremental approach is that in the latter case all state variables ® effectively become additional, or internal, design responses, whose sensitivities must be computed and updated at the end of each time increment to proceed with the DSA in the next increment. The number of such internal responses may be significant with obvious effects both on the computational time and memory requirement.
The DSA solution procedure is similar to that in the total displacement approach. After the equilibrium computations are complete, the arrays of explicit design derivatives D \phi _ { r } / D h _ { p } , D F ^ { M } / D h _ { p } (the pseudoload), and the derivatives with respect to displacements \partial \phi _ { r } / \partial \Delta u ^ { M } are assembled in the element loop. The set of design responses \phi _ { r } , r = 1 , \ldots , R , includes in this case all the scalars and tensor components of ®. In the direct differentiation method the following system of equations is solved for each design parameter h _ { p } :
K _ {M N} \frac {d \Delta u ^ {N}}{d h _ {p}} = \frac {D F ^ {M}}{D h _ {p}},
and the solution vectors are substituted into Equation 2.17.1-5.
Computational approach
Procedures
The derivatives required for DSA can be computed analytically or numerically. In the analytical approach the finite element equations are differentiated exactly, following the theory described in the previous sections. This approach is difficult to implement, but it is efficient and yields exact sensitivities. In the numerical approach some or all of the required derivatives are computed using the finite difference technique. The numerical approach can be further subdivided into the overall or global finite difference approach and the semi-analytic approach. In the global finite difference approach the response sensitivities with respect to a particular design parameter are obtained by perturbing that design parameter a number of times (depending on the finite difference technique) and performing an entire equilibrium analysis for each perturbation. The responses are retained for each analysis and then differenced to obtain the response sensitivities. This approach is computationally expensive since an entire equilibrium problem must be solved for each perturbation, but it is easily implemented. The semi-analytic approach is used in ABAQUS and can be viewed as a compromise between the analytic and global finite difference approaches. In the semi-analytic approach the DSA element vectors are obtained by differencing; but, like the analytic approach, the DSA solution is obtained by back-substitution against K _ { M N } . The advantage of the semi-analytic approach is that it is much easier to implement than the analytic approach and much more efficient than the global finite difference approach. The details of this method are described in the following paragraphs.
The objective of the semi-analytic approach is to compute the DSA vectors D F ^ { M } / D h _ { p } and d \phi _ { r } / d h _ { p } numerically by finite differencing. For simplicity, assume that the finite difference technique is central difference such that for a given function A ( x ) , the derivative of A with respect to x is
{\frac {d A}{d x}} = {\frac {A (x + \delta x) - A (x - \delta x)}{2 \delta x}},
where \delta x is the perturbation of x . .
For generality, consider the history-dependent case. To approximate the explicit design derivatives of F ^ { M } , the incremental displacement is held constant while a positive perturbation \delta h _ { p } is applied to each design parameter h _ { p } . In this way perturbed values of F ^ { M } are obtained as
F ^ {M} + \delta F ^ {M} = F ^ {M} (\pmb {\alpha} ^ {t} (h _ {p} + \delta h _ {p}), \Delta u (h _ {p}), h _ {p} + \delta h _ {p}).
The change in the state corresponding to a perturbation in the design parameters is approximated by
\pmb {\alpha} ^ {t} (h _ {p} + \delta h _ {p}) = \pmb {\alpha} ^ {t} (h _ {p}) + \frac {d \pmb {\alpha} ^ {t}}{d h _ {p}} \delta h _ {p}.
The above process is repeated for a negative perturbation ( - \delta h _ { p } ) , after which the results are differenced to arrive at the explicit design derivative D F ^ { M } / D h _ { p } .
Once the (incremental) displacement sensitivities are found, the response sensitivities d \phi _ { r } / d h _ { p } can be obtained using
\phi_ {r} + \delta \phi_ {r} = \phi_ {r} (\pmb {\alpha} ^ {t} (h _ {p} + \delta h _ {p}), \Delta u ^ {N} (h _ {p} + \delta h _ {p}), h _ {p} + \delta h _ {p}),
Procedures
where
\Delta u ^ {N} (h _ {p} + \delta h _ {p}) = \Delta u ^ {N} (h _ {p}) + \frac {d \Delta u ^ {N}}{d h _ {p}} \delta h _ {p}.
The process is repeated for a negative perturbation of h _ { p } , and the results are differenced.
The finite difference interval must be chosen carefully. If the interval is too small, round-off or cancellation errors occur due to loss of precision during the differencing operations. On the other hand, if the interval is too large, truncation errors may occur. Truncation errors arise from the fact that differencing formulas are based on truncated Taylor series expansions. ABAQUS will automatically choose a perturbation size that provides the best compromise between cancellation and truncation errors.
3. Elements
3.1 Overview
3.1.1 Element library: overview
The ABAQUS element library provides a complete geometric modeling capability. For this reason any combination of elements can be used to make up the model. Sometimes multi-point constraints are required for application of the necessary kinematic relations to form the model (for example, to model part of a shell surface with solid elements and part with shell elements or to model a pipe elbow with a mixture of beam and shell elements).
All elements use numerical integration to allow complete generality in material behavior. Shell and beam element properties can be defined as general section behaviors, or each cross-section of the element can be integrated numerically, so that nonlinear response can be tracked accurately when needed. A composite layered section can be specified, with different materials at different heights through the section. Some special elements (such as line springs) use an approximate analytical solution to model nonlinear behavior.
All of the elements in ABAQUS are formulated in a global Cartesian coordinate system except the axisymmetric elements, which are formulated in terms of r-z coordinates. In almost all elements, primary vector quantities (such as displacements u and rotations \phi ) are defined in terms of nodal values with scalar interpolation functions. For example, in elements with a two-dimensional topology the interpolation can be written as
\mathbf {u} (g, h) = N ^ {N} (g, h) \mathbf {u} ^ {N},
where the interpolation functions N ^ { N } ( g , h ) are written in terms of the parametric coordinates g and h. In most element types the same parametric interpolation is used for the coordinate vector x:
\mathbf {x} (g, h) = N ^ {N} (g, h) \mathbf {x} ^ {N}.
Such isoparametric elements are guaranteed to be able to represent all rigid body modes and homogeneous deformation modes exactly, a necessary condition for convergence to the exact solution as the mesh is refined.
All elements in ABAQUS are integrated numerically. Hence, the virtual work integral as described in ``Nonlinear solution methods in ABAQUS/Standard,'' Section 2.2.1, will be replaced by a summation:
\int_ {V} \boldsymbol {\sigma}: \delta \mathbf {D} d V \rightarrow \sum_ {i = 1} ^ {n} \boldsymbol {\sigma} _ {i}: \delta \mathbf {D} _ {i} V _ {i},
where n is the number of integration points in the element and V _ { i } is the volume associated with integration point i. ABAQUS will use either "full" or "reduced" integration. For full integration the
Elements
number of integration points is sufficient to integrate the virtual work expression exactly, at least for linear material behavior. All triangular and tetrahedral elements in ABAQUS use full integration. Reduced integration can be used for quadrilateral and hexahedral elements; in this procedure the number of integration points is sufficient to integrate exactly the contributions of the strain field that are one order less than the order of interpolation. The (incomplete) higher-order contributions to the strain field present in these elements will not be integrated.
The advantage of the reduced integration elements is that the strains and stresses are calculated at the locations that provide optimal accuracy, the so-called Barlow points ( Barlow, 1976). A second advantage is that the reduced number of integration points decreases CPU time and storage requirements. The disadvantage is that the reduced integration procedure can admit deformation modes that cause no straining at the integration points. These zero-energy modes make the element rank-deficient and cause a phenomenon called "hourglassing," where the zero energy mode starts propagating through the mesh, leading to inaccurate solutions. This problem is particularly severe in first-order quadrilaterals and hexahedra. To prevent these excessive deformations, an additional artificial stiffness is added to the element. In this so-called hourglass control procedure, a small artificial stiffness is associated with the zero-energy deformation modes. This procedure is used in many of the solid and shell elements in ABAQUS.
Most fully integrated solid elements are unsuitable for the analysis of (approximately) incompressible material behavior. The reason for this is that the material behavior forces the material to deform (approximately) without volume changes. Fully integrated solid element meshes, and in particular lower-order element meshes, do not allow such deformations (other than purely homogeneousdo deformation). For that reason ABAQUS uses "selectively reduced" integration in these elements: reduced integration is used for the volume strain and full integration for the deviatoric strains. As a consequence the lower-order elements give an acceptable performance for approximately incompressible behavior. For fully incompressible material behavior, another complication occurs: the bulk modulus and, hence, the stiffness matrix becomes infinitely large. For this case a mixed (hybrid) formulation is required, where the displacement field is augmented with a hydrostatic pressure field. In this formulation only the inverse of the bulk modulus appears, and, consequently, the contribution to the operator matrix vanishes. The hydrostatic pressure field plays the role of a Lagrange multiplier field enforcing the incompressibility constraints.
ABAQUS/Standard also provides elements for multifield problems. Examples are the pore pressure elements used for the analysis of porous solids with fluid diffusion, thermally coupled elements that couple heat transfer with stress analysis, and piezoelectric elements that couple electrical conduction with stress analysis. In these multifield elements the scalar variable (such as the temperature) is usually interpolated with different scalar functions as the displacement field; i.e.,
T (g, h) = M ^ {N} (g, h) T ^ {N},
where M ^ { N } ( g , h ) may differ from N ^ { N } ( g , h ) . The coupling of the fields will generally occur at the integration points; for example, in thermally coupled elements the coupling is due to temperature-dependent mechanical properties and heat generation is due to inelastic work. Finally, ABAQUS offers a complete set of diffusion elements to analyze conductive and convective heat
transfer. In these elements only temperatures appear as nodal degrees of freedom. The temperatures are interpolated with essentially the same interpolation function, M ^ { N } ( g , h ) , as used in the thermally coupled elements.
3.2 Continuum elements
3.2.1 Solid element overview
ABAQUS contains a library of solid elements for two-dimensional and three-dimensional applications. The two-dimensional elements allow modeling of plane and axisymmetric problems and include extensions to generalized plane strain (when the model exists between two planes that may move with respect to each other, providing thickness direction strain that may vary with position in the plane of the model but is constant with respect to thickness position). The material description of three-dimensional solid elements may include several layers of different materials, in different orientations, for the analysis of laminated composite solids. A set of nonlinear elements for asymmetric loading of axisymmetric models is also available, and linear infinite elements in two and three dimensions can be used to model unbounded domains.
The solid element library includes isoparametric elements: quadrilaterals in two dimensions and "bricks" (hexahedra) in three dimensions. These isoparametric elements are generally preferred for most cases because they are usually the more cost-effective of the elements that are provided in ABAQUS. They are offered with first- and second-order interpolation and are described in detail in ``Solid isoparametric quadrilaterals and hexahedra, '' Section 3.2.4. For practical reasons it is sometimes not possible to use isoparametric elements throughout a model; for example, some commercial mesh generators use automatic meshing techniques that rely on triangulation to fill arbitrarily shaped regions. Because of these needs ABAQUS includes triangular, tetrahedron, and wedge elements. For most cases it is recommended that these elements be only used to fill in awkward parts of the mesh and, in particular, that well-shaped isoparametric elements be used in any critical region (such as an area where the strain must be predicted accurately). The isoparametric elements can also be degenerated to make simpler shapes. Generally the elements written for those particular geometries are preferred to this method. The exception to this rule occurs in cases where singularities are to be modeled (such as in fracture mechanics applications), since the degenerate second-order isoparametric elements can provide a 1 / \sqrt { r } singularity through the use of the "quarter point" technique (placing the midside nodes 1/4 of the distance along the side from the node at the singularity instead of at the middle point of the side).
Solid elements are provided with first-order (linear) and second-order (quadratic) interpolation, and the user must decide which approach is more appropriate for the application. Some guidelines are as follows. Standard first-order elements are essentially constant strain elements: the isoparametric forms can provide more than constant strain response, but the higher-order content of the solutions they give is generally not accurate and, thus, of little value. The "incompatible mode" elements, described in ``Continuum elements with incompatible modes,'' Section 3.2.5, are from the user's perspective lower-order elements but have internal degrees of freedom that enable the element to represent almost all linear strain patterns. These elements can represent certain important linear strain fields exactly: the most important field is the one due to bending. The second-order elements are capable of representing
Elements
all possible linear strain fields. Thus, in the case of elliptic problems--problems for which the governing partial differential equations are elliptic in character, such as elasticity, heat conduction, acoustics, in which smoothness of the solution is assured--much higher solution accuracy per degree of freedom is usually available with the higher-order elements. Therefore, it is generally recommended that the highest-order elements available be used for such cases: in ABAQUS this means second-order elements. This observation logically leads to the use of the "hierarchical" finite element technique or "p"-method--refining the model by increasing the interpolation order in the elements in critical regions:"p"-method--refining the model by increasing the interpolation order in the elements in critical regions: this approach is as yet not available in ABAQUS.
A case where both incompatible mode elements and second-order elements can be used effectively isA case where both incompatible mode elements and second-order elements can be used effectively is the stress analysis of relatively thin members subjected to bending: such problems are often encountered in practical applications. In such cases the strain variation through the thickness must be at least linear, and constant strain (first-order) elements do a poor job of representing this variation. Fully integrated first-order isoparametric elements also suffer from "shear locking" in these geometries: they cannot provide the pure bending solution because they must shear at the numerical integration points to respond with an appropriate kinematic behavior corresponding to the bending. This shearing then locks the element--the response is far too stiff. For the isoparametric elements reduced integration provides a cure for these problems, but at the cost of allowing spurious singular modes ("hourglassing"). The use of second-order elements is a more reliable alternative, because the second-order interpolation naturally contains the linear strain field--one element through the thickness is enough to represent the behavior of a thin component subjected to bending loads quite accurately. Another alternative is formed by the incompatible mode elements: the linear strain field in these elements contains the modes required to solve the bending problem exactly if the elements are rectangular in shape. For a detailed discussion of the performance of ABAQUS continuum elements in bending problems, see ``Performance of continuum and shell elements for linear analysis of bending problems,'' Section 2.3.5 of the ABAQUS Benchmarks Manual. (It should be remembered, however, that ABAQUS offers shell and beam elements that are specifically written for thin geometries: the use of solid elements for such cases should only be considered when beam or shell elements are not practical.)
For all of these reasons the second-order elements are preferred in elliptic applications. The argument is readily extended to higher-order interpolation (cubic, quartic, etc), but the rapid increase in cost peris readily extended to higher-order interpolation (cubic, quartic, etc), but the rapid increase in cost per element for higher-order forms means that--even though the accuracy per degree of freedom is higher--the accuracy per computational cost may not be increasing. Practical experience suggests that--except in special cases--little is gained by going beyond the second-order elements, so ABAQUSthat--except in special cases--little is gained by going beyond the second-order elements, so ABAQUS does not offer any higher-order forms.
Many problems of practical interest are not elliptic: localizations arise in one form or another. Plasticity applications are an example--as the solution approaches the limit load, most plasticity models tend toward hyperbolic behavior. This allows discontinuities to occur in the solution--the slip line solutions of classical perfect plasticity theory are plots of the characteristic lines of velocity discontinuities in the hyperbolic equations of the problem. If the finite element solution is to exhibit accuracy, these discontinuities in the gradient field of the solution should be reasonably well modeled. With a fixed mesh that does not use special elements that admit discontinuities in their formulation, this suggests that the lowest-order elements--the first-order elements--are likely to be the most