Note

Go to the end to download the full example code.

Flexure design

Example of the design of a flexure

It performs the topology optimization with:

  1. Maximum stiffness for the constraint modes

  2. Limited stiffness for the mechanism mode

  3. (Optional) constrained maximum use of material (volume constraint)

The example uses pymoto.SystemOfEquations to solve the mechanism and constraint modes.

References:

Koppen, S., Langelaar, M., & van Keulen, F. (2022). A simple and versatile topology optimization formulation for flexure synthesis. Mechanism and Machine Theory, 172, 104743. DOI: http://dx.doi.org/10.1016/j.mechmachtheory.2022.104743

 20 import numpy as np
 21 import pymoto as pym
 22
 23
 24 class Symmetry(pym.Module):
 25     """Enforce symmetry in the x and y direction """
 26     def __init__(self, domain: pym.VoxelDomain):
 27         self.domain = domain
 28
 29     def __call__(self, x):
 30         x = np.reshape(x, (self.domain.nely, self.domain.nelx))
 31         x = (x + np.flip(x, 0)) / 2
 32         x = (x + np.flip(x, 1)) / 2
 33         return x.flatten()
 34
 35     def _sensitivity(self, dfdv):
 36         dfdv = np.reshape(dfdv, (self.domain.nely, self.domain.nelx))
 37         dfdv = (dfdv + np.flip(dfdv, 1)) / 2
 38         dfdv = (dfdv + np.flip(dfdv, 0)) / 2
 39         return dfdv
 40
 41
 42 def flexure(nx: int = 20,
 43             ny: int = 20,
 44             doc: str = 'tx',
 45             dof: str = 'ty',
 46             emax: float = 1.0,
 47             filter_radius: float = 2.0,
 48             E: float = 100.0,
 49             nu: float = 0.3,
 50             xmin: float = 1e-9,
 51             volfrac=0.3,
 52             initial_volfrac: float = 0.2,
 53             use_symmetry=False):
 54     """Run a flexure optimization
 55     The compliance of the mechanism mode (`dof`) is maximized and the constraint mode is constrained (`doc`).
 56
 57     Args:
 58         nx (int, optional): Number of elements in x-direction. Defaults to 20.
 59         ny (int, optional): Number of elements in y-direction. Defaults to 20.
 60         doc (str, optional): Direction of constraint. Defaults to 'tx'.
 61         dof (str, optional): Direction of freedom (mechanism mode). Defaults to 'ty'.
 62         emax (float, optional): Maximum compliance value. Defaults to 1.0.
 63         filter_radius (float, optional): Filter radius. Defaults to 2.0.
 64         E (float, optional): Young's modulus. Defaults to 100.0.
 65         nu (float, optional): Poisson ratio. Defaults to 0.3.
 66         xmin (float, optional): Minimum pseudo-density value. Defaults to 1e-9.
 67         volfrac (float, optional): Volume fraction. Defaults to 0.3.
 68         initial_volfrac (float, optional): Volume fraction in the initial design. Defaults to 0.2.
 69         use_symmetry (bool, optional): Enforce symmetry in the design. Defaults to False.
 70     """
 71
 72
 73     # Set up the domain
 74     domain = pym.VoxelDomain(nx, ny)
 75
 76     # Node and dof groups
 77     nodes_top = domain.nodes[:, -1]
 78     nodes_bottom = domain.nodes[:, 0]
 79
 80     dofs_top = domain.get_dofnumber(nodes_top, [0, 1], ndof=2)
 81     dofs_bottom = domain.get_dofnumber(nodes_bottom, [0, 1], ndof=2)
 82
 83     all_dofs = np.arange(0, 2 * domain.nnodes)
 84     prescribed_dofs = np.unique(np.hstack([dofs_bottom.flatten(), dofs_top.flatten()]))
 85     free_dofs = np.setdiff1d(all_dofs, prescribed_dofs)
 86
 87     dofs_top_x = dofs_top[..., 0].flatten()
 88     dofs_top_y = dofs_top[..., 1].flatten()
 89
 90     # Construct deformation modes
 91     v = np.zeros((2 * domain.nnodes, 3), dtype=float)
 92     v[dofs_top_x, 0] = 1.0  # tx
 93     v[dofs_top_y, 1] = 1.0  # ty
 94     v[dofs_top_x, 2] = 1.0 * ny / nx  # rz
 95     v[dofs_top_y, 2] = np.linspace(1, -1, nx + 1)
 96
 97     # Select directions for constraint/mechanism
 98     degrees = ('tx', 'ty', 'rz')
 99     u = v[:, [degrees.index(doc), degrees.index(dof)]]
