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 (809) 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 mafillsmmain_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}$$ (934)

  
  

  is calculated.

• For static linear calculations

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

  
  

  is determined.

• For frequency and Green calculations

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

  
  

  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 sens_coor.c.

For the objective $G$ the total sensitivity $\frac{d G(s,U(s))}{d s}$ is written as $\frac{\partial G}{\partial s} + \frac{\partial G}{\partial U} \cdot \frac{\partial U}{\partial s}$. So the objective function is used in the terms $\frac{\partial G}{\partial s}$ and $\frac{\partial G}{\partial U}$. The routine objectivemain_se.c 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 STRAIN ENERGY as objective function a distinction has to be made whether the calculation is geometrically linear or nonlinear. For a linear geometric calculation Equation (810) 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).$$ (937)

  
  

  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 (810) 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).$$ (938)

  
  

  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.$$ (939)

  
  

  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 (814), which can also be written as:

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

  
  

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

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

  
  

  leads to the following expressions for $c_j$:

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

  
  

  and $c_i=0$. The latter equation results from the differentiation of the mass normalization condition $U_i^T M U_i=1$. Indeed, one obtains:

  $$2 \left ( \frac{\partial U_i}{\partial s} \right) ^T M U_i =0,$$ (943)

  
  

  since $M$ is symmetric and not dependent on the orientation of the material. This means that $\left ( \frac{\partial U_i}{\partial s} \right) ^T$ is orthogonal to $U_i$. Mathematically, this can also be interpreted as follows: the solution of Equation (940) can be seen as a particular solution complemented by an arbitrary multiple of the homogeoneous solution. The homogeneous solution is $U_i$, a particular solution is obtained by a linear combination of the other eigenmodes, which are orthogonal to $U_i$ (since K and M are symmetric). Now, Equation (943) expresses that the solution is orthogonal to $U_i$, so only the particular solution with coefficients from Equation (942) remains.

  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$$ (944)

  
  

  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.$$ (945)

  
  

  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 sensi_coor.c.

• The ALL-DISP 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),$$ (946)

  
  

  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)$$ (947)

  
  

  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 MISESSTRESS, PS1STRESS or PS3STRESS as objective are immediately based on Equation (810). Here too, one can write an equation of the form:

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

  
  

  i.e. the stress S is in general 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.

  For geometrical design variables both terms $\partial S/\partial s$ and $\partial S/\partial U$ have to be evaluated, i.e. keeping the displacements at the nodes constant while changing the geometry will lead to stress changes, as well as changing the displacements while keeping the geometry constant.

  For orientation design variables $\partial S/\partial s=0$ since the stress change only comes through the material law.

  For geometrical design variables the stress in the unperturbed state is calculated in resultsstr.c, while the derivative w.r.t. s is done in stress_sen_dx.f and the derivative w.r.t. U in stress_sen_dv.f. The calculation of Equation (810) is done exactly as explained in Section 6.9.23 from left to right, i.e. first calculating $(\partial G/\partial s) K^{-1}$ before continuing with the expression in brackets.

  For orientation design variables the sensivity 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},$$ (949)

  
  

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

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

  
  

  So the stress sensitivities require the knowledge of the displacement sensitivities. This was treated in the previous item. 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.

• Results for MODALSTRESS as objective are also based on Equation (804), but now Equation (941) is used for $\partial U_i/\partial s$. Right now MODALSTRESS can only be used for geometric sensitivities. Therefore, ${\partial M/\partial s \ne 0}$ and Equation (943) has to be replaced by:

  $$2 \frac{\partial U_i}{\partial s} ^ T M U_i = - U_i^T \frac{\partial M}{\partial s} U_i,$$ (951)

  
  

  and $c_i$ now satisfies:

  $$c_i = - \frac{1}{2} U_i^T \frac{\partial M}{\partial s} U_i,$$ (952)

  
  

  instead of being zero. Defining

  $$\alpha_j = \frac{\partial S}{\partial U_i} U_j,$$ (953)

  
  

  one obtains for the stress sensitivity for mode $i$:

  $$\frac{d S}{d s} = \frac{\partial S}{\partial s} + c_i \alpha_i + ... ...bda_i -\lambda_j} \right) \left[ U_j^T \frac{\partial F_i}{\partial s} \right].$$ (954)