2.18.1 Design sensitivity analysis

Product: ABAQUS/Design  

ABAQUS/Design supports design sensitivity analysis (DSA) for static stress and frequency problems. 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. The DSA theory is presented from the perspective of computing the required derivatives analytically, first for static stress analysis and then for frequency analysis. At the end of each section an alternative numerical approach based on this theory is discussed.

DSA for static stress analysis

The equations pertaining to DSA can be derived based on a total displacement formulation or an incremental displacement formulation. The total displacement formulation is intended for history-independent problems, where the current state of the problem depends only on the total displacements. The incremental formulation is intended for history-dependent problems, where the current state of the problem depends on the state at the beginning of the increment and the incremental displacements.

Total displacement DSA formulation for nonlinear equilibrium problems

Let and be the numbers of design responses and design parameters, respectively. Let each response , , be a function of design parameters , and depend on them both explicitly and via the displacement field represented here by the nodal displacement vector (see the definition of finite element interpolation in Procedures: overview and basic equations, Section 2.1.1),

The dependence is only implicit; i.e., it is implied only by the design dependence of coefficients in the equilibrium equation system whose solution is .

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 as well as values of all responses. Sensitivity of a response with respect to design parameter is defined as

All but one quantity in the above equation can be determined explicitly given the equilibrium solution. The only unknown is ; to compute it, an additional system of equations has to be solved.

Rewrite Equation 2.1.1–2 in the form

where

All the quantities in the above equation are assumed to depend on design parameters explicitly or via displacement field . Differentiation of the above two equations with respect to design parameters leads to the following equation:

in which

is the tangent stiffness (Jacobian) matrix defined in Equation 2.1.1–4 and is an explicitly determinable quantity. Substituting Equation 2.18.1–3 into Equation 2.18.1–1, we obtain

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 , , and have to be computed in an element-by-element manner. 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.18.1–3 for each with respect to the unknown vectors of nodal displacement sensitivity . The displacement sensitivities are then substituted into Equation 2.18.1–1 to compute .

The coefficient matrix 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 in Equation 2.18.1–2 or to define design responses do not lend themselves to be expressed as functions of total displacement . Rather, at each time increment they are functions of certain state variables at the beginning of the increment (referred to as the time instant ) and of the incremental displacements, :

see, for example, Kleiber et al. (1997). The notation 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 . Some responses may also depend directly on the displacement , and the beginning-of-the-increment value of will, generally, also enter into the set .

In such a case Equation 2.18.1–4 takes the following form:

where

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 , (the pseudoload), and the derivatives with respect to displacements are assembled in the element loop. The set of design responses , , 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 :

and the solution vectors are substituted into Equation 2.18.1–5.

Computational approach

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 . 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 and numerically by finite differencing. For simplicity, assume that the finite difference technique is central difference such that for a given function , the derivative of with respect to is

where is the perturbation of .

For generality, consider the history-dependent case. To approximate the explicit design derivatives of , the incremental displacement is held constant while a positive perturbation is applied to each design parameter . In this way perturbed values of are obtained as

The change in the state corresponding to a perturbation in the design parameters is approximated by

The above process is repeated for a negative perturbation (), after which the results are differenced to arrive at the explicit design derivative .

Once the (incremental) displacement sensitivities are found, the response sensitivities can be obtained using

where

The process is repeated for a negative perturbation of , 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.

DSA for frequency analysis

This discussion will build upon the concepts and terminology described above, so it is recommended that the previous section be read first. The discussion below is divided into two sections depending on the characteristics of the eigenvalue problem: distinct eigenvalues and repeated eigenvalues.

Distinct eigenvalue case

The theory presented below assumes that all the eigenvalues are distinct (i.e., no repeated eigenvalues). If this is not the case, further manipulations are required to obtain the eigenvalue and eigenvector sensitivities corresponding to the repeated eigenvalues, and the following equations for the sensitivities will be incorrect. The repeated eigenvalue case is considered in the next section.

Performing a frequency analysis means solving the following eigenvalue problem (see Eigenvalue extraction, Section 2.5.1):

where represents the eigenvalues, represents the eigenvectors, and is the mass matrix. In addition, the eigenvectors are scaled such that either

or

for each mode. The default is the first scaling scheme. To obtain eigenvalue and eigenvector sensitivities, first differentiate Equation 2.18.1–6 with respect to design parameter to obtain the following equation:

where represents a particular mode number.

Pre-multiplying by , making use of Equation 2.18.1–6, and manipulating the result gives the eigenvalue sensitivities:

