Note

Go to the end to download the full example code.

Self-weight gravity

Example of the design of an arch; compliance minimization under self-weight.

Note: This is a difficult problem to optimize (in terms of convergence), because the load is dependent on the design.

With modifications on: (i) interpolation of Young’s modulus (SIMPLIN)

Consequently, in contrast to Bruyneel et al. (2005), no special treatment of the sequential approximate optimization algorithm is required.

References:

Bruyneel, M., & Duysinx, P. (2005). Note on topology optimization of continuum structures including self-weight. Structural and Multidisciplinary Optimization, 29, 245-256. DOI: https://doi.org/10.1007/s00158-004-0484-y

Zhu, J., Zhang, W., & Beckers, P. (2009). Integrated layout design of multi-component system. International Journal for Numerical Methods in Engineering, 78(6), 631-651. DOI: https://doi.org/10.1002/nme.2499

 27 import numpy as np
 28 from scipy.sparse import diags
 29 import pymoto as pym
 30
 31 # Problem settings
 32 nx, ny = 100, 50  # Domain size
 33 xmin, filter_radius = 1e-6, 1.5
 34 initial_volfrac = 1.0
 35
 36 # Material properties (Aluminium)
 37 E = 70e9  # Young's modulus (Pa)
 38 nu = 0.3
 39 rho = 2700  # kg/m^3
 40
 41 # Loadcase
 42 load = 0.0 # 100.0  # Constant point load on the top left of the domain, in negative y direction (N)
 43 gravity = np.array([0.0, -1.0]) * 9.81  # Gravity (m/s^2)
 44 # Note: Adding a point-load makes convergence easier
 45
 46 """ Choose boundary condition type
 47 1: arch
 48 2: mbb-beam
 49 3: fully-clamped
 50 4: double-arch
 51 """
 52 bc = 1
 53
 54 # Constraint on maximum volume
 55 use_volume_constraint = True
 56 max_volume = 0.2  # Maximum allowed volume fraction
 57
 58 # Constraint on maximum stress
 59 use_stress_constraint = True
 60 maximum_vm_stress = 0.2e+6  # Maximum stress value (Pa)
 61
 62
 63 if __name__ == "__main__":
 64     # Set up the domain of size (1, 1*ny/nx, 0.1)m
 65     domain = pym.VoxelDomain(nx, ny, unitx=1/nx, unity=1/nx, unitz=0.1)
 66
 67     # Node and dof groups for boundary conditions
 68     if bc == 1:  # Arch
 69         nodes_right = domain.nodes[-6:-1, 0].flatten()
 70         dofs_right = domain.get_dofnumber(nodes_right, [0, 1], 2)
 71     elif bc == 2:  # MBB-beam
 72         nodes_right = domain.nodes[-6:-1, 0].flatten()
 73         dofs_right = domain.get_dofnumber(nodes_right, 1, 2)
 74     elif bc == 3:  # Fully clamped
 75         nodes_right = domain.nodes[-1, ...].flatten()
 76         dofs_right = domain.get_dofnumber(nodes_right, [0, 1], 2)
 77     elif bc == 4:  # Double arch
 78         nodes_right = domain.nodes[-1, ...].flatten()
 79         dofs_right = domain.get_dofnumber(nodes_right, 1, 2)
 80     else:
 81         raise NotImplementedError("'bc' must be 1, 2, 3, or 4")
 82
 83     # Symmetric bc on the left (rollers)
 84     nodes_left = domain.nodes[0, ...]
 85     dofs_left_x = nodes_left*2
 86
 87     all_dofs = np.arange(0, 2 * domain.nnodes)
 88     fixed_dofs = np.unique(np.hstack([dofs_left_x.flatten(), dofs_right.flatten()]))
 89     free_dofs = np.setdiff1d(all_dofs, fixed_dofs)
 90
 91     # Setup rhs for constant load contribution
 92     f_const = np.zeros(domain.nnodes * 2)  # Generate a force vector
 93     f_const[2 * domain.nodes[0:5, -1] + 1] = -load/5
 94
 95     # Design variables
 96     s_variables = pym.Signal('x', state=initial_volfrac * np.ones(domain.nel))
 97
 98     # Setup optimization problem
 99     with pym.Network() as fn:
