Processing the sensitivity

The sensitivity is calculated in routines results_se.c, mafillsmmain_se.c and objectivemain_se.c.

In routine results_se.c $\frac{\partial F_{\text{int}}}{\partial s}$ is determined. For geometrically nonlinear calculations (parameter NLGEOM on the *STEP card) the unperturbed displacements leading to the internal force correspond to the displacements at the end of the previous static step, if any, augmented by the actual prescribed displacements. For linear geometric calculations the unperturbed displacements correspond to zero augmented by the actual prescribed displacements. Indeed, nonzero prescribed displacements lead to internal forces in linear calculations. Therefore, the term $\frac{\partial F}{\partial s}$ in Equation (435) can be replaced by $\frac{\partial F_{\text{ext}}}{\partial s}-\frac{\partial F_{\text{int}}}{\partial s}$ for linear calculations, noting that only nonzero initial displacement boundary conditions lead to internal forces (and not any previous displacements). The second reason why results_se.c has to be called for linear calculations too is that the material tangent $\frac{\partial S}{\partial E}$ at each integration point, which is needed in mafillsmmain_se.c for the set up of the stiffness matrix, is also determined in results_se.c.

Routine objectivemain_se.c calculates the derivative of the external forces and of the stiffness matrix (and similar matrices):

• For static nonlinear geometric calculations

  $$\frac{\partial F_{\text{ext}}}{\partial s}$$ (536)

  
  

  is calculated.

• For static linear calculations

  $$\frac{\partial F_{\text{ext}}}{\partial s}-\frac{\partial K}{\partial s} \cdot U$$ (537)

  
  

  is determined.

• For frequency and Green calculations

  $$-\frac{\partial (K- \sigma M)}{\partial s} \cdot U$$ (538)

  
  

  is calculated, where $\sigma$ is an appropriately defined scalar.

Out of computational efficiency the latter terms are calculated at the element level and assembled into a global matrix thereupon.

The last major routine, objectivemain_se.c assembles the previous information to obtain the final sensitivity. For the orientation as design variable these sensitivities are immediately stored in the .dat or the .frd file. The sensitivity for the geometry (normal directions of nodes on the external surface) as design variable, however, is kept for further postprocessing in sensitivity.c.

For the objective $G$ the total sensitivity $\frac{\partial G(s,U(s))}{\partial s}$ is written as $\frac{\partial G}{\partial s} + \frac{\partial G}{\partial U}$. The routine is split according to the objective function:

• The MASS objective function does not depend on the displacements, i.e. the deformation of a body does not change its mass. So only the first term in the above equation is needed. This term examines how the change of the design variables directly changes the mass. For the orientation as design variable the mass does not change at all. For coordinates as design variables, however, $\frac{\partial G}{\partial s}$ is appropriately calculated. This is done by determining $G$ for the unperturbed geometry and for the geometry in which one of the design variables (the geometric change in normal direction in a node on the external surface) is changed by a small amount (finite difference approximation). The routine in which this is done is objective_mass_dx.f.

  In general, the objective function does not have to apply to the total structure, e.g. one can define the mass of part of the structure as design variable. In that case all other elements are deactivated. This is done in routine actideacti.f. This applies to all objective functions, for which only part of the structure is included.

• For the SHAPE ENERGY as objective function a distinction has to be made whether the calculation is geometrically linear or nonlinear. For a linear geometric calculation Equation (436) reduces to:

  $$\frac{DG}{Ds}=\frac{\partial G}{\partial s} + U \left( \frac{\partial F}{\partial s}-\frac{\partial K}{\partial s} \cdot U \right).$$ (539)

  
  

  The first term on the right hand side is calculated in a similar way as for the MASS in routine objective_shapeener_dx.f. The term in brackets on the right hand side was already determined in results_se.f and mafillsmmain_se.f. Premultiplying it with the displacements from the previous static step and adding the first term yields the sensitivities (objective_shapeener_tot.f).

  For a nonlinear geometric calculation Equation (436) reduces to:

  $$\frac{DG}{Ds}=\frac{\partial G}{\partial s} + F_{\text{int}}^T K^... ...F_{\text{ext}}}{\partial s}-\frac{\partial F_{\text{int}}}{\partial s} \right).$$ (540)

  
  

  Now, $Y \equiv F_{\text{int}}^T K^{-1}$ is calculated by solving the set of equations $K Y = F_{\text{int}}$. The remaining operations are similar to the linear case.

