PDF

Topology Optimization of a Loaded Knee Structure
Introduction
Imagine that you are designing a bracket that must carry a specified load offset down and to the side of its anchoring. As a first attempt, you try a fully solid construction like the one below.
Figure 1: Geometry of the structure with loads and constraints.
This design is clearly overdesigned and does not make efficient use of material. As a direct simulation reveals, stresses are unevenly distributed, meaning that some areas carry significantly less load than others. Provided that the stiffness of the part can be reduced to, for example, 40% of the value computed from the fully solid design without violating the specification, it is tempting to remove some material — saving both weight and costs. But what is the minimum amount of material necessary, and how should it be distributed?
This model demonstrates how you can apply the SIMP model  (Solid Isotropic Material with Penalization) for structural topology optimization with COMSOL Multiphysics. Note that, in contrast to the way SIMP is usually applied, this example minimizes weight for a specified minimum stiffness, instead of maximizing stiffness for a specified maximum weight.
Model Definition
The area that may contain solid material is bounded by the initial L-shaped design; see Figure 1. The structure’s uppermost boundary is fixed while the load is carried evenly distributed over the outermost 5 cm of the boundary indicated in the figure. The material from which the structure is made is a structural steel with Young’s modulus, E0 = 200 GPa, and Poisson’s ratio, ν = 0.33. The part is to be cut from a sheet of 2 cm thickness.
In the SIMP model (see Ref. 1), the stress tensor is considered to be a function of the true Young’s modulus E0 and the material volume factor, θp, which is a function of the control variable, 0 ≤ θc ≤ 1:
Because the stiffness must for numerical reasons not vanish completely anywhere in the model, the penalized material volume factor, θp, is constrained such that 10-3 = θmin ≤ θp 1. The exponent p ≥ 1 is a penalty factor that makes intermediate densities provide less stiffness compared to what they cost in weight. The material volume factor is related to the control variable via a Helmholtz filter and a projection function based on the hyperbolic tangent function:
Here θf, β, and θβ are the filtered material volume factor, projection slope and projection point. The solution can be made independent of the mesh resolution by using a fixed value of the filter radius, Rmin, but in this example we will not change the mesh and thus apply the default filter radius based on the local mesh element size.
When there is a lower limit on the stiffness of any acceptable design, expressed equivalently as an upper bound on the total strain energy
(1)
the SIMP penalization forces θc toward either of its bounds. Increasing p therefore leads to a sharper solution. The maximum strain energy Wsmax is conveniently expressed as a factor times the lowest possible strain energy, as computed for the fully solid design in Figure 1.
Conceptually, the optimization problem is simply minimizing the amount of material used, measured as the solid fraction of the initial design domain
(2)
without violating the stiffness constraint in Equation 1. Here, A is the area of the design domain.
The Finite Element Mesh
Without regularization, problems of this type are intrinsically mesh-size dependent. For this model, generate an initial free triangular mesh using a maximum element size of 15 mm. This mesh, shown in Figure 2, is rather coarse and effectively limits the smallest possible size of details in the topology solution. When the mesh is later refined and the optimization solver restarted on the finer mesh, proper scaling of the regularization term ensures that the updated solution retains the same topology, only sharpening the details.
Figure 2: Finite element mesh of the structure.
Results and Discussion
Figure 3 shows one possible material distribution, optimal for the chosen regularization strategy. Black areas represent material and white areas represent void. Unless composite materials are considered, gray areas are unphysical so a useful optimized design should be essentially black and white.
Figure 3: Distribution of the control variable after optimization.
As can be seen in Figure 4, below, the stress is relatively evenly distributed over the solid material. Each individual member in the truss-like structure is subjected to roughly the same maximum stress, indicating that the design makes full use of the material’s stiffness. At the same time, the peak stress is not excessive: only about 3 times the mean stress in the solid areas.
Figure 4: The von Mises stress of the optimal distribution of the control variable.
Reference
1. M.P. Bendsøe and O. Sigmund, Topology Optimization Theory, Methods, and Applications, Springer, 2004.
Notes About the COMSOL Implementation
A Solid Mechanics interface is combined with a Density Model to represent the structural properties of the bracket, while the objective and constraint are defined directly in a Topology Optimization study step. The internal strain energy is a predefined variable, solid.Ws, readily available to use as the objective function for the optimization problem.
The optimization solver is selected and controlled from the study step. For topology optimization, either the MMA, IPOPT or the SNOPT solver can be used. This example demonstrates a strategy based on using MMA with a limited number of iterations on successively finer meshes. It quickly and reliably produces trustworthy topologies showing good improvement of the objective function value. These solutions are not necessarily globally optimal, which may, however, be of less importance in practice.
Application Library path: Optimization_Module/Topology_Optimization/loaded_knee
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select Structural Mechanics>Solid Mechanics (solid).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies>Stationary.
6
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Note that the reference strain energy is initially set equal to 1 and must be updated before optimization is started.
Geometry 1
Create the geometry. To simplify this step, insert a prepared geometry sequence.
1
In the Geometry toolbar, click Insert Sequence and choose Insert Sequence.
2
3
In the Geometry toolbar, click  Build All.
4
Click the  Zoom Extents button in the Graphics toolbar.
The geometry should now look like that in Figure 1. Note that the inserted geometry is parameterized and that the parameters used are automatically added to the list of global parameters in the model.
Add Material
1
In the Home toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Built-in>Structural steel.
4
Click Add to Global Materials in the window toolbar.
5
In the Home toolbar, click  Add Material to close the Add Material window.
Materials
Material Link 1 (matlnk1)
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose More Materials>Material Link.
Add a Density Model feature, which can be used to distinguish between void and solid regions. This variable is coupled back to the Solid Mechanics interface via the Young’s modulus. The filter radius should not be smaller than the mesh element size, so the default will work, but a fixed value has to be chosen to make the result of the optimization mesh independent.
Definitions
In the Definitions toolbar, click  Optimization and choose Topology Optimization>Density Model.
Topology Optimization
Density Model 1 (dtopo1)
1
2
In the Settings window for Density Model, locate the Projection section.
3
From the Projection type list, choose Hyperbolic tangent projection.
4
In the β text field, type beta.
Add a Mass Properties node to compute total mass and area.
Definitions
Mass Properties 1 (mass1)
1
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose Physics Utilities>Mass Properties.
2
In the Settings window for Mass Properties, locate the Density section.
3
From the Density source list, choose From physics interface.
Solid Mechanics (solid)
1
In the Model Builder window, under Component 1 (comp1) click Solid Mechanics (solid).
2
In the Settings window for Solid Mechanics, locate the 2D Approximation section.
3
4
Locate the Thickness section. In the d text field, type d.
Materials
Topology Link 1 (toplnk1)
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose More Materials>Topology Link.
2
In the Settings window for Topology Link, locate the Link Settings section.
3
From the Topology source list, choose Density Model 1 (dtopo1).
4
Locate the Geometric Entity Selection section. From the Selection list, choose All domains.
In the next steps, you fix the structure for rigid body translation and rotation by prescribing zero displacements along one of the sides, then apply the load.
Solid Mechanics (solid)
Fixed Constraint 1
1
In the Physics toolbar, click  Boundaries and choose Fixed Constraint.
2
Boundary Load 1
1
In the Physics toolbar, click  Boundaries and choose Boundary Load.
2
3
In the Settings window for Boundary Load, locate the Force section.
4
From the Load type list, choose Total force.
5
Specify the Ftot vector as
Mesh 1
Free Triangular 1
In the Mesh toolbar, click  Free Triangular.
Size
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
Click the Custom button.
4
Locate the Element Size Parameters section. In the Maximum element size text field, type h0.
5
Click  Build All.
The resulting mesh is shown in Figure 2.
Study 1
First compute only a stationary solution on which you can evaluate a reference value for the objective, as well as set up suitable plots.
1
In the Home toolbar, click  Compute.
Results
Stress (solid)
Click the  Zoom Extents button in the Graphics toolbar.
Penalized Material Volume Factor
1
In the Model Builder window, expand the Results>Topology Optimization node, then click Output material volume factor.
2
In the Settings window for 2D Plot Group, type Penalized Material Volume Factor in the Label text field.
Surface 1
Plotting the Penalized Material Volume Factor gives a sharper and fairer picture of the topology, as seen by the Solid Mechanics interface.
1
In the Model Builder window, expand the Penalized Material Volume Factor node, then click Surface 1.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Definitions>Density Model 1>Auxiliary variables>dtopo1.theta_p - Penalized material volume factor.
3
Locate the Coloring and Style section. From the Color table list, choose GrayScale.
4
From the Color table transformation list, choose Reverse.
5
In the Penalized Material Volume Factor toolbar, click  Plot.
Global Evaluation 1
1
In the Model Builder window, expand the Results>Derived Values node.
2
Right-click Results>Derived Values and choose Global Evaluation.
3
In the Settings window for Global Evaluation, click Add Expression in the upper-right corner of the Expressions section. From the menu, choose Component 1 (comp1)>Solid Mechanics>Global>solid.Ws_tot - Total elastic strain energy - J.
4
Click  Evaluate.
Use the value from the Table 1 window to update the value of the parameter WsRef. This will normalize the compliance constraint in such a way that a maximum allowed relative compliance increase can be easily prescribed.
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Set up probes for the objective functions and constraint.
Definitions
Global Variable Probe 1 (var1)
1
In the Definitions toolbar, click  Probes and choose Global Variable Probe.
2
In the Settings window for Global Variable Probe, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Definitions>Density Model 1>Global>dtopo1.theta_avg - Average material volume factor.
3
Locate the Expression section. Select the Description check box.
4
Global Variable Probe 2 (var2)
1
In the Definitions toolbar, click  Probes and choose Global Variable Probe.
2
In the Settings window for Global Variable Probe, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1)>Solid Mechanics>Global>solid.Ws_tot - Total elastic strain energy - J.
3
Locate the Expression section. In the Expression text field, type solid.Ws_tot/WsRef.
4
Select the Description check box.
5
Study 1
Add the Topology Optimization study step, select a solver, and limit the number of iterations. Switch on plotting while solving.
Topology Optimization
1
In the Study toolbar, click  Optimization and choose Topology Optimization.
2
In the Settings window for Topology Optimization, locate the Optimization Solver section.
3
In the Maximum number of iterations text field, type 50.
This value is low enough to terminate the optimization process before the optimality tolerance limit has been reached, resulting in a warning from the solver. However, as you will be able to verify from the probes you just set up, it is sufficiently large for the objective function value to reach a stable level.
4
Click Add Expression in the upper-right corner of the Objective Function section. From the menu, choose Component 1 (comp1)>Definitions>Density Model 1>Global>comp1.dtopo1.theta_avg - Average material volume factor.
5
Click Add Expression in the upper-right corner of the Constraints section. From the menu, choose Component 1 (comp1)>Solid Mechanics>Global>comp1.solid.Ws_tot - Total elastic strain energy - J.
6
Locate the Constraints section. In the table, enter the following settings:
7
Locate the Output While Solving section. Select the Plot check box.
8
From the Plot group list, choose Penalized Material Volume Factor.
9
In the Model Builder window, click Study 1.
10
In the Settings window for Study, locate the Study Settings section.
11
Clear the Generate default plots check box.
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
4
5
Locate the Output While Solving section. From the Probes list, choose None.
6
Click to expand the Advanced Settings section. Select the Reuse solution from previous step check box.
7
In the Study toolbar, click  Compute.
Results
Penalized Material Volume Factor
Click the  Zoom Extents button in the Graphics toolbar.
The default plot shows the stress distribution in the optimized structure, clearly indicating that the load paths follow the optimized topology (Figure 4). Switch to the other plot group to display a sharper picture of the final design (Figure 3).