Source code for openaerostruct.structures.failure_ks
importnumpyasnpimportopenmdao.apiasom
[docs]classFailureKS(om.ExplicitComponent):""" Aggregate failure constraints from the structure. To simplify the optimization problem, we aggregate the individual elemental failure constraints using a Kreisselmeier-Steinhauser (KS) function. The KS function produces a smoother constraint than using a max() function to find the maximum point of failure, which produces a better-posed optimization problem. The rho inputeter controls how conservatively the KS function aggregates the failure constraints. A lower value is more conservative while a greater value is more aggressive (closer approximation to the max() function). parameters ---------- vonmises[ny-1, 2] : numpy array von Mises stress magnitudes for each FEM element. Returns ------- failure : float KS aggregation quantity obtained by combining the failure criteria for each FEM node. Used to simplify the optimization problem by reducing the number of constraints. """definitialize(self):self.options.declare("surface",types=dict)self.options.declare("rho",types=float,default=100.0)defsetup(self):surface=self.options["surface"]rho=self.options["rho"]ifsurface["fem_model_type"]=="tube":num_failure_criteria=2elifsurface["fem_model_type"]=="wingbox":num_failure_criteria=4self.ny=surface["mesh"].shape[1]self.add_input("vonmises",val=np.zeros((self.ny-1,num_failure_criteria)),units="N/m**2")self.add_output("failure",val=0.0)self.sigma=surface["yield"]self.rho=rhoself.declare_partials("*","*")defcompute(self,inputs,outputs):sigma=self.sigmarho=self.rhovonmises=inputs["vonmises"]fmax=np.max(vonmises/sigma-1)nlog,nsum,nexp=np.log,np.sum,np.expks=1/rho*nlog(nsum(nexp(rho*(vonmises/sigma-1-fmax))))outputs["failure"]=fmax+ksdefcompute_partials(self,inputs,partials):vonmises=inputs["vonmises"]sigma=self.sigmarho=self.rho# Find the location of the max stress constraintfmax=np.max(vonmises/sigma-1)i,j=np.where((vonmises/sigma-1)==fmax)i,j=i[0],j[0]# Set incoming seed as 1 so we simply get the jacobian entriesksb=1.0# Use results from the AD code to compute the jacobian entriestempb0=ksb/(rho*np.sum(np.exp(rho*(vonmises/sigma-fmax-1))))tempb=np.exp(rho*(vonmises/sigma-fmax-1))*rho*tempb0fmaxb=ksb-np.sum(tempb)# Populate the entriesderivs=tempb/sigmaderivs[i,j]+=fmaxb/sigma# Reshape and save them to the jac dictpartials["failure","vonmises"]=derivs.reshape(1,-1)