100         # Density filtering
101         s_filtered_variables = pym.DensityFilter(domain, radius=filter_radius)(s_variables)
102         s_filtered_variables.tag = "Filtered design"
103
104         # SIMP penalization
105         """
106         Note the use of SIMPLIN: y = xmin + (1-xmin) * (alpha * x + (alpha - 1) * x^p)
107
108         References:
109
110         Zhu, J., Zhang, W., & Beckers, P. (2009). Integrated layout design of multi-component system.
111         International Journal for Numerical Methods in Engineering, 78(6), 631-651. https://doi.org/10.1002/nme.2499
112         """
113         s_penalized_variables = pym.MathExpression(f"{xmin} + {1 - xmin}*(0.01*inp0 + 0.99*inp0^3)")(s_filtered_variables)
114
115         # Assemble stiffness matrix
116         # Apply boundary conditions by adding relatively stiff springs. This helps reduce edge stress concentrations.
117         bc_diag = np.zeros_like(f_const)
118         bc_diag[fixed_dofs] = 100 * E * np.prod(domain.element_size)
119         if hasattr(fixed_dofs, '__len__') and fixed_dofs.size > 2:
120             bc_diag[fixed_dofs[[0, -1]]] /= 2 # edge effect
121         Kfix = diags(bc_diag)
122         s_K = pym.AssembleStiffness(domain, e_modulus=E, poisson_ratio=nu, add_constant=Kfix)(s_penalized_variables)
123
124         # Determine gravity load
125         mass_el = rho * np.prod(domain.element_size)  # Mass of one element
126         s_mass = pym.MathExpression(f"{mass_el} * inp0")(s_filtered_variables)  # Scaled element mass
127         s_gravity = pym.NodalOperation(domain, np.repeat(gravity/4, 4))(s_mass)
128
129         # Add gravity load and constant load together
130         s_load0 = pym.MathExpression("inp0 + inp1")(f_const, s_gravity)
131
132         # Set the load on the boundary condition to zero
133         s_load = pym.SetValue(indices=fixed_dofs, value=0.0)(s_load0)
134
135         # Solve linear system of equations
136         s_u = pym.LinSolve()(s_K, s_load)
137
138         # Compliance
139         s_compliance = pym.EinSum('i,i->')(s_u, s_load)
140
141         # Objective
142         s_objective = pym.Scaling(scaling=1000)(s_compliance)
143         s_objective.tag = "Objective"
144
145         responses = [s_objective]
146
147         # Plotting
148         pym.PlotDomain(domain, saveto="out/design")(s_filtered_variables)
149
150         if use_stress_constraint:
151             # Calculate stress
152             s_stress = pym.Stress(domain, e_modulus=E, poisson_ratio=nu)(s_u)
153             V = np.array([[1, -0.5, 0], [-0.5, 1, 0], [0, 0, 3]])  # Vandermonde matrix
154             s_stress_vm2 = pym.EinSum("ij,ik,kj->j")(s_stress, V, s_stress)
155             s_stress_vm = pym.MathExpression('sqrt(inp0)*inp1')(s_stress_vm2, s_penalized_variables)
156
157             # Approximate maximum stress with aggregation
158             s_stress_agg = pym.PNorm(p=2,
159                                      scaling=pym.AggScaling('max', 0.3),
160                                      active_set=pym.AggActiveSet(lower_rel=0.5))(s_stress_vm)
161             s_stress_constraint = pym.Scaling(scaling=10, maxval=maximum_vm_stress)(s_stress_agg)
162             s_stress_constraint.tag = "Stress constraint"
163             responses.append(s_stress_constraint)
164
165             # Plotting
166             s_stress_scaled = pym.MathExpression("inp0/1e+6")(s_stress_vm)
167             s_stress_scaled.tag = "Von-Mises stress (MPa)"
168             pym.PlotDomain(domain, cmap='jet')(s_stress_scaled)
169
170         if use_volume_constraint:
171             # Volume
172             s_total_volume = pym.EinSum("i->")(s_filtered_variables)
173
174             # Volume constraint
175             s_volume_constraint = pym.Scaling(scaling=10, maxval=max_volume*domain.nel)(s_total_volume)
176             s_volume_constraint.tag = "Volume constraint"
177             responses.append(s_volume_constraint)
178
179         pym.PlotIter()(*responses)
180
181     # pym.finite_difference(s_variables, responses, dx=1e-4)
182
183     # Optimization (try different algorithms by uncommenting)
184     pym.minimize_slp(s_variables, responses)
185     # pym.minimize_oc(s_variables, responses[0], verbosity=2, maxit=300)  # Only unconstrained
186     # pym.minimize_mma(s_variables, responses, verbosity=2, maxit=300, move=0.1)

Gallery generated by Sphinx-Gallery