Note

Go to the end to download the full example code.

Dynamic compliance

Example of the design of cantilever for minimum dynamic compliance and shows usage of complex values

This example contains some specific modules used in dynamic problems

  • pymoto.AssembleMass To assemble the mass matrix

  • pymoto.AddMatrix For addition of mass and stiffness matrices (with complex components) to form the dynamic stiffness matrix

  • pymoto.ComplexNorm To calculate the norm of a complex value, corresponding to amplitude of vibration

References:

Silva, O. M., Neves, M. M., & Lenzi, A. (2019). A critical analysis of using the dynamic compliance as objective function in topology optimization of one-material structures considering steady-state forced vibration problems. Journal of Sound and Vibration, 444, 1-20. DOI: https://doi.org/10.1016/j.jsv.2018.12.030

 20 import numpy as np
 21 import pymoto as pym
 22
 23 # Problem settings
 24 lx, ly = 1, 0.5  # Domain size
 25 ny = 50
 26 nx = int(lx/ly)*ny
 27 unitx, unity = lx / nx, ly / ny
 28 unitz = 1.0  # out-of-plane thickness
 29
 30 xmin, filter_radius, volfrac = 1e-6, 2, 0.49  # Density settings
 31
 32 E, nu, rho = 210e9, 0.3, 7860  # Properties of steel in SI units [Pa], [-], [kg/m3]
 33
 34 # fundamental eigenfreq is 372.6 Hz, above this frequency the optimization will not result in a nice design
 35 omega = 370 * (2 * np.pi)
 36
 37 # force
 38 force_magnitude = -9000
 39
 40 # Rayleigh damping parameters
 41 alpha, beta = 1e-3, 1e-8
 42
 43
 44 if __name__ == "__main__":
 45     # Set up the domain
 46     domain = pym.VoxelDomain(nx, ny, unitx=unitx, unity=unity, unitz=unitz)
 47
 48     # Node and dof groups
 49     nodes_left = domain.nodes[0, :]
 50     dofs_left = domain.get_dofnumber(nodes_left, [0, 1], ndof=2)
 51
 52     # Setup rhs for loadcase
 53     f = np.zeros(domain.nnodes * 2)  # Generate a force vector
 54     f[2 * domain.get_nodenumber(nx, ny // 2) + 1] = force_magnitude
 55
 56     # Design vector signal
 57     s_x = pym.Signal('x', state=volfrac * np.ones(domain.nel))
 58
 59     # Setup optimization problem
 60     with pym.Network() as fn:
 61         # Density filtering
 62         s_xfilt = pym.DensityFilter(domain, radius=filter_radius)(s_x)
 63         s_xfilt.tag = 'Filtered density'
 64
 65         # Plot the design
 66         pym.PlotDomain(domain, saveto="out/design")(s_xfilt)
 67
 68         # SIMP penalization
 69         s_xsimp = pym.MathExpression(f"{xmin} + {1 - xmin}*inp0^3")(s_xfilt)
 70
 71         # Assembly of stiffness matrix
 72         s_K = pym.AssembleStiffness(domain, e_modulus=E, poisson_ratio=nu, bc=dofs_left)(s_xsimp)
 73
 74         # Assemble mass matrix
 75         s_M = pym.AssembleMass(domain, bc=dofs_left, material_property=rho, ndof=2)(s_xsimp)
 76
 77         # Calculate the eigenfrequencies only once
 78         calculate_eigenfrequencies = True
 79         if calculate_eigenfrequencies:
 80             s_eig, _ = pym.EigenSolve(hermitian=True, nmodes=3)(s_K, s_M)
 81
 82             eigfreq = np.sqrt(s_eig.state)
 83             print(f"Eigenvalues are {eigfreq} rad/s or {eigfreq / (2 * np.pi)} Hz")
 84
 85         # Build dynamic stiffness matrix Z = K + iω(αM + βK) - ω^2 M
 86         s_Z = pym.AddMatrix()(1 + 1j*omega*beta, s_K, -omega**2 + 1j*omega*alpha, s_M)
 87
 88         # Solve linear system of equations
 89         s_u = pym.LinSolve()(s_Z, f)
 90
 91         # Output displacement (is a complex value)
 92         s_cdyn = pym.EinSum('i,i->')(s_u, f)
 93
 94         # Absolute value (amplitude of response)
 95         s_ampl = pym.ComplexNorm()(s_cdyn)
 96
 97         # Objective
 98         s_objective = pym.Scaling(scaling=100.0)(s_ampl)
 99         s_objective.tag = "Objective"
100
101         # Volume
102         s_volume = pym.EinSum('i->')(s_xfilt)
103
104         # Volume constraint
105         s_volume_constraint = pym.Scaling(scaling=10.0, maxval=volfrac * domain.nel)(s_volume)
106         s_volume_constraint.tag = "Volume constraint"
107
108         # List of optimization responses: first the objective and all others the constraints
109         responses = [s_objective, s_volume_constraint]
110
111         # Show iteration history
112         pym.PlotIter()(*responses)
113
114     # Optimization
115     pym.minimize_mma(s_x, responses, verbosity=2)

Gallery generated by Sphinx-Gallery