COMSOL 6.3 - Topology Optimization of an MBB Beam

Topology optimization in a structural mechanics context can answer the question: Given that you know the loads on the structure, which distribution of the available material maximizes stiffness? Or, conversely, how much material is necessary to obtain a predefined stiffness, and how should it be distributed? Such investigations typically occur during the concept design stages.

The conflicting goals of stiffness maximization and mass minimization lead to a continuum of possible optimal solutions, depending on how you balance the goals against each other. This topology optimization example demonstrates how to use a penalization method (SIMP) to obtain an optimal distribution of a fixed amount of material such that stiffness is maximized. Changing the amount of material available leads to a different solution which is also Pareto optimal, representing a different balance between the conflicting objectives.

The geometry and load pattern is similar to the Messerschmitt–Bölkow–Blohm beam (MBB) used as a validation test for topology optimization.

This application is also used in Introduction to the Optimization Module.

Model Definition

The model studies optimal material distribution in the beam, which consists of a linear elastic material, structural steel. The dimensions of the beam region — 6-meters by 1-meter by 0.5 meters — means an original total weight of 23,550 kg. The beam is symmetric about the plane x = 3. Both its ends rest on rollers while an edge load acts on the top middle part (Figure 1).

Figure 1: Geometry of the beam with loads and constraints.

The optimality criterion is defined by the objective function, which is chosen to be the total strain energy in this example. Note that the strain energy exactly balances the work done by the applied load, so minimizing the strain energy minimizes the displacement induced at the points where loads are applied, effectively minimizing the compliance of the structure — maximizing its stiffness. The other, conflicting, objective is minimization of total mass, which is implemented as an upper bound on the mass of the optimized structure.

The design can be described by an artificial material volume factor, θc, which acts as control variable in the optimization problem. In structural mechanics, convention dictates that it is 0 in void domains and 1 in solid domains. In order to make the problem formulation well-posed one has to limit the complexity of the design; otherwise the optimization will just tell us that the optimal structure is a composite with infinitesimal features, which cannot be resolved with any mesh. The detail of the design can be limited using a filter that produces another field, θf, which is guaranteed to be free of features smaller than some set value, Rmin because it is computed as a solution to a Helmholtz equation (see Ref. 1):

The equation essentially just smears out the design as illustrated below with and θc and θf to the left and right, respectively.

Figure 2: The material volume factor (right) has a lot of gray scale.

This filtered design variable can then be projected to reduce the extent of the gray scale region. COMSOL supports projection based on the hyperbolic tangent function (see Ref. 2):

This is illustrated for θβ = 0.5 and β = 8 below with θf and θ to the left and right.

Figure 3: The (projected) material volume factor shown to the right has no visible gray scale.

High values of β are associated with strong projection and less gray scale, but also slows down the optimization process, so initially we will not use it and thus set θ = θf.

To prevent a completely solid design, we impose an overall limit Vfrac on the fraction θavg of the domain area occupied by solid material

To guarantee a binary solution we use the SIMP penalization method (see Ref. 3) to define the penalized field

Since the stiffness must for numerical reasons not vanish completely anywhere in the model, the minimum penalized volume fraction, θmin, has been introduced. The exponent p ≥ 1 is a penalty factor that makes intermediate densities provide less stiffness compared to what they cost in weight. Increasing p therefore forces θc toward either of its bounds leading to a sharper solution.

Ideally, the topology of the resulting designs should be mesh independent, and for designs with moderate complexity this can be achieved, but if the optimal design is very complex, other designs with slightly different topologies perform only marginally worse, so the optimization problem has many local minima, and it is likely that a slightly suboptimal design is identified.

Finally, it is worth noting that the problem of finding a stiff design can be quite different from finding a design that does not break and it is generally advised to use shape or parameter optimization to ensure that the design is free of excessive stress concentrations.

Results

The following plot shows the stiffness distribution for the optimized solution. The resulting design is an approximation of a truss structure, which is expected for this beam size. Note that the optimal design for a longer beam is quite different.

Figure 4: The penalized stiffness for the optimized design.

