Note

Go to the end to download the full example code.

Eigenfrequency maximization

Minimal example for an eigenfrequency topology optimization

This example contains some specific modules used in dynamic problems

References:

Du, J., Olhoff, N. (2007) Topological design of freely vibrating continuum structures for maximum values of simple and multiple eigenfrequencies and frequency gaps. Structural and Multidisciplinary Optimization 34, 91-110 (2007). DOI: https://doi.org/10.1007/s00158-007-0101-y

 20 import numpy as np
 21 import pymoto as pym
 22
 23 nx, ny, nz = 50, 30, 0  # Set nz to zero for the 2D problem
 24 xmin = 1e-9
 25 filter_radius = 2.0
 26 volfrac = 0.5
 27
 28 rho = 2700.0
 29 E = 68.9e+9
 30
 31
 32 class MassInterpolation(pym.Module):
 33     """ Two-range material interpolation based on Du and Olhoff's paper
 34     For x >= threshold:
 35         y = rho * x^p1
 36     For x < threshold:
 37         y = rho * x^p0 / (t^(p0-p1))
 38     """
 39
 40     def __init__(self, rhoval=1.0, threshold=0.1, p0=6.0, p1=1.0):
 41         self.rhoval = rhoval
 42         self.threshold = threshold
 43         self.p0, self.p1 = p0, p1
 44
 45     def __call__(self, x):
 46         xx = x ** self.p1
 47         xx[x < self.threshold] = x[x < self.threshold] ** self.p0 / (self.threshold ** (self.p0 - self.p1))
 48         return self.rhoval * xx
 49
 50     def _sensitivity(self, drho):
 51         x = self.sig_in[0].state
 52         dx = self.p1 * x ** (self.p1 - 1)
 53         dx[x < self.threshold] = self.p0 * x[x < self.threshold]**(self.p0 - 1) / (self.threshold**(self.p0 - self.p1))
 54         return self.rhoval * dx * drho
 55
 56
 57 if __name__ == "__main__":
 58     print(__doc__)
 59
 60     if nz == 0:  # 2D structural eigenfrequency analysis
 61         # Generate a grid
 62         domain = pym.VoxelDomain(nx, ny)
 63
 64         # Generate a non-design area that has mass
 65         nondesign_area = domain.elements[3*nx//4:, ny//4:ny*3//4].flatten()
 66     else:  # 3D
 67         # Generate a grid
 68         domain = pym.VoxelDomain(nx, ny, nz)
 69
 70         # Generate a non-design area that has mass
 71         nondesign_area = domain.elements[3*nx//4:, ny//4:ny*3//4, nz//4:nz*3//4].flatten()
 72
 73     # Calculate boundary dof indices
 74     boundary_nodes = domain.nodes[0, ...]
 75     boundary_dofs = domain.get_dofnumber(boundary_nodes, ndof=domain.dim)
 76
 77     # Safety, to prevent overloading the memory in your machine; no problem to remove this
 78     if domain.nnodes > 1e+6:
 79         print("Too many nodes :(")
 80         exit()
 81
 82     # Make design vector, and fill with initial values
 83     sx = pym.Signal('x', state=np.ones(domain.nel) * volfrac)
 84
 85     # Start building the modular network
 86     with pym.Network() as fn:
 87         # Density filter
 88         sxfilt = pym.DensityFilter(domain, radius=filter_radius)(sx)
 89
 90         # Set the non-design domain
 91         sxndd = pym.SetValue(indices=nondesign_area, value=1.0)(sxfilt)
 92
 93         sx_analysis = sxndd  # Alias for the design to perform the analysis on
 94
 95         # Show the design on the screen as it optimizes
 96         pym.PlotDomain(domain, saveto="out/design", clim=[0, 1])(sx_analysis)
 97
 98         # SIMP material interpolation
 99         # Note: Material properties can either be set in the scaling variables (sSIMP and sDENS; as is done here), or in
100         # the assembly modules AssembleStiffness and AssembleMass by providing the relevant keyword arguments.
101         sSIMP = pym.MathExpression(f"{E}*({xmin} + {1.0 - xmin}*inp0^3)")(sx_analysis)
102         sDENS = MassInterpolation(rhoval=rho)(sx_analysis)
103
104         # Assemble mass and stiffness matrix
105         sK = pym.AssembleStiffness(domain, bc=boundary_dofs)(sSIMP)
106         sM = pym.AssembleMass(domain, bc=boundary_dofs, ndof=domain.dim)(sDENS)
107
108         # Eigenvalue solver
109         slams, seigvec = pym.EigenSolve(hermitian=True, nmodes=3)(sK, sM)
110         slams.tag = "eigenvalues"
111         seigvec.tag = "eigenvectors"
112
113         # Output the design and eigenmodes to a Paraview file
114         pym.WriteToVTI(domain, saveto='out/dat.vti')(sx_analysis, seigvec)
115
116         # Get harmonic mean of three lowest eigenvalues
117         sharm = pym.MathExpression('1/inp0 + 1/inp1 + 1/inp2')(slams[0], slams[1], slams[2])
118         sharm.tag = "harmonic mean"
119
120         # MMA needs correct scaling of the objective
121         sg0 = pym.Scaling(scaling=100.0)(sharm)
122         sg0.tag = "objective"
123
124         # Calculate the volume of the domain by adding all design densities together
125         svol = pym.EinSum('i->')(sx_analysis)
126         svol.tag = 'volume'
127
128         # Volume constraint; note that also pym.Scaling(scaling=10.0, maxval=domain.nel*volfrac)(svol) could be used
129         sg1 = pym.MathExpression(f'10*(inp0/{domain.nel} - {volfrac})')(svol)
130         sg1.tag = "volume constraint"
131
132         # Maybe you want to check the design-sensitivities with finite difference?
133         do_finite_difference = False
134         if do_finite_difference:
135             pym.finite_difference(sx, [sg0, sg1], dx=1e-4)
136             exit()
137
138         pym.PlotIter()(sg0, sg1)  # Plot iteration history
139         pym.ScalarToFile('out/log.csv')(sharm, slams, svol)
140
141     # Do the optimization with MMA
142     pym.minimize_mma(sx, [sg0, sg1], verbosity=2)
143
144     # Here you can do some post processing
145     print("The optimization has finished!")
146     print(f"The final eigenfrequencies obtained are {np.sqrt(slams.state)}")

Gallery generated by Sphinx-Gallery