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]
[MeshGenerator]
type = GeneratedMeshGenerator
dim = 2
nx = 150
ny = 50
xmin = 0
xmax = 30
ymin = 0
ymax = 10
[]
[node]
type = ExtraNodesetGenerator
input = MeshGenerator
new_boundary = pull
nodes = 0
[]
[push]
type = ExtraNodesetGenerator
input = node
new_boundary = push
coord = '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]
[Emin]
family = MONOMIAL
order = CONSTANT
initial_condition = ${Emin}
[]
[power]
family = MONOMIAL
order = CONSTANT
initial_condition = ${power}
[]
[E0]
family = MONOMIAL
order = CONSTANT
initial_condition = ${E0}
[]
[Dc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${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]
[elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
youngs_modulus = E_phys
poissons_ratio = poissons_ratio
args = 'Emin mat_den power E0'
[]
[E_phys]
type = DerivativeParsedMaterial
# Emin + (density^penal) * (E0 - Emin)
expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
coupled_variables = 'mat_den'
property_name = E_phys
[]
[poissons_ratio]
type = GenericConstantMaterial
prop_names = poissons_ratio
prop_values = 0.3
[]
[stress]
type = ComputeLinearElasticStress
[]
[dc]
type = ComplianceSensitivity
design_density = mat_den
youngs_modulus = E_phys
incremental = 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]
[rad_avg]
type = RadialAverage
radius = 1.2
weights = linear
prop_name = sensitivity
execute_on = TIMESTEP_END
force_preaux = true
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = TIMESTEP_BEGIN
[]
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
execute_on = TIMESTEP_END
force_postaux = true
[]
[]
(modules/combined/examples/optimization/2d_mbb.i)(modules/combined/examples/optimization/2d_mbb.i)
vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 2
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[MeshGenerator]
type = GeneratedMeshGenerator
dim = 2
nx = 150
ny = 50
xmin = 0
xmax = 30
ymin = 0
ymax = 10
[]
[node]
type = ExtraNodesetGenerator
input = MeshGenerator
new_boundary = pull
nodes = 0
[]
[push]
type = ExtraNodesetGenerator
input = node
new_boundary = push
coord = '30 10 0'
[]
[]
[AuxVariables]
[Emin]
family = MONOMIAL
order = CONSTANT
initial_condition = ${Emin}
[]
[power]
family = MONOMIAL
order = CONSTANT
initial_condition = ${power}
[]
[E0]
family = MONOMIAL
order = CONSTANT
initial_condition = ${E0}
[]
[Dc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${vol_frac}
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = SMALL
add_variables = true
incremental = false
[]
[]
[BCs]
[no_x]
type = DirichletBC
variable = disp_y
boundary = pull
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.0
[]
[]
[NodalKernels]
[pull]
type = NodalGravity
variable = disp_y
boundary = push
gravity_value = -1
mass = 1
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
youngs_modulus = E_phys
poissons_ratio = poissons_ratio
args = 'Emin mat_den power E0'
[]
[E_phys]
type = DerivativeParsedMaterial
# Emin + (density^penal) * (E0 - Emin)
expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
coupled_variables = 'mat_den'
property_name = E_phys
[]
[poissons_ratio]
type = GenericConstantMaterial
prop_names = poissons_ratio
prop_values = 0.3
[]
[stress]
type = ComputeLinearElasticStress
[]
[dc]
type = ComplianceSensitivity
design_density = mat_den
youngs_modulus = E_phys
incremental = false
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
[rad_avg]
type = RadialAverage
radius = 1.2
weights = linear
prop_name = sensitivity
execute_on = TIMESTEP_END
force_preaux = true
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = TIMESTEP_BEGIN
[]
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
execute_on = TIMESTEP_END
force_postaux = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
nl_abs_tol = 1e-8
dt = 1.0
num_steps = 70
[]
[Outputs]
[out]
type = Exodus
execute_on = 'TIMESTEP_END'
[]
print_linear_residuals = false
[]
[Postprocessors]
[total_vol]
type = ElementIntegralVariablePostprocessor
variable = mat_den
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
(modules/combined/examples/optimization/2d_mbb.i)
vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 2
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[MeshGenerator]
type = GeneratedMeshGenerator
dim = 2
nx = 150
ny = 50
xmin = 0
xmax = 30
ymin = 0
ymax = 10
[]
[node]
type = ExtraNodesetGenerator
input = MeshGenerator
new_boundary = pull
nodes = 0
[]
[push]
type = ExtraNodesetGenerator
input = node
new_boundary = push
coord = '30 10 0'
[]
[]
[AuxVariables]
[Emin]
family = MONOMIAL
order = CONSTANT
initial_condition = ${Emin}
[]
[power]
family = MONOMIAL
order = CONSTANT
initial_condition = ${power}
[]
[E0]
family = MONOMIAL
order = CONSTANT
initial_condition = ${E0}
[]
[Dc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${vol_frac}
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = SMALL
add_variables = true
incremental = false
[]
[]
[BCs]
[no_x]
type = DirichletBC
variable = disp_y
boundary = pull
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.0
[]
[]
[NodalKernels]
[pull]
type = NodalGravity
variable = disp_y
boundary = push
gravity_value = -1
mass = 1
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
youngs_modulus = E_phys
poissons_ratio = poissons_ratio
args = 'Emin mat_den power E0'
[]
[E_phys]
type = DerivativeParsedMaterial
# Emin + (density^penal) * (E0 - Emin)
expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
coupled_variables = 'mat_den'
property_name = E_phys
[]
[poissons_ratio]
type = GenericConstantMaterial
prop_names = poissons_ratio
prop_values = 0.3
[]
[stress]
type = ComputeLinearElasticStress
[]
[dc]
type = ComplianceSensitivity
design_density = mat_den
youngs_modulus = E_phys
incremental = false
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
[rad_avg]
type = RadialAverage
radius = 1.2
weights = linear
prop_name = sensitivity
execute_on = TIMESTEP_END
force_preaux = true
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = TIMESTEP_BEGIN
[]
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
execute_on = TIMESTEP_END
force_postaux = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
nl_abs_tol = 1e-8
dt = 1.0
num_steps = 70
[]
[Outputs]
[out]
type = Exodus
execute_on = 'TIMESTEP_END'
[]
print_linear_residuals = false
[]
[Postprocessors]
[total_vol]
type = ElementIntegralVariablePostprocessor
variable = mat_den
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
(modules/combined/examples/optimization/2d_mbb.i)
vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 2
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[MeshGenerator]
type = GeneratedMeshGenerator
dim = 2
nx = 150
ny = 50
xmin = 0
xmax = 30
ymin = 0
ymax = 10
[]
[node]
type = ExtraNodesetGenerator
input = MeshGenerator
new_boundary = pull
nodes = 0
[]
[push]
type = ExtraNodesetGenerator
input = node
new_boundary = push
coord = '30 10 0'
[]
[]
[AuxVariables]
[Emin]
family = MONOMIAL
order = CONSTANT
initial_condition = ${Emin}
[]
[power]
family = MONOMIAL
order = CONSTANT
initial_condition = ${power}
[]
[E0]
family = MONOMIAL
order = CONSTANT
initial_condition = ${E0}
[]
[Dc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${vol_frac}
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = SMALL
add_variables = true
incremental = false
[]
[]
[BCs]
[no_x]
type = DirichletBC
variable = disp_y
boundary = pull
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.0
[]
[]
[NodalKernels]
[pull]
type = NodalGravity
variable = disp_y
boundary = push
gravity_value = -1
mass = 1
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
youngs_modulus = E_phys
poissons_ratio = poissons_ratio
args = 'Emin mat_den power E0'
[]
[E_phys]
type = DerivativeParsedMaterial
# Emin + (density^penal) * (E0 - Emin)
expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
coupled_variables = 'mat_den'
property_name = E_phys
[]
[poissons_ratio]
type = GenericConstantMaterial
prop_names = poissons_ratio
prop_values = 0.3
[]
[stress]
type = ComputeLinearElasticStress
[]
[dc]
type = ComplianceSensitivity
design_density = mat_den
youngs_modulus = E_phys
incremental = false
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
[rad_avg]
type = RadialAverage
radius = 1.2
weights = linear
prop_name = sensitivity
execute_on = TIMESTEP_END
force_preaux = true
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = TIMESTEP_BEGIN
[]
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
execute_on = TIMESTEP_END
force_postaux = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
nl_abs_tol = 1e-8
dt = 1.0
num_steps = 70
[]
[Outputs]
[out]
type = Exodus
execute_on = 'TIMESTEP_END'
[]
print_linear_residuals = false
[]
[Postprocessors]
[total_vol]
type = ElementIntegralVariablePostprocessor
variable = mat_den
execute_on = 'INITIAL TIMESTEP_END'
[]
[]
(modules/combined/examples/optimization/2d_mbb.i)
vol_frac = 0.5
E0 = 1
Emin = 1e-8
power = 2
[GlobalParams]
displacements = 'disp_x disp_y'
[]
[Mesh]
[MeshGenerator]
type = GeneratedMeshGenerator
dim = 2
nx = 150
ny = 50
xmin = 0
xmax = 30
ymin = 0
ymax = 10
[]
[node]
type = ExtraNodesetGenerator
input = MeshGenerator
new_boundary = pull
nodes = 0
[]
[push]
type = ExtraNodesetGenerator
input = node
new_boundary = push
coord = '30 10 0'
[]
[]
[AuxVariables]
[Emin]
family = MONOMIAL
order = CONSTANT
initial_condition = ${Emin}
[]
[power]
family = MONOMIAL
order = CONSTANT
initial_condition = ${power}
[]
[E0]
family = MONOMIAL
order = CONSTANT
initial_condition = ${E0}
[]
[Dc]
family = MONOMIAL
order = CONSTANT
initial_condition = -1.0
[]
[mat_den]
family = MONOMIAL
order = CONSTANT
initial_condition = ${vol_frac}
[]
[]
[Physics/SolidMechanics/QuasiStatic]
[all]
strain = SMALL
add_variables = true
incremental = false
[]
[]
[BCs]
[no_x]
type = DirichletBC
variable = disp_y
boundary = pull
value = 0.0
[]
[no_y]
type = DirichletBC
variable = disp_x
boundary = right
value = 0.0
[]
[]
[NodalKernels]
[pull]
type = NodalGravity
variable = disp_y
boundary = push
gravity_value = -1
mass = 1
[]
[]
[Materials]
[elasticity_tensor]
type = ComputeVariableIsotropicElasticityTensor
youngs_modulus = E_phys
poissons_ratio = poissons_ratio
args = 'Emin mat_den power E0'
[]
[E_phys]
type = DerivativeParsedMaterial
# Emin + (density^penal) * (E0 - Emin)
expression = '${Emin} + (mat_den ^ ${power}) * (${E0}-${Emin})'
coupled_variables = 'mat_den'
property_name = E_phys
[]
[poissons_ratio]
type = GenericConstantMaterial
prop_names = poissons_ratio
prop_values = 0.3
[]
[stress]
type = ComputeLinearElasticStress
[]
[dc]
type = ComplianceSensitivity
design_density = mat_den
youngs_modulus = E_phys
incremental = false
[]
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[UserObjects]
[rad_avg]
type = RadialAverage
radius = 1.2
weights = linear
prop_name = sensitivity
execute_on = TIMESTEP_END
force_preaux = true
[]
[update]
type = DensityUpdate
density_sensitivity = Dc
design_density = mat_den
volume_fraction = ${vol_frac}
execute_on = TIMESTEP_BEGIN
[]
[calc_sense]
type = SensitivityFilter
density_sensitivity = Dc
design_density = mat_den
filter_UO = rad_avg
execute_on = TIMESTEP_END
force_postaux = true
[]
[]
[Executioner]
type = Transient
solve_type = NEWTON
petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
petsc_options_value = 'lu superlu_dist'
nl_abs_tol = 1e-8
dt = 1.0
num_steps = 70
[]
[Outputs]
[out]
type = Exodus
execute_on = 'TIMESTEP_END'
[]
print_linear_residuals = false
[]
[Postprocessors]
[total_vol]
type = ElementIntegralVariablePostprocessor
variable = mat_den
execute_on = 'INITIAL TIMESTEP_END'
[]
[]