Notes About the COMSOL Implementation

The control variable field, Helmholtz filter and SIMP penalization are defined through the use of a feature. A Solid Mechanics interface represents the structural properties of the beam, while the objective and the constraint is added directly in the study step. The elastic strain energy density is a predefined variable, solid.Ws, available to use as the objective function for the optimization problem.

Only the left half of the beam geometry needs to be modeled because of the assumed symmetry, which is implemented as a Roller condition on the symmetry plane. Plots for postprocessing and inspection of the solution while solving are set up using a Mirror dataset, to show the complete solution. The geometry is parameterized, making it easy to experiment with different beam sizes, but keep in mind that the mesh size is taken as the default filter radius.

The optimization solver is selected and controlled from an Optimization study step. For topology optimization, either the MMA or the SNOPT solver can be used. They have different merits and weaknesses. MMA tends to be braver in the beginning, proceeding quickly toward an approximate optimum, while SNOPT is more cautious but also converges more efficiently close to the final solution thanks to its second-order approximation of the objective.

This example demonstrates a strategy based on using MMA with a limited number of iterations on successively finer meshes (see the option in the ). 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.

References

1. B.S. Lazarov and O. Sigmund, “Filters in topology optimization based on Helmholtz-type differential equations,” Int. J. Numer. Methods Eng., vol. 86, pp. 765–781, 2011.

2. F. Wang, B.S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” Struct. Multidiscipl. Optim., vol. 43, pp. 767–784, 2011.

3. M.P. Bendsøe, “Optimal shape design as a material distribution problem,” Structural Optimization, vol. 1, pp. 193–202, 1989.

Optimization_Module/Topology_Optimization/mbb_beam_optimization

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

Click

In the tree, select > > .

Click

Global Definitions

Parameters 1

In the window, expand the > node, then click > .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
a	3 [m]	3 m	Beam half width
b	1 [m]	1 m	Beam height
L1	0.1 [m]	0.1 m	Support width
meshsz	6 [cm]	0.06 m	Mesh size
volfrac	0.5	0.5	Maximum volume fraction
beta	1	1	Projection slope

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type a.

In the text field, type b.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type L1.

Point 2 (pt2)

In the toolbar, click

In the window for , locate the section.

In the text field, type a-L1/2.

In the text field, type b.

In the toolbar, click

Load Boundary

In the toolbar, click

and choose .

In the window for , type Load Boundary in the text field.

Locate the section. From the list, choose .

Locate the section. In the text field, type a-L1/1.999.

In the text field, type b*0.999.

Locate the section. From the list, choose .

Symmetry y

In the toolbar, click

and choose .

In the window for , type Symmetry y in the text field.

Locate the section. From the list, choose .

Locate the section. In the text field, type a*0.001.

In the text field, type L1*1.001.

Locate the section. From the list, choose .

Symmetry x

In the toolbar, click

and choose .

In the window for , type Symmetry x in the text field.

Locate the section. From the list, choose .

Locate the section. In the text field, type a*0.999.

Locate the section. From the list, choose .

In the toolbar, click

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click in the window toolbar.

In the toolbar, click

to close the window.

Materials

Topology Link 1 (toplnk1)

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Solid Mechanics (solid)

Roller 1

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

Prescribed Displacement 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

This is effectively a roller condition along the x-axis, but it is applied on a vertical boundary to avoid bending stiffness.

Boundary Load 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Since this is a linear problem, the magnitude of the applied force does not affect the optimal topology.

Specify the Ftot vector as


0	x
-100 [kN]	y

Mesh 1

In the window, under right-click and choose .

Size

In the window, under > click .

In the window for , locate the section.

Click the button.

Locate the section. In the text field, type meshsz.

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
volfrac	0.5	0.5	Maximum volume fraction
beta	1	1	Projection slope

Here we have also defined the maximum volume fraction.

Add a density topology feature, which can be used to distinguish between void and solid regions. This variable is coupled back to the 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.

Topology Optimization

Density Model 1 (dtopo1)

In the window, under > click .

In the window for , locate the section.

From the list, choose .

In the β text field, type beta.