2D MBB Beam with a PDE Filter

In this example we will go over using a PDE filter instead of a convolution type filter (see Wallin et al. (2020)). For larger problems this method may scale better depending on processor counts and filter radius size. Only new material not covered in the previous example will be covered here 2D Topology Optimization with a Convolution Filter.

First there is a new variable Dc that will be the filtered sensitivity.

Listing 1: MBB Variables block

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [Dc]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = -1.0
  []
[]
(modules/combined/examples/optimization/2d_mbb_pde.i)

The AuxVariables block sensitivity is now used as a source term for the PDE filter. There is also now a Dc_elem variable that will be used for the density update.

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}
  []
  [sensitivity]
    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)"}>>> = FIRST
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = -1.0
    [AuxKernel]
      type = MaterialRealAux
      variable = sensitivity
      property = sensitivity
      execute_on = LINEAR
    []
  []
  [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}
  []
  [Dc_elem]
    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
    [AuxKernel]
      type = SelfAux
      variable = Dc_elem
      v = Dc
      execute_on = 'TIMESTEP_END'
    []
  []
[]
(modules/combined/examples/optimization/2d_mbb_pde.i)

In the Kernel block the filtering is done using a FunctionDiffusion kernel, Reaction Kernel, and a CoupledForce kernel. The function coefficient () in the FunctionDiffusion kernel is related to the radius () of the RadialAverage filter by the equation .

Listing 3: MBB Kernels block

[Kernels<<<{"href": "../../../../syntax/Kernels/index.html"}>>>]
  [diffusion]
    type = FunctionDiffusion<<<{"description": "Diffusion with a function coefficient.", "href": "../../../../source/kernels/FunctionDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = Dc
    function<<<{"description": "Function multiplier for diffusion term."}>>> = 0.15 # radius coeff
  []
  [potential]
    type = Reaction<<<{"description": "Implements a simple consuming reaction term with weak form $(\\psi_i, \\lambda u_h)$.", "href": "../../../../source/kernels/Reaction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = Dc
  []
  [source]
    type = CoupledForce<<<{"description": "Implements a source term proportional to the value of a coupled variable. Weak form: $(\\psi_i, -\\sigma v)$.", "href": "../../../../source/kernels/CoupledForce.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = Dc
    v<<<{"description": "The coupled variable which provides the force"}>>> = sensitivity
  []
[]
(modules/combined/examples/optimization/2d_mbb_pde.i)

One advantage of using a PDE filter is that by applying boundary conditions on the sensitivity variable on the boundary of the domain the filter will prevent "sticking" of the material commonly seen in topology optimization. That penalty condition is applied using the ADRobinBC where the coef controls how much the sensitivity is penalized on the boundary.

Listing 4: MBB BCs block

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [no_x]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = hold_y
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [no_y]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    value<<<{"description": "Value of the BC"}>>> = 0.0
  []
  [boundary_penalty]
    type = ADRobinBC<<<{"description": "Imposes the Robin integrated boundary condition $\\frac{\\partial u}{\\partial n}=u$.", "href": "../../../../source/bcs/ADRobinBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = Dc
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left top'
    coefficient<<<{"description": "Coefficient multiplier for the Robin boundary condition term."}>>> = 10
  []
  [boundary_penalty_right]
    type = ADRobinBC<<<{"description": "Imposes the Robin integrated boundary condition $\\frac{\\partial u}{\\partial n}=u$.", "href": "../../../../source/bcs/ADRobinBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = Dc
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'right'
    coefficient<<<{"description": "Coefficient multiplier for the Robin boundary condition term."}>>> = 10
  []
[]
(modules/combined/examples/optimization/2d_mbb_pde.i)

The UserObjects block now only contains the DensityUpdate object.

Listing 5: MBB UserObjects block

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [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_elem
    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
    force_postaux<<<{"description": "Forces the UserObject to be executed in POSTAUX"}>>> = true
  []
[]
(modules/combined/examples/optimization/2d_mbb_pde.i)

References

  1. Mathias Wallin, Niklas Ivarsson, Oded Amir, and Daniel Tortorelli. Consistent boundary conditions for pde filter regularization in topology optimization. Struct. Multidiscip. Optim., 62(3):1299–1311, sep 2020. URL: https://doi.org/10.1007/s00158-020-02556-w, doi:10.1007/s00158-020-02556-w.[BibTeX]