2D MBB Beam with a Convolution Filter

In this example we will go through the general setup of a topology optimization problem using solid isotropic material penalization (SIMP), see Sigmund (2001). The problem is to find the optimal material distribution in a 2D domain that minimizes the compliance of the structure. We will first define the problem parameters. Below is a list of the parameters that we will use in this example corresponding to the volume fraction, Young's modulus of the material, and the penalization power. We will go over the material that is needed for the setup and running the optimization problem, but skip over any information that is covered in other MOOSE tutorials.


vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 3

Next we define the mesh and add the necessary nodesets.

Listing 1: MBB Mesh block

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  [MeshGenerator]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
    nx<<<{"description": "Number of elements in the X direction"}>>> = 150
    ny<<<{"description": "Number of elements in the Y direction"}>>> = 50
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = 0
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 30
    ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = 0
    ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 10
  []
  [node]
    type = ExtraNodesetGenerator<<<{"description": "Creates a new node set and a new boundary made with the nodes the user provides.", "href": "../../../../source/meshgenerators/ExtraNodesetGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = MeshGenerator
    new_boundary<<<{"description": "The names of the boundaries to create"}>>> = pull
    nodes<<<{"description": "The nodes you want to be in the nodeset (Either this parameter or \"coord\" must be supplied)."}>>> = 0
  []
  [push]
    type = ExtraNodesetGenerator<<<{"description": "Creates a new node set and a new boundary made with the nodes the user provides.", "href": "../../../../source/meshgenerators/ExtraNodesetGenerator.html"}>>>
    input<<<{"description": "The mesh we want to modify"}>>> = node
    new_boundary<<<{"description": "The names of the boundaries to create"}>>> = push
    coord<<<{"description": "The nodes with coordinates you want to be in the nodeset. Separate multple coords with ';' (Either this parameter or \"nodes\" must be supplied)."}>>> = '30 10 0'
  []
[]
(modules/combined/examples/optimization/2d_mbb.i)

In the AuxVariables block there are two initial conditions. The first is a constant, negative value that is needed for the sensitivity variable Dc. It needs to be negative for the first density update. The second initial condition is setting the material density to the initial value of vol_frac.

Listing 2: MBB AuxVariables block

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]

  [Emin]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = ${Emin}
  []
  [power]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = ${power}
  []
  [E0]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = ${E0}
  []
  [Dc]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = -1.0
  []
  [mat_den]
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = ${vol_frac}
  []
[]
(modules/combined/examples/optimization/2d_mbb.i)

The next block is the Materials block. The first material is the "SIMP" density altered young's modulus material. The material follows the form E = Emin + (density^penal) * (E0 - Emin). The second material is the compliance sensitivity, which is used for updating the density field.

Listing 3: MBB Materials block

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor<<<{"description": "Compute an isotropic elasticity tensor for elastic constants that change as a function of material properties", "href": "../../../../source/materials/ComputeVariableIsotropicElasticityTensor.html"}>>>
    youngs_modulus<<<{"description": "Name of material property defining the Young's Modulus"}>>> = E_phys
    poissons_ratio<<<{"description": "Name of material property defining the Poisson's Ratio"}>>> = poissons_ratio
    args<<<{"description": "Variable dependence for the Young's Modulus and Poisson's Ratio materials"}>>> = 'Emin mat_den power E0'
  []
  [E_phys]
    type = DerivativeParsedMaterial<<<{"description": "Parsed Function Material with automatic derivatives.", "href": "../../../../source/materials/DerivativeParsedMaterial.html"}>>>
    # Emin + (density^penal) * (E0 - Emin)
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'mat_den'
    property_name<<<{"description": "Name of the parsed material property"}>>> = E_phys
  []
  [poissons_ratio]
    type = GenericConstantMaterial<<<{"description": "Declares material properties based on names and values prescribed by input parameters.", "href": "../../../../source/materials/GenericConstantMaterial.html"}>>>
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = poissons_ratio
    prop_values<<<{"description": "The values associated with the named properties"}>>> = 0.3
  []
  [stress]
    type = ComputeLinearElasticStress<<<{"description": "Compute stress using elasticity for small strains", "href": "../../../../source/materials/ComputeLinearElasticStress.html"}>>>
  []
  [dc]
    type = ComplianceSensitivity<<<{"description": "Computes compliance sensitivity needed for SIMP method.", "href": "../../../../source/materials/ComplianceSensitivity.html"}>>>
    design_density<<<{"description": "Design density variable name."}>>> = mat_den
    youngs_modulus<<<{"description": "DerivativeParsedMaterial for Youngs modulus."}>>> = E_phys
    incremental<<<{"description": "Optional flag for error checking if an incremental or total model should be used."}>>> = false
  []
[]
(modules/combined/examples/optimization/2d_mbb.i)

The final block is the UserObjects block. This block contains the main optimization functionality. First is the RadialAverage and SensitivityFilter objects that filter the sensitivity to prevent checkerboarding (see Sigmund (2007)). The radius of the filter sets the minimum size of a feature in the structure. Finally is the DensityUpdate object that updates the density field based on the sensitivity and the keeps the volume constraint satisfied.

Listing 4: MBB UserObjects block

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [rad_avg]
    type = RadialAverage<<<{"description": "Perform a radial average of a material property", "href": "../../../../source/userobjects/RadialAverage.html"}>>>
    radius<<<{"description": "Cut-off radius for the averaging"}>>> = 1.2
    weights<<<{"description": "Distance based weight function"}>>> = linear
    prop_name<<<{"description": "Name of the material property to average"}>>> = sensitivity
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = TIMESTEP_END
    force_preaux<<<{"description": "Forces the UserObject to be executed in PREAUX"}>>> = true
  []
  [update]
    type = DensityUpdate<<<{"description": "Compute updated densities based on sensitivities using an optimality criteria method to keep the volume constraint satisified.", "href": "../../../../source/userobjects/DensityUpdate.html"}>>>
    density_sensitivity<<<{"description": "Name of the density_sensitivity variable."}>>> = Dc
    design_density<<<{"description": "Design density variable name."}>>> = mat_den
    volume_fraction<<<{"description": "Volume Fraction"}>>> = ${vol_frac}
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = TIMESTEP_BEGIN
  []
  [calc_sense]
    type = SensitivityFilter<<<{"description": "Computes the filtered sensitivity using a radial average user object.", "href": "../../../../source/userobjects/SensitivityFilter.html"}>>>
    density_sensitivity<<<{"description": "Name of the density_sensitivity variable."}>>> = Dc
    design_density<<<{"description": "Design density variable name."}>>> = mat_den
    filter_UO<<<{"description": "Radial Average user object"}>>> = rad_avg
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = TIMESTEP_END
    force_postaux<<<{"description": "Forces the UserObject to be executed in POSTAUX"}>>> = true
  []
[]
(modules/combined/examples/optimization/2d_mbb.i)

References

  1. Ole Sigmund. A 99 line topology optimization code written in matlab. Structural and multidisciplinary optimization, 21:120–127, 2001.[BibTeX]
  2. Ole Sigmund. Morphology-based black and white filters for topology optimization. Structural and Multidisciplinary Optimization, 33:401–424, 04 2007. doi:10.1007/s00158-006-0087-x.[BibTeX]