• For the EIGENFREQUENCY as objective function the sensitivity reduces to:

  $$\frac{\partial \lambda_i}{\partial s} = U_i^T \cdot \left( \frac{... ...ial K}{\partial s} - \lambda_i \frac{\partial M}{\partial s} \right) \cdot U_i.$$ (541)

  
  

  The part starting with the brackets on the right hand side has been determined in mafillsmmain_se.f. Consequently,the sensitivity of the eigenfrequencies only requires the premultiplication with the eigenmodes. This is done in objective_freq.f and objective_freq_cs.f (cyclic symmetry).

  For the sensitivity of the eigenmodes (only calculated for the orientation as design variable) the relevant equation is Equation (440), which can also be written as:

  $$(K - \lambda_i M) \frac{\partial U_i}{\partial s} = F_i.$$ (542)

  
  

  Assuming the sensitivity to be a linear combination of the eigenmodes:

  $$\frac{\partial U_i}{\partial s} = \sum_j c_j U_j,$$ (543)

  
  

  leads to the following expressions for $c_j$:

  $$c_j=\frac{U_j^T F_i}{\lambda_j - \lambda_i}, i \ne j$$ (544)

  
  

  and $c_i=0$. The latter equation results from the differentiation of the mass normalization condition $U_i^T M U_i=1$. The determination of

  $$F_i= - \left( \frac{\partial K}{\partial s} - \lambda_i \frac{\partial M}{\partial s} - \frac{\partial \lambda_i}{\partial s} M \right) \cdot U_i$$ (545)

  
  

  is straightforward and is based on the expression calculated in mafillsmmain_se.f and the sensitivity of the eigenfrequencies.

• For the sensitivity of the Green functions (only calculated for the orientation as design variables) the relevant equation reads:

  $$(K - \omega_0^2 M) \frac{\partial U_i}{\partial s} = \frac{\partial K}{\partial s} \cdot U_i.$$ (546)

  
  

  and requires the solution of a system of equations for each design variable. The system matrix, however, does not change, so the LU-decomposition of the matrix has only to be done once.

  For the orientation as design variable the frequency sensitivities are stored in the .dat file, whereas the sensitivities of the eigenmodes and/or Green functions are stored in the .frd file (frd_sen.c, called from objectivemain_se.c). For the geometry as design variable only the frequency sensitivities are determined. They are not stored in objectivemain_se.c since they may need further postprocessing in sensitivity.c.

• The DISPLACEMENT objective function is differently defined for orientation design variables than for geometric design variables. The processing, however, is similar. The relevant equation is

  $$\frac{\partial U}{\partial s} = K^{-1} \left( \frac{\partial F_{\text{ext}}}{\partial s}-\frac{\partial F_{\text{int}}}{\partial s} \right),$$ (547)

  
  

  for geometrically nonlinear calculations and

  $$\frac{\partial U}{\partial s} = K^{-1} \left( \frac{\partial F}{\partial s}-\frac{\partial K}{\partial s} \cdot U \right)$$ (548)

  
  

  for geometrically linear calculations. In both cases the term in brackets on the right hand side (let us call it $F$) has been calculated before. Therefore, the complete right hand side is determined by solving $K Y = F$ for each design variable. Since the matrix of the system does not depend on the design variable, it has only to be LU-decomposed once.

  For orientation design variables the result is transferred from the degrees of freedom to the (node,direction) representation in resultsnoddir.f and stored in the .frd file in frd_sen.c. For geometrical design variables the result is processed in objective_disp_dx.f. This is due to the fact that the displacement geometric function for geometrical design variables is defined as the square root of the sum of the square of the displacements in all design nodes. After leaving objective_disp_dx the result is kept for further postprocessing.

• Results for the STRESS as objective are immediately based on the displacement sensitivities. Indeed, one can write:

  $$S=f(U(s),s),$$ (549)

  
  

  i.e. the stress S is a direct function of the design variables and an indirect function through the displacements. Indeed, the stress is the result of the ``multiplication'' of the material contants with the derivative of the displacements with respect to the geometry. Therefore, the indirect dependence through the displacements is clear. For geometrical design variables the direct dependence follow from ``WITH RESPECT TO THE GEOMETRY'' in the previous sentence, for the orientation design variables it follows from ``THE MATERIAL CONSTANTS''. Therefore, the sensitivity of the stress can be written as:

  $$\frac{\partial S}{\partial s} \approx \frac{S(U(s+\Delta s),s+\Delta s)-S(U(s),s)}{\Delta s},$$ (550)

  
  

  where $U(s+\Delta s)$ is approximated by

  $$U(s+\Delta s) \approx U(s) + \frac{\partial{U}}{\partial{s}} \Delta s.$$ (551)

  
  

  So the stress sensitivities generally require the knowledge of the displacement sensitivities. For orientation design variables the above operations require the routines resultsnoddir.f and resultsstr.c (and their subroutines). The results (i.e. the sensitivity of the von Mises stress at all nodes w.r.t. a change in an anisotropic orientation) are stored in the frd-file.

  For geometrical design variables stress_sen.c and its subroutines are needed. The evaluation of the Kreisselmeier-Steinhauser function based on the von Mises Stress is done in subroutine objective_stress_dx.f. These results are kept for further postprocessing.

• The THICKNESS objective can only be used as a constraint to another objective function. For a LE constraint the sensitivity is 1 if the actual thickness exceeds the reference thickness, else it is 0. For a GE constraint the sensitivity is 1 if the actual thickness is less than the reference thickness, else is is 0.