Except for the mass and stiffness derivatives, all quantities in this equation are known once the eigenvalue problem has been solved.

Repeated eigenvalue case

This section outlines the formulation used to obtain eigenvalue sensitivities for repeated eigenvalues. Further information can be found in the papers by Mills-Curran (1988) and Shaw (1991). When an eigenvalue repeats times, the eigenvectors associated with are linearly independent but not unique—any linear combination of these eigenvectors is also an eigenvector. Because of this non-uniqueness, the eigenvectors may not be continuous or differentiable in the design parameter. However, a set of eigenvectors that are continuous and differentiable can be obtained by an orthogonal transformation:

where , and is to be determined. Replacing the eigenvectors with the eigenvectors in Equation 2.18.1–9 and premultiplying by yields in matrix notation:

where , is a diagonal matrix,

and

Equation 2.18.1–12 is recognized as an reduced eigenvalue problem whose eigenvalues are the eigenvalue sensitivities associated with the repeated eigenvalue and whose eigenvectors are . sensitivities are obtained for the single eigenvalue . The physical interpertation of this is that a perturbation in the design parameter may cause the original repeated mode to branch into as many as distinct modes (imagine a beam with a circular cross-section perturbed into an elliptical section; the repeated modes associated with the original symmetry of the section now split into distinct modes associated with the minor and major axes of the ellipse).

Computational approach

A semi-analytic approach is used to compute the eigenvalue sensitivities. The basic idea of this approach, as outlined in the section on static DSA, is to compute some of the required intermediate derivatives using finite differencing. In the context of DSA for frequency procedures this means that the derivatives of the mass and stiffness matrices are computed using finite differencing.

Choosing the finite difference interval

ABAQUS uses a heuristic algorithm to automatically determine the perturbation sizes to be used in the differencing scheme. The objective of the algorithm is to select the perturbation sizes that lead to accurately computed derivatives. This is done on an element to element basis, so it is possible for a given design parameter that the perturbation size for one element will be different from that of another element. The selection of the perturbation sizes is based on the behavior of a scalar metric for each design parameter . This metric must be both convenient and representative of the terms that are numerically differenced. It is chosen as follows:

  • Static steps. For static steps the scalar metric is chosen as the norm of the element pseudoload:

    Since an accurate pseudoload is necessary to obtain an accurate sensitivity solution, the norm of the pseudoload is an appropriate choice for . This choice is also convenient since the pseudoload has to be computed in any case.

  • Frequency steps. For frequency steps the scalar metric is chosen as the projection of the matrix onto an eigenvector :

    The choice of depends on whether a given mode has a distinct eigenvalue or is associated with a repeated eigenvalue.

    If mode has a distinct eigenvalue, is taken as . Consequently, becomes simply the numerator of Equation 2.18.1–10. Therefore, is a direct measure of the magnitude of the eigenvalue sensitivity and is also convenient since this term already must be calculated to obtain the eigenvalue sensitivity.

    Unlike the distinct eigenvalue case, the sensitivities of a repeated eigenvalue cannot be treated independently. The sensitivities of a repeated eigenvalue are extracted simultaneously from the same reduced eigenvalue system, and this system is obtained by numerical differencing (recall Equation 2.18.1–12 and Equation 2.18.1–13). Consequently a single perturbation size (for each design parameter) must be used to obtain all sensitivities of a repeated eigenvalue. To calculate the single perturbation size, a single scalar metric is obtained by taking as the sum of the eigenvectors associated with the repeated eigenvalue. The calculation of is similar to the calculation of the matrix (the only term in the reduced eigenvalue system obtained by numerical differencing); therefore, this choice is both representative of the repeated eigenvalue case and involves differencing calculations that are already being done to obtain the reduced eigenvalue system.

The basic idea of the algorithm is compute the scalar metric for a range of perturbation sizes , , which vary by orders of magnitude. For each , a corresponding is computed. The relative change in between consecutive perturbation sizes, calculated as , is used to measure how close a perturbation size is to optimum. In this range of perturbation sizes the one that yields the smallest relative change, denoted as , is identified as the best perturbation size, . If is not less than a specified tolerance, the range of perturbation sizes is expanded (up to a certain limit), and the testing continues. It is important to realize that is not used directly in assessing the accuracy of the numerical differentiation but rather is intended as a means for determining the optimum perturbation size. Thus, a tighter tolerance on causes the algorithm to expend more effort in finding an optimum perturbation size but does not guarantee more accurate sensitivities.

Reference