100
101     # Signal with design variables
102     s_x = pym.Signal('x', state=initial_volfrac * np.ones(domain.nel))
103
104     # Setup optimization problem
105     with pym.Network() as fn:
106         # Force symmetry
107         if use_symmetry:
108             s_x1 = Symmetry(domain)(s_x)
109         else:
110             s_x1 = s_x
111
112         # Density filtering
113         s_xfilt = pym.DensityFilter(domain, radius=filter_radius)(s_x1)
114         s_xfilt.tag = 'Filtered density'
115
116         # Plot the design
117         pym.PlotDomain(domain, saveto="out/design")(s_xfilt)
118
119         # SIMP penalization
120         s_xsimp = pym.MathExpression(f"{xmin} + {1 - xmin}*inp0^3")(s_xfilt)
121
122         # Assembly of stiffness matrix
123         s_K = pym.AssembleStiffness(domain, e_modulus=E, poisson_ratio=nu)(s_xsimp)
124
125         # Solve system of equations for the two loadcases
126         up = u[prescribed_dofs, :]
127         ff = np.zeros((free_dofs.size, 2))
128         s_u, s_f = pym.SystemOfEquations(free=free_dofs, prescribed=prescribed_dofs)(s_K, ff, up)
129
130         # Calculate compliances for each loadcase [constraint, mechanism]
131         s_compl = pym.EinSum('ij,ij->j')(s_u, s_f)
132
133         # Objective function: maximize compliance of mechanism mode
134         s_objective = pym.Scaling(scaling=-10)(s_compl[0])
135         s_objective.tag = "Objective"
136
137         # Compliance constraint on the constraint mode
138         s_compliance_constraint = pym.Scaling(scaling=10, maxval=emax)(s_compl[1])
139         s_compliance_constraint.tag = "Compliance constraint"
140
141         # List of optimization responses: first the objective and all others the constraints
142         responses = [s_objective, s_compliance_constraint]
143
144         # Add volume constraint if requested
145         if volfrac:
146             # Volume
147             s_volume = pym.EinSum('i->')(s_xfilt)
148
149             # Volume constraint
150             s_volume_constraint = pym.Scaling(scaling=10, maxval=volfrac * domain.nel)(s_volume)
151             s_volume_constraint.tag = "Volume constraint"
152
153             responses.append(s_volume_constraint)
154
155         pym.PlotIter()(*responses)
156
157     # Optimization
158     pym.minimize_mma(s_x, responses, function=fn, verbosity=2, maxit=200)
159
160
161 if __name__ == "__main__":
162     # flexure(100, 120, 'rz', 'ty', 0.1)  # axial spring, no rotation (flexible in y, stiff in rotation)
163     flexure(100, 100, 'tx', 'ty', 1, use_symmetry=True)  # axial spring, no shear (flexible in y, stiff in x)
164     # flexure(100, 100, 'ty', 'tx', 0.01)  # parallel guiding system (flexible in x, stiff in y)
165     # flexure(100, 100, 'rz', 'tx', 0.01)  # parallel guiding system 2 (flexible in x, stiff in rotation)
166     # flexure(100, 100, 'ty', 'rz', 0.01, volfrac=0.3)  # notch hinge (flexible in rotation, stiff in y)
167     # flexure(100, 50, 'tx', 'rz', 0.01, use_symmetry=True)  # notch hinge 2 (flexible in rotation, stiff in x)

Gallery generated by Sphinx-Gallery