Ordered SIMP: Multimaterial

In addition to establish a maximum volume fraction and optimize the distribution of material over the domain, one can optimize multiple materials within the same domain and assign, additionally, individual cost to each of them (see Tavakoli and Mohseni (2014)). This allows to enforce another constraint: The total cost over the domain.

This approach uses the object DensityUpdateTwoConstraints, which allows multiple material optimization to be applied to blocks individually. The materials in this example have pseudo-densities of 0.4, 0.7, and 1.0; in addition to the void material of 0.0.

Cost and Young's modulus are defined as interpolated quantities across the three materials:

Listing 1: MBB Material interpolation for multimaterial optimization

[Materials]
  [elasticity_tensor]
    type = ComputeVariableIsotropicElasticityTensor
    youngs_modulus = E_phys
    poissons_ratio = poissons_ratio
    args = 'mat_den'
  []
  [E_phys]
    type = DerivativeParsedMaterial
    # ordered multimaterial simp
    expression = "A1:=(${E0}-${E1})/(${rho0}^${power}-${rho1}^${power}); B1:=${E0}-A1*${rho0}^${power}; E1:=A1*mat_den^${power}+B1; A2:=(${E1}-${E2})/(${rho1}^${power}-${rho2}^${power}); B2:=${E1}-A2*${rho1}^${power}; E2:=A2*mat_den^${power}+B2; A3:=(${E2}-${E3})/(${rho2}^${power}-${rho3}^${power}); B3:=${E2}-A3*${rho2}^${power}; E3:=A3*mat_den^${power}+B3; if(mat_den<${rho1},E1,if(mat_den<${rho2},E2,E3))"
    coupled_variables = 'mat_den'
    property_name = E_phys
  []

  [Cost_mat]
    type = DerivativeParsedMaterial
    # ordered multimaterial simp
    expression = "A1:=(${C0}-${C1})/(${rho0}^(1/${power})-${rho1}^(1/${power})); B1:=${C0}-A1*${rho0}^(1/${power}); C1:=A1*mat_den^(1/${power})+B1; A2:=(${C1}-${C2})/(${rho1}^(1/${power})-${rho2}^(1/${power})); B2:=${C1}-A2*${rho1}^(1/${power}); C2:=A2*mat_den^(1/${power})+B2; A3:=(${C2}-${C3})/(${rho2}^(1/${power})-${rho3}^(1/${power})); B3:=${C2}-A3*${rho2}^(1/${power}); C3:=A3*mat_den^(1/${power})+B3; if(mat_den<${rho1},C1,if(mat_den<${rho2},C2,C3))"
    coupled_variables = 'mat_den'
    property_name = Cost_mat
  []
  [CostDensity]
    type = ParsedMaterial
    property_name = CostDensity
    coupled_variables = 'mat_den Cost'
    expression = 'mat_den*Cost'
  []
  [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
  []
  [cc]
    type = CostSensitivity
    design_density = mat_den
    cost = Cost_mat
    outputs = 'exodus'
  []
[]
(modules/combined/examples/optimization/three_materials.i)

The optimization process and the enforcement of volume and cost constraints are driven by the DensityUpdateTwoConstraints, radial average, and sensitivity filter user objects:

Listing 2: MBB User objects for multimaterial optimization

[UserObjects]
  [rad_avg]
    type = RadialAverage
    radius = 3
    weights = linear
    prop_name = sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [rad_avg_cost]
    type = RadialAverage
    radius = 3
    weights = linear
    prop_name = cost_sensitivity
    execute_on = TIMESTEP_END
    force_preaux = true
  []
  [update]
    type = DensityUpdateTwoConstraints
    # This is
    density_sensitivity = Dc
    cost_density_sensitivity = Cc
    cost = Cost
    cost_fraction = ${cost_frac}
    design_density = mat_den
    volume_fraction = ${vol_frac}
    bisection_lower_bound = 0
    bisection_upper_bound = 1.0e16 # 100
    bisection_move = 0.05
    adaptive_move = true
    relative_tolerance = 1.0e-3
    execute_on = TIMESTEP_BEGIN
  []
  # Provides Dc
  [calc_sense]
    type = SensitivityFilter
    density_sensitivity = Dc
    design_density = mat_den
    filter_UO = rad_avg
    execute_on = TIMESTEP_END
    force_postaux = true
  []
  # Provides Cc
  [calc_sense_cost]
    type = SensitivityFilter
    density_sensitivity = Cc
    design_density = mat_den
    filter_UO = rad_avg_cost
    execute_on = TIMESTEP_END
    force_postaux = true
  []
[]
(modules/combined/examples/optimization/three_materials.i)

The final distribution of materials for the bridge structure (note that symmetry is being used) looks as follows:

References

  1. Rouhollah Tavakoli and Seyyed Mohammad Mohseni. Alternating active-phase algorithm for multimaterial topology optimization problems: a 115-line matlab implementation. Structural and Multidisciplinary Optimization, 49:621–642, 2014.[BibTeX]