Note

Go to the end to download the full example code.

Stress constraint

Example of the design of cantilever for minimum volume subjected to stress constraints

It uses the following specific modules:

References:

Verbart, A., Langelaar, M., & Keulen, F. V. (2017). A unified aggregation and relaxation approach for stress-constrained topology optimization. Structural and Multidisciplinary Optimization, 55, 663-679. DOI: https://doi.org/10.1007/s00158-016-1524-0

 19 import numpy as np
 20 import pymoto as pym
 21
 22 # Problem settings
 23 nx, ny = 50, 100  # Domain size
 24 xmin, filter_radius, volfrac = 1e-9, 2, 1.0
 25
 26 maximum_vm_stress = 0.4  # Maximum allowed von-Mises stress
 27
 28 # Extra displacement constraint
 29 displacement_constraint = True
 30 max_displacement = 20.0
 31
 32
 33 class ConstraintAggregation(pym.Module):
 34     """Unified aggregation and relaxation.
 35
 36     Implemented by @artofscience (s.koppen@tudelft.nl), based on:
 37
 38     Verbart, A., Langelaar, M., & Keulen, F. V. (2017).
 39     A unified aggregation and relaxation approach for stress-constrained topology optimization.
 40     Structural and Multidisciplinary Optimization, 55, 663-679.
 41     DOI: https://doi.org/10.1007/s00158-016-1524-0
 42     """
 43
 44     def __init__(self, P=10):
 45         self.P = P
 46
 47     def __call__(self, x):
 48         """
 49         a = x + 1
 50         b = aggregation(a)
 51         c = b - 1
 52         """
 53         self.n = len(x)
 54         self.x = x
 55         self.y = self.x + 1
 56         self.z = self.y ** self.P
 57         z = (np.sum(self.z) / self.n) ** (1 / self.P)  # P-mean aggregation function
 58         return z - 1
 59
 60     def _sensitivity(self, dfdc):
 61         return (dfdc / self.n) * (np.sum(self.z) / self.n) ** (1 / self.P - 1) * self.y ** (self.P - 1)
 62
 63
 64 if __name__ == "__main__":
 65     # Set up the domain
 66     domain = pym.VoxelDomain(nx, ny)
 67
 68     # Node and dof groups
 69     nodes_left = domain.nodes[0, :]
 70     dofs_left = domain.get_dofnumber(nodes_left, [0, 1], ndof=2)
 71
 72     # Setup rhs for loadcase
 73     f = np.zeros(domain.nnodes * 2)  # Generate a force vector
 74     f[2 * domain.nodes[nx, ny//2] + 1] = 1.0
 75
 76     # Signal with design vector
 77     s_x = pym.Signal('x', state=volfrac * np.ones(domain.nel))
 78
 79     # Setup optimization problem
 80     with pym.Network() as fn:
 81         # Density filter
 82         s_xfilt = pym.DensityFilter(domain, radius=filter_radius)(s_x)
 83         s_xfilt.tag = "Filtered density"
 84
 85         # Plot the design
 86         pym.PlotDomain(domain, saveto="out/design")(s_xfilt)
 87
 88         # SIMP penalization
 89         s_xsimp = pym.MathExpression(f"{xmin} + {1 - xmin}*inp0^3")(s_xfilt)
 90
 91         # Assembly of stiffness matrix
 92         s_K = pym.AssembleStiffness(domain, e_modulus=1.0, poisson_ratio=0.3, bc=dofs_left)(s_xsimp)
 93
 94         # Solve
 95         s_displacement = pym.LinSolve()(s_K, f)
 96
 97         # Calculate stress components
 98         s_stress = pym.Stress(domain, e_modulus=1.0, poisson_ratio=0.3)(s_displacement)
 99
100         # Calculate Von-Mises stress
101         V = np.array([[1, -0.5, 0], [-0.5, 1, 0], [0, 0, 3]])  # Vandermonde matrix
102         s_stress_vm2 = pym.EinSum('ij,ik,kj->j')(s_stress, V, s_stress)
103         s_stress_vm = pym.MathExpression('sqrt(inp0)')(s_stress_vm2)
104
105         # Stress constraint
106         s_stress_constraints = pym.Scaling(maxval=maximum_vm_stress, scaling=1.0)(s_stress_vm)
107         s_stress_constraints_scaled = pym.EinSum('i,i->i')(s_xfilt, s_stress_constraints)
108         s_stress_constraint = ConstraintAggregation(P=10)(s_stress_constraints_scaled)
109         s_stress_constraint.tag = "Stress constraint"
110
111         # Volume
112         s_volume = pym.EinSum('i->')(s_xfilt)
113
114         # Objective
115         s_objective = pym.Scaling(scaling=100)(s_volume)
116         s_objective.tag = "Objective (volume)"
117
118         # Plot the stress values
119         s_stress_scaled = pym.EinSum('i,i->i')(s_xfilt, s_stress_vm)
120         s_stress_scaled.tag = 'Scaled von-Mises stress'
121         pym.PlotDomain(domain, cmap='jet')(s_stress_scaled)
122
123         # List of optimization responses: first the objective and all others the constraints
124         responses = [s_objective, s_stress_constraint]
125
126         if displacement_constraint:
127             # Output displacement
128             s_uout = pym.EinSum('i,i->')(s_displacement, f)
129
130             # Displacement constraint
131             s_uout_constraint = pym.Scaling(scaling=10, maxval=max_displacement)(s_uout)
132             s_uout_constraint.tag = "Displacement constraint"
133
134             responses.append(s_uout_constraint)
135
136         pym.PlotIter()(*responses)
137
138     # Optimization
139     pym.minimize_mma(s_x, responses, verbosity=2, maxit=300)

Gallery generated by Sphinx-Gallery