Checkerboard-free topology optimization for compliance minimization of continuum elastic structures based on the generalized finite-volume theory

Abstract

Topology optimization is a well-suited method to establish the best material distribution inside an analysis domain. It is common to observe some numerical instabilities in its gradient-based version, such as the checkerboard pattern, mesh dependence, and local minima. This research demonstrates the generalized finite-volume theory's checkerboard-free property by performing topology optimization algorithms without filtering techniques. The formation of checkerboard regions is associated with the finite element method's displacement field assumptions, where the equilibrium and continuity conditions are satisfied through the element nodes. On the other hand, the generalized finite-volume theory satisfies the continuity conditions between common faces of adjacent subvolumes, which is more likely from the continuum mechanics point of view. Also, the topology optimization algorithms based on the generalized finite-volume theory are performed using a mesh independent filter that regularizes the subvolume sensitivities, providing optimum topologies that avoid the mesh dependence and length scale issues.

Keywords:
Topology optimization; checkerboard-free approach; generalized finite-volume theory; continuum elastic structures; and finite element method

Graphical Abstract

1 INTRODUCTION

In structural design, engineers seek to find the best project that attends all the design restrictions and optimizes structural performance. Currently, the best project is accomplished based on the engineer experience, causing dependence on their work. Therefore, structural optimization techniques have been developed to help engineers find the optimal configuration for structural designs, without the need to base their designs on past experiences. In general, structural optimization problems are divided into two main categories: material optimization and material distribution optimization. The first category intends to establish the best material properties to a design, while the second seeks to find the best material distribution inside an analysis domain.

The material distribution-based optimization problems include sizing optimization, which seeks to find the optimal size in terms of length, thickness, and highness; shape optimization, which introduces shape changes on the design to find the optimal solution; and topology optimization, which seeks to find the best material distribution inside the analysis domain for the given objective function and constraints. The structural topology optimization is one of the most important structural optimization problems, becoming one of the fastest-growing research fields in the structural analysis due to its applications in different areas, such as solid mechanics, physics, multi-material modeling, and computer sciences.

The topology optimization problem was proposed initially by Michell (1904), who derived the Optimality Criteria (OC) method for the least weight layout of trusses. However, this method is typically used for compliance minimization or stiffness maximization problems, usually combined with the so-called SIMP (Solid Isotropic Material with Penalization) approach. In this approach, the interest is in determining the best solid isotropic material distribution on the analysis domain (Bendsøe and Sigmund, 2003). Therefore, the material properties are modeled by the relative material density raised to a given power to penalize the intermediate values. In topology optimization, the SIMP method has been extensively used due to its versatility, convergence, and ease implementation (Rozvany, 2009).

Topology optimization has raised as a powerful technique for structural design, although there are some problems related to numerical instabilities. The main numerical problems are the checkerboard pattern effect, which refers to the formation of regions alternating solid and void elements in a checkerboard shape; mesh dependence, which refers to the problem of not having qualitatively the same solution for different discretizations; and local minima, which refers to the problem of having different solutions for the same discretizations when different input parameters are employed. It is undesirable to have any of these instabilities in the optimal solution.

Since the pioneering work of Bendsøe and Kikuchi (1988) in the homogenization method, the finite element-based strategy for structural topology optimization has received full attention and experienced considerable progress (Wang and Wang, 2006). Therefore, the advantages and disadvantages are well-known. For instance, according to Díaz and Sigmund (1995), the checkerboard pattern is directly associated with the finite element method numerical assumptions, which leads to some artificial stiffness. Different approaches lead efficiently with this problem, as the adoption of higher-order finite elements (Díaz and Sigmund, 1995; Sigmund and Petersson, 1998), filtering techniques based on image processing or perimeter control (Sigmund, 2007; Haber et al., 1996) and the employment of modified finite elements (Rozvany et al., 2003; Pomezanski et al., 2005; Poulsen, 2002).

An alternative technique to the finite element method is the finite-volume theory, which employs the volume-average of the different fields that define the material behavior and imposes the boundary and continuity conditions in an averaged sense. This technique has shown to be a well suitable method for elastic stress analysis in solid mechanics, investigations of its numerical efficiency can be found in Cavalcante et al. (2007a, b and 2008) and Cavalcante and Pindera (2012a, b). The satisfaction of equilibrium equations at the subvolume level, concomitant to kinematic and static continuities established in a surface-averaged sense between common faces of adjacent subvolumes, are features that distinguish the finite-volume theory from the finite element method. More recently, Araujo et al. (2020) have employed a topology optimization technique for compliance minimization based on the standard finite-volume theory to obtain checkerboard free topologies in the absence of filtering techniques, also providing more efficient topologies when a mesh-independent filter is employed.

The checkerboard instability mentioned later is related to the finite element method's assumptions, such as the satisfaction of equilibrium and continuity conditions at the element nodes. Also, the equilibrium equations are not satisfied at the element level, only when a sufficiently refined mesh is employed. Differently, the finite-volume theory satisfies the equilibrium equations at the subvolume level, and the compatibility conditions are established through the subvolume interfaces. Thus, in the finite-volume theory, the connections between adjacent subvolumes occur through subvolumes' faces, which is more likely from the continuum mechanics point of view. In the finite element method, the connections between neighboring elements occur through the nodes, leading to optimum topologies with checkerboard regions in the absence of regularization techniques for trilateral or quadrilateral elements.

This paper addresses a new approach for topology optimization based on the generalized finite-volume theory for continuum elastic structures in the context of compliance minimization problems, showing that the checkerboard pattern is a problem related to the conventional finite element analysis. Two different ways to evaluate the objective function in the context of the higher-order versions of the finite-volume theory are investigated as essential guidance for this technique's employment. Comparison results between the three versions of the generalized finite-volume theory and similar approaches based on the finite element method are provided, demonstrating the efficiency of the new topology optimization technique, with competitive processing time, even when the higher-order versions of the theory are employed.

2 GENERALIZED FINITE-VOLUME THEORY

The finite-volume method is a well-suited numerical technique for boundary-value problems in fluid mechanics governed by parabolic and hyperbolic equations (Versteeg and Malalasekera, 2007). The formulations of this method in solid mechanics are characterized by differences in the subvolume displacement representation and the domain discretization and the local satisfaction of differential equilibrium equations (Cavalcante and Pindera, 2012a).

The structural finite-volume theory has its origins in the so-called higher-order theory for functionally graded materials, developed in a sequence of papers in the 1990s, and summarized in Aboudi et al. (1999). A reconstruction of this theory is firstly suggested by Bansal and Pindera (2003 and 2005) and Zhong et al. (2004). They have simplified the design domain discretization and implemented an efficient local/global stiffness matrix approach. Therefore, this reconstruction has revealed the new higher-order approach as, in fact, a finite-volume method, motivating the nomenclature changing to reflect the aspects of the reconstructed theory fundamentally. After that, Cavalcante et al. (2007a, b) introduced a parametric mapping in the elasticity-based version of the finite-volume theory, enabling the modeling of curved structures. Following Cavalcante et al. (2007a, b), Gattu et al. (2008) and Khatam and Pindera (2009 and 2010) suggested a parametric mapping of the homogenized version of the finite-volume theory, known as FVDAM (Finite-Volume Direct Averaging Micromechanics).

However, the second-order displacement field representation inside the subvolumes and the enforcement of tractions and displacements in a surface-averaged sense leads to interpenetrations between common faces of adjacent subvolumes (Cavalcante and Pindera, 2012a). As a result, Cavalcante and Pindera (2012a) suggested a generalization of the finite-volume theory, based on a higher-order displacement field representation. They have introduced new surface-averaged kinematic and static variables, inspired on the linear elasticity theory assumptions, preserving the finite-volume framework, as the local satisfaction of equilibrium equations and the establishment of continuity conditions in a surface-averaged sense. Thus, the additional coefficients of the displacement field can be expressed in terms of the new surface-averaged kinematic variables, which enforces continuity across adjacent subvolumes, avoiding undesirable interfacial interpenetrations.

The generalization proposed by Cavalcante and Pindera (2012a, b) is applicable for rectangular analysis domains discretized in rectangular subvolumes. This generalization is accomplished by adding systematically different orders to the zeroth-order (standard) finite-volume theory, which corresponds to the original version presented by Bansal and Pindera (2003). Each order corresponds to an increase in the displacement field complexity, followed by the addition of kinematic quantities evaluated in an average sense at the subvolume faces. Thus, the first order finite-volume theory incorporates rotations to the original version, while the second-order finite-volume theory includes rotations and curvatures. Cavalcante and Pindera (2014a, b) proposed a generalization of the homogenized version of the finite-volume theory for periodic materials under finite deformations.

Recently, Chen et al. (2018) proposed a three-dimensional parametric formulation of the FVDAM theory for multiphase heterogeneous materials with periodic microstructure. Similarly, Vieira and Marques (2019) have proposed a parametric three-dimensional extension of the finite-volume theory to evaluate the thermal conductivity of periodic multiphase composites. Summarily, the finite-volume theory is quite a new numerical approach, mainly employed for heterogeneous materials with periodic microstructures, which is an excellent solution for the checkerboard pattern issue usually presented in topology optimization for compliance minimization based on the finite element method.

Different versions of the finite-volume method can be found in the literature, as the cell-centered and vertex-centered approaches (Cavalcante and Pindera, 2012a). They can share similar features with the finite-volume theory, as the satisfaction of the equilibrium equations locally, and the continuity conditions imposed through the faces, as expected from a continuum mechanics point of view. These features can also be found in the Discrete Element Method (DEM) for continuous medium. See, for example, one of the most recent applications of the DEM approach to the multi-scale modeling of heterogeneous materials with periodic microstructure (Ferretti, 2020). These approaches can also be explored to solve the checkerboard pattern of topology optimization.

2.1 Theoretical framework

The presented formulation has its roots in the second-order version of the generalized finite-volume theory presented in Cavalcante and Pindera (2012a). This technique approximates the displacement field by second-order Legendre polynomials expressed as a function of the local coordinates inside each subvolume (Cavalcante et al., 2007a). Besides, the boundary and continuity conditions are imposed in a surface-averaged sense, and the equilibrium equations are satisfied at the subvolume level.

Figure 1 shows the adopted rectangular domain in $x_1-x_2$ plane with $0\le x_1\le L$ and $0\le x_2\le H$, which is discretized in $N_{\beta }$ horizontal subvolumes and $N_{\gamma }$ vertical subvolumes. The subvolume dimensions are $l_q$ and $h_q$ for $q=\mathrm{1,}\dots ,N_q$, where $N_q=N_{\beta }\cdot N_{\gamma }$ is the total number of subvolumes. In this cartesian formulation of the generalized finite-volume theory, the components of the displacement field can be approximated by the Legendre polynomial expansion in the local coordinate system, shown in Figure 1, Cavalcante and Pindera (2012a):

$u_i^{\left(q\right)}=W_{i\left(00\right)}^{\left(q\right)}+x_1^{\left(q\right)}{W}_{i\left(10\right)}^{\left(q\right)}+x_2^{\left(q\right)}{W}_{i\left(01\right)}^{\left(q\right)}+x_1^{\left(q\right)}{x}_2^{\left(q\right)}{W}_{i\left(11\right)}^{\left(q\right)}+\frac{1}{2}\left(3{x_1^{\left(q\right)}}^2-\frac{l_q^2}{4}\right)W_{i\left(20\right)}^{\left(q\right)}+\frac{1}{2}\left(3{x_2^{\left(q\right)}}^2-\frac{h_q^2}{4}\right)W_{i\left(02\right)}^{\left(q\right)}+\frac{1}{2}\left(3{x_1^{\left(q\right)}}^2-\frac{l_q^2}{4}\right)x_2^{\left(q\right)}{W}_{i\left(21\right)}^{\left(q\right)}+\frac{1}{2}\left(3{x_2^{\left(q\right)}}^2-\frac{h_q^2}{4}\right)x_1^{\left(q\right)}{W}_{i\left(12\right)}^{\left(q\right)}+\frac{1}{4}\left(3{x_1^{\left(q\right)}}^2-\frac{l_q^2}{4}\right)\left(3{x_2^{\left(q\right)}}^2-\frac{h_q^2}{4}\right)W_{i\left(22\right)}^{\left(q\right)}$ (1)

where $i=\mathrm{1,2}$ and $W_{i\left(mn\right)}^{\left(q\right)}$ are the unknown coefficients of the displacement field. These coefficients are expressed as a function of the following kinematic quantities: surface-averaged displacements, rotations, and curvatures, which are responsible for determining the generalized stiffness matrices (Cavalcante and Pindera, 2012a).

Figure 1:
Discretized analysis domain and local coordinate system of a generic subvolume $q$.

Figure 2:
Surface-averaged quantities: (a) surface-averaged kinematic quantities and (b) surface-averaged static quantities.

Figure 2a illustrates the kinematic quantities associated with each face of a generic subvolume $q$. Therefore, these quantities in terms of surface-averaged displacements, rotations, and curvatures can be defined, respectively, as follows

$\begin{array}{c}{\stackrel{-}{u}}_i^{\left(\mathrm{1,3}\right)}=\frac{1}{l_q}\underset{-\frac{l_q}{2}}{\overset{\frac{l_q}{2}}{\int }}{u}_i\left(x_1^{\left(q\right)},\mp \frac{h_q}{2}\right)d{x}_1^{\left(q\right)}\\ {\stackrel{-}{u}}_i^{\left(\mathrm{2,4}\right)}=\frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}{u}_i\left(\pm \frac{l_q}{2},x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}\\ {\stackrel{-}{\theta }}_{21}^{\left(\mathrm{1,3}\right)}=\frac{1}{l_q}\underset{-\frac{l_q}{2}}{\overset{\frac{l_q}{2}}{\int }}\frac{\partial u_2\left(x_1^{\left(q\right)},\mp \frac{h_q}{2}\right)}{\partial x_1^{\left(q\right)}}d{x}_1^{\left(q\right)}\\ {\stackrel{-}{\theta }}_{12}^{\left(\mathrm{2,4}\right)}=-\frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}\frac{\partial u_1\left(\pm \frac{l_q}{2},x_2^{\left(q\right)}\right)}{\partial x_2^{\left(q\right)}}d{x}_2^{\left(q\right)}\\ {\stackrel{-}{\kappa }}_{21}^{\left(\mathrm{1,3}\right)}=\frac{1}{l_q}\underset{-\frac{l_q}{2}}{\overset{\frac{l_q}{2}}{\int }}\frac{{\partial }^2{u}_2\left(x_1^{\left(q\right)},\mp \frac{h_q}{2}\right)}{\partial x_1^{2\left(q\right)}}d{x}_1^{\left(q\right)}\\ {\stackrel{-}{\kappa }}_{12}^{\left(\mathrm{2,4}\right)}=-\frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}\frac{{\partial }^2{u}_1\left(\pm \frac{l_q}{2},x_2^{\left(q\right)}\right)}{\partial x_2^{2\left(q\right)}}d{x}_2^{\left(q\right)}\end{array}$ (2)

where the superscript indicates the subvolume face number, indexed as illustrated in Figure 2a.

Thereafter, the compatibility conditions are established in terms of the surface-averaged variables, which is motivated by the satisfaction of point-wise continuity conditions between adjacent subvolumes (Cavalcante and Pindera, 2012a). Thus, the kinematic compatibilization between the third and first faces of adjacent subvolumes is established as

$\begin{array}{c}{\left.{\stackrel{-}{u}}_i^{\left(3\right)}\right|}^{\left(q\right)}={\left.{\stackrel{-}{u}}_i^{\left(1\right)}\right|}^{\left(q+1\right)}\\ {\left.{\stackrel{-}{\theta }}_{21}^{\left(3\right)}\right|}^{\left(q\right)}={\left.{\stackrel{-}{\theta }}_{21}^{\left(1\right)}\right|}^{\left(q+1\right)}\\ {\left.{\stackrel{-}{\kappa }}_{21}^{\left(3\right)}\right|}^{\left(q\right)}={\left.{\stackrel{-}{\kappa }}_{21}^{\left(1\right)}\right|}^{\left(q+1\right)}\end{array}$ (3)

Similarly, the kinematic variables must be also compatibilized between the second and fourth faces of adjacent subvolumes. Substituting the polynomial representation of the displacement field, Eq. (1), in Eq. (2), 16 expressions are obtained for the surface-averaged displacements, rotations, and curvatures, which can be represented in matrix notation as follows

$\left[\begin{array}{c}{\stackrel{-}{\mathit{u}}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\theta }}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\kappa }}}^{\left(q\right)}\end{array}\right]={\mathit{A}}_{\left(16\times 16\right)}^{\left(q\right)}\left[\begin{array}{c}{\mathit{W}}^{\left(q\right)}\\ {\mathit{W}}_{\nabla }^{\left(q\right)}\\ {\mathit{W}}_{{\nabla }^2}^{\left(q\right)}\end{array}\right]+{\mathit{a}}_{\left(16\times 2\right)}^{\left(q\right)}{\mathit{W}}_{\left(00\right)}^{\left(q\right)}$ (4)

where ${\stackrel{-}{\mathit{u}}}^{\left(q\right)}$ is the local surface-averaged displacement vector, ${\stackrel{-}{\mathit{\theta }}}^{\left(q\right)}$ is the local surface-averaged rotation vector, ${\stackrel{-}{\mathit{\kappa }}}^{\left(q\right)}$ is the local surface-averaged curvature vector, ${\mathit{W}}^{\left(q\right)}$ is the vector containing the unknown coefficients related to the zeroth-order finite-volume theory, ${\mathit{W}}_{\nabla }^{\left(q\right)}$ is the vector formed by the unknown coefficients related to the first-order finite-volume theory, ${\mathit{W}}_{{\nabla }^2}^{\left(q\right)}$ is the vector composed by the unknown coefficients related exclusively to the second-order finite-volume theory and ${\mathit{W}}_{\left(00\right)}^{\left(q\right)}$ is the vector containing the zeroth-order unknown coefficients. ${\mathit{A}}_{\left(16\times 16\right)}^{\left(q\right)}$ and ${\mathit{a}}_{\left(16\times 2\right)}^{\left(q\right)}$ are matrices that depend only on the dimensions of the subvolume $q$.

Similarly, the surface-averaged static quantities, shown in Figure 2b, can be defined in terms of averaged tractions, first and second derivative of normal tractions acting on the faces of a generic subvolume. Thus, these surface-averaged static quantities are respectively defined as

$\begin{array}{c}{\stackrel{-}{t}}_i^{\left(\mathrm{1,3}\right)}=\mp \frac{1}{l_q}\underset{-\frac{l_q}{2}}{\overset{\frac{l_q}{2}}{\int }}{\sigma }_{2i}\left(x_1^{\left(q\right)},\mp \frac{h_q}{2}\right)d{x}_1^{\left(q\right)}\\ {\stackrel{-}{t}}_i^{\left(\mathrm{2,4}\right)}=\pm \frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}{\sigma }_{1i}\left(\pm \frac{l_q}{2},x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}\\ {\stackrel{-}{t}}_{2/1}^{\left(\mathrm{1,3}\right)}=\mp \frac{1}{l_q}\underset{-\frac{l_q}{2}}{\overset{\frac{l_q}{2}}{\int }}\frac{\partial {\sigma }_{22}\left(x_1^{\left(q\right)},\mp \frac{h_q}{2}\right)}{\partial x_1^{\left(q\right)}}d{x}_1^{\left(q\right)}\\ {\stackrel{-}{t}}_{1/2}^{\left(\mathrm{2,4}\right)}=\mp \frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}\frac{\partial {\sigma }_{11}\left(\pm \frac{l_q}{2},x_2^{\left(q\right)}\right)}{\partial x_2^{\left(q\right)}}d{x}_2^{\left(q\right)}\\ {\stackrel{-}{t}}_{2/11}^{\left(\mathrm{1,3}\right)}=\mp \frac{1}{l_q}\underset{-\frac{l_q}{2}}{\overset{\frac{l_q}{2}}{\int }}\frac{{\partial }^2{\sigma }_{22}\left(x_1^{\left(q\right)},\mp \frac{h_q}{2}\right)}{\partial x_1^{2\left(q\right)}}d{x}_1^{\left(q\right)}\\ {\stackrel{-}{t}}_{1/22}^{\left(\mathrm{2,4}\right)}=\mp \frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}\frac{{\partial }^2{\sigma }_{11}\left(\pm \frac{l_q}{2},x_2^{\left(q\right)}\right)}{\partial x_2^{2\left(q\right)}}d{x}_2^{\left(q\right)}\end{array}$ (5)

In the static analysis, the continuity conditions between the third and first faces of adjacent subvolumes are imposed as follows

$\begin{array}{c}{\left.{\stackrel{-}{t}}_i^{\left(3\right)}\right|}^{\left(q\right)}-{\left.{\stackrel{-}{t}}_i^{\left(1\right)}\right|}^{\left(q+1\right)}=0\\ {\left.{\stackrel{-}{t}}_{2/1}^{\left(3\right)}\right|}^{\left(q\right)}-{\left.{\stackrel{-}{t}}_{1/2}^{\left(1\right)}\right|}^{\left(q+1\right)}=0\\ {\left.{\stackrel{-}{t}}_{2/11}^{\left(3\right)}\right|}^{\left(q\right)}-{\left.{\stackrel{-}{t}}_{2/11}^{\left(1\right)}\right|}^{\left(q+1\right)}=0\end{array}$ (6)

These continuities must also be established between the fourth and second faces of adjacent subvolumes.

Considering linear elastic isotropic materials, where the generalized Hooke’s law is valid, 16 expressions are obtained for the surface-averaged static variables in terms of the unknown coefficients. These expressions can be arranged in matrix notation as follows

$\left[\begin{array}{c}{\stackrel{-}{\mathit{t}}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{t}}}_{\nabla }^{\left(q\right)}\\ {\stackrel{-}{\mathit{t}}}_{{\nabla }^2}^{\left(q\right)}\end{array}\right]={\mathit{B}}_{\left(16\times 16\right)}^{\left(q\right)}\left[\begin{array}{c}{\mathit{W}}^{\left(q\right)}\\ {\mathit{W}}_{\nabla }^{\left(q\right)}\\ {\mathit{W}}_{{\nabla }^2}^{\left(q\right)}\end{array}\right]$ (7)

where ${\stackrel{-}{\mathit{t}}}^{\left(q\right)}$ is the local surface-averaged traction vector, ${\stackrel{-}{\mathit{t}}}_{\nabla }^{\left(q\right)}$ is the local surface-averaged normal traction first derivative vector and ${\stackrel{-}{\mathit{t}}}_{{\nabla }^2}^{\left(q\right)}$ is the local surface-averaged normal traction second derivative vector. ${\mathit{B}}_{\left(16\times 16\right)}^{\left(q\right)}$ is a matrix that depends on the subvolume dimensions and the material elastic properties.

In the absence of body forces, the equilibrium conditions at the subvolume level are established as

$\sum _{p=1}^4{\stackrel{-}{\mathit{t}}}_{\left(p\right)}^{\left(q\right)}{L}_p^{\left(q\right)}=0$ (8)

where $L_1^{\left(q\right)}=l_q$, $L_2^{\left(q\right)}=h_q$, $L_3^{\left(q\right)}=l_q$ and $L_4^{\left(q\right)}=h_q$ are the subvolume edges lengths and ${\stackrel{-}{\mathit{t}}}_{\left(p\right)}^{\left(q\right)}$ is taken from Eq. (7) and can be expressed as

${\stackrel{-}{\mathit{t}}}_{\left(p\right)}^{\left(q\right)}={\mathit{B}}_{\left(2\times 16\right)}^{\left(q,p\right)}{\left({\mathit{A}}_{\left(16\times 16\right)}^{\left(q\right)}\right)}^{-1}\left({\stackrel{-}{\mathit{u}}}^{\left(q\right)}-{\mathit{a}}_{\left(16\times 2\right)}^{\left(q\right)}{\mathit{W}}_{\left(00\right)}^{\left(q\right)}\right)$ (9)

where ${\mathit{B}}_{\left(2\times 16\right)}^{\left(q,p\right)}$ are submatrices of selected components of the matrix ${\mathit{B}}_{\left(16\times 16\right)}^{\left(q\right)}$ related to the surface-averaged tractions acting on a face $p$ of the subvolume $q$. Replacing Eqs. (4) and (7) in Eqs. (9) and (8), the following expression is obtained

$\left(\sum _{p=1}^4{\mathit{B}}_{\left(2\times 16\right)}^{\left(q,p\right)}{L}_p^{\left(q\right)}\right){\left({\mathit{A}}_{\left(16\times 16\right)}^{\left(q\right)}\right)}^{-1}\left[\begin{array}{c}{\stackrel{-}{\mathit{u}}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\theta }}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\kappa }}}^{\left(q\right)}\end{array}\right]=\left(\sum _{p=1}^4{\mathit{B}}_{\left(2\times 16\right)}^{\left(q,p\right)}{L}_p^{\left(q\right)}\right){\left({\mathit{A}}_{\left(16\times 16\right)}^{\left(q\right)}\right)}^{-1}{\mathit{a}}_{\left(16\times 2\right)}^{\left(q\right)}{\mathit{W}}_{\left(00\right)}^{\left(q\right)}$ (10)

The vector ${\mathit{W}}_{\left(00\right)}^{\left(q\right)}$ can be obtained from Eq. (10), which is given by

${\mathit{W}}_{\left(00\right)}^{\left(q\right)}={\stackrel{-}{\mathit{a}}}_{\left(2\times 16\right)}^{\left(q\right)}\left[\begin{array}{c}{\stackrel{-}{\mathit{u}}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\theta }}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\kappa }}}^{\left(q\right)}\end{array}\right]$ (11)

where ${\stackrel{-}{\mathit{a}}}_{\left(2\times 16\right)}^{\left(q\right)}={\left[\left(\sum _{p=1}^4{\mathit{B}}_{\left(2\times 16\right)}^{\left(q,p\right)}{L}_p^{\left(q\right)}\right){\left({\mathit{A}}_{\left(16\times 16\right)}^{\left(q\right)}\right)}^{-1}{\mathit{a}}_{\left(16\times 2\right)}^{\left(q\right)}\right]}^{-1}\left(\sum _{p=1}^4{\mathit{B}}_{\left(2\times 16\right)}^{\left(q,p\right)}{L}_p^{\left(q\right)}\right){\left({\mathit{A}}_{\left(16\times 16\right)}^{\left(q\right)}\right)}^{-1}$.

Replacing Eq. (11) in Eq. (4), the following expression can be obtained:

$\left[\begin{array}{c}{\mathit{W}}^{\left(q\right)}\\ {\mathit{W}}_{\nabla }^{\left(q\right)}\\ {\mathit{W}}_{{\nabla }^2}^{\left(q\right)}\end{array}\right]={\stackrel{-}{\mathit{A}}}_{\left(16\times 16\right)}^{\left(q\right)}\left[\begin{array}{c}{\stackrel{-}{\mathit{u}}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\theta }}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\kappa }}}^{\left(q\right)}\end{array}\right]$ (12)

where ${\stackrel{-}{\mathit{A}}}_{\left(16\times 16\right)}^{\left(q\right)}={\left({\mathit{A}}_{\left(16\times 16\right)}^{\left(q\right)}\right)}^{-1}-{\left({\mathit{A}}_{\left(16\times 16\right)}^{\left(q\right)}\right)}^{-1}{\mathit{a}}_{\left(16\times 2\right)}^{\left(q\right)}{\stackrel{-}{\mathit{a}}}_{\left(2\times 16\right)}^{\left(q\right)}$. Thus, the local system of equations for a generic subvolume can be obtained by replacing Eq. (12) in Eq. (7), which gives

$\left[\begin{array}{c}{\stackrel{-}{\mathit{t}}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{t}}}_{\nabla }^{\left(q\right)}\\ {\stackrel{-}{\mathit{t}}}_{{\nabla }^2}^{\left(q\right)}\end{array}\right]={\mathit{K}}_{\left(16\times 16\right)}^{\left(q\right)}\left[\begin{array}{c}{\stackrel{-}{\mathit{u}}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\theta }}}^{\left(q\right)}\\ {\stackrel{-}{\mathit{\kappa }}}^{\left(q\right)}\end{array}\right]$ (13)

where ${\mathit{K}}_{\left(16\times 16\right)}^{\left(q\right)}={\mathit{B}}_{\left(16\times 16\right)}^{\left(q\right)}{\stackrel{-}{\mathit{A}}}_{\left(16\times 16\right)}^{\left(q\right)}$ is the local stiffness matrix.

For the global stiffness matrix assemblage, it is considered the individual contribution of each subvolume on the discretized structure. Therefore, the global system of equations can be defined as

${\mathit{T}}_{\left(ndof\times 1\right)}={\mathit{K}}_{\left(ndof\times ndof\right)}{\mathit{U}}_{\left(ndof\times 1\right)}$ (14)

where $ndof$ is the total number of degrees of freedom, ${\mathit{T}}_{\left(ndof\times 1\right)}$ and ${\mathit{U}}_{\left(ndof\times 1\right)}$ are the global surface-averaged static and kinematic vectors, respectively, and ${\mathit{K}}_{\left(ndof\times ndof\right)}$ is the global stiffness matrix evaluated by

${\mathit{K}}_{\left(ndof\times Ndf\right)}=\sum _{q=1}^{N_q}\left[{\left({\mathit{L}}_{\left(16\times ndof\right)}^{\left(q\right)}\right)}^T{\mathit{K}}_{\left(16\times 16\right)}^{\left(q\right)}{\mathit{L}}_{\left(16\times ndof\right)}^{\left(q\right)}\right]$ (15)

where ${\mathit{L}}_{\left(16\times ndof\right)}^{\left(q\right)}$ is the kinematic and static incidence matrix.

The previous theoretical development corresponds to the formulation of the second-order version of the generalized finite-volume theory for continuum elastic structures. For the lower order versions of the generalized finite-volume theory, the framework can be obtained by uncoupling curvatures, in the case of the first-order version, and curvatures and rotations, in the case of the zeroth-order version. The vectors composed by the unknown coefficients must also be uncoupled following the corresponding version of the generalized finite-volume theory.

3 TOPOLOGY OPTIMIZATION PROBLEM

In general, the topology optimization problem is formulated as an algorithm that seeks to find the best material distribution inside a reference domain. Since Bendsøe and Kikuchi (1988), a significant part of the advances in topology optimization has been obtained through methodologies based on the compliance minimization problem, whose concepts are well-established (Collet et al., 2017). Some examples of applications using this type of optimization problem are presented in Rozvany (2009), Lopes et al. (2015), Shobeiri (2016), and Wang et al. (2017). Here, it is used the topology optimization problem for compliance minimization, where the design domain is assumed to be rectangular and discretized in rectangular elements or subvolumes.

The topology optimization problem based on the power-law approach applied in the context of the finite element method, where the objective is to minimize the compliance structural function under volume constraint, can be described as

$\left\{\begin{array}{c}\min c\left(\rho \right)=\sum _{e=1}^N{\left({\rho }_e\right)}^p{\mathit{d}}_e^T{\mathit{K}}_e^0{\mathit{d}}_e\\ \text{s}\text{u}\text{b}\text{j}\text{e}\text{c}\text{t}\text{t}\text{o}\text{:}\\ \frac{V\left(\mathit{\rho }\right)}{\stackrel{-}{V}}=f\\ 0<{\rho }_{min}\le {\rho }_e\le 1\end{array}\right.$ (16)

where $V\left(\mathit{\rho }\right)$ and $\stackrel{-}{V}$ are the material and the reference domain volumes, respectively, ${\mathit{d}}_e$ is the local displacement vector, ${\mathit{K}}_e^0$ is the local stiffness matrix for a unitary relative density, $\mathit{\rho }$ is the relative density vector, $p$ is the penalty factor, $f$ is the prescribed volume fraction, ${\rho }_{min}$ is the minimum relative density, ${\rho }_e$ is the relative density associated with each element and $N$ is the total number of elements.

The optimization problem presented in Eq. (16) is solved using the classical approach denoted by optimality criteria (OC) method. Therefore, following the procedure suggested by Bendsøe and Sigmund (2003), a heuristic update for the design variables is established as

${\rho }_e^{k+1}=\left\{\begin{array}{c}\max \left({\rho }_{min},{\rho }_e-m\right),\\ if{\rho }_e^k{B}_e^{\eta }\le \max \left({\rho }_{min},{\rho }_e-m\right),\\ {\rho }_e^k{B}_e^{\eta },\\ if\max \left({\rho }_{min},{\rho }_e-m\right)<{\rho }_e^k{B}_e^{\eta }\le \min \left(1,{\rho }_e+m\right),\\ \min \left(1,{\rho }_e+m\right),\\ if\min \left(1,{\rho }_e+m\right)\le {\rho }_e^k{B}_e^{\eta }\end{array}\right.$ (17)

where $k$ is the current iteration, $m$ is the move-limit, $\eta$ is the damping factor and $B_e$ is given by

$B_e=-\frac{\frac{\partial c}{\partial {\rho }_e}}{\lambda \frac{\partial f}{\partial {\rho }_e}}$ (18)

where $\lambda$ is the Lagrangian multiplier for the constrained volume, which is determined by a bisection method.

The damping factor can be used to regularize possible oscillations during the optimization, mainly when no filtering techniques are employed. The parameter $\eta$ is directly related to the method performance, once this affects the speed variation of $B_e^{\eta }$ (Montes, 2016). A high value for $\eta$ can accelerate the optimization convergence process, which may cause oscillations in the displacement field for the low-density regions (Ma et al., 1993). Also, the adoption of minor values of $\eta$ can prevent divergence in the topology optimization algorithm; however, this results in small changes in the design variables, which leads to a slower convergence process (Ma et al., 1993). The value of $\eta$ that provides the faster convergence for the overall process is 1/2, so it is recommended to maintain the damping factor as close as possible of this value.

3.1 Mesh-independency filter

To avoid the occurrence of mesh dependency, it is suggested the modification of the elements’ sensitivities by the following expression:

$\frac{\partial c}{\partial {\rho }_e}=\frac{1}{{\rho }_e\sum _{f=1}^N{\widehat{H}}_f}\sum _{f=1}^N{\widehat{H}}_f{\rho }_f\frac{\partial c}{\partial {\rho }_f}$ (19)

where ${\widehat{H}}_f$ is the convolution operator (weighting function) given as

${\widehat{H}}_f=R-\mathrm{dist}\left(e,f\right)\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{dist}\left(e,f\right)\le R\mathrm{a}\mathrm{n}\mathrm{d}{\widehat{H}}_f=0\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},$ (20)

where $\mathrm{dist}\left(e,f\right)$ is the distance between the element center of $e$ and the element center of $f$ (Sigmund, 2001). To consider the contribution only of the neighbor elements (with shared nodes), it is adopted a filter radius of $R=1.01\sqrt{{\left(l_e\right)}^2+{\left(h_e\right)}^2}$.

4 COMPLIANCE FUNCTION FOR THE GENERALIZED FINITE-VOLUME THEORY

In general, the total strain energy of a deforming material and the work done by external loadings are equivalent to a conservative internal force system in a quasi-static analysis. Therefore, in structural analysis, this principle is mostly observed on energy-based numerical methods, as the finite element method. However, in the finite-volume theory, this feature is observed only for the zeroth-order version of the generalized finite-volume theory, since the local equilibrium are established only in terms of the surface-averaged tractions. As a result, for the first and second-order versions, the equivalence between the total strain energy and the work done by external forces is observed only when a sufficiently refined mesh is employed. One of the main objectives of this contribution is to define whether the total strain energy or the total work done by external forces produce the best results for the proposed optimization problem.

The compliance function can be defined as twice the total strain energy produced by a displacement field $\mathit{u}$; thus, this function can be expressed as

$c\left(\mathit{u},\mathit{\rho }\right)=2U\left(\mathit{u},\mathit{\rho }\right)=\underset{\Omega }{\overset{}{\iiint }}2\stackrel{-}{U}\left(\mathit{u},\mathit{\rho }\right)d\Omega =\underset{\Omega }{\overset{}{\iiint }}{\sigma }_{ij}\left(\mathit{u},\mathit{\rho }\right){\epsilon }_{ij}\left(\mathit{u}\right)d\Omega =\underset{\Omega }{\overset{}{\iiint }}{C}_{ijkl}\left(\mathit{\rho }\right){\epsilon }_{kl}\left(\mathit{u}\right){\epsilon }_{ij}\left(\mathit{u}\right)d\Omega$ (21)

where $U\left(\mathit{u},\mathit{\rho }\right)$ is the total strain energy, $\stackrel{-}{U}\left(\mathit{u},\mathit{\rho }\right)$ is the specific strain energy, ${\sigma }_{ij}\left(\mathit{u},\mathit{\rho }\right)$ is the stress tensor, ${\epsilon }_{ij}\left(\mathit{u}\right)$ is the strain tensor, $C_{ijkl}\left(\mathit{\rho }\right)$ is the stiffness tensor and $\Omega$ is the reference domain.

In the absence of body forces and considering $u_i=0$ in $S_u$, the work done by external loadings can be defined as

$W=\frac{1}{2}{\iint }_{S_{\sigma }}^{}{t}_i{u}_{i}dS$ (22)

where $t_i$ are the traction components, $u_i$ are the displacement components, $S_{\sigma }$ is the external surface subjected to external loadings, $S_u$ is the external surface with predicted displacements and $S=S_u\cup S_{\sigma }$.

Applying Cauchy’s law and the divergence theorem to Eq. (22), it follows

$2W={\iiint }_{\Omega }^{}\frac{\partial {\sigma }_{ji}}{\partial x_j}{u}_{i}d\Omega +{\iiint }_{\Omega }^{}{\sigma }_{ij}{\epsilon }_{ij}d\Omega +{\iiint }_{\Omega }^{}{\sigma }_{ji}{\omega }_{ij}d\Omega$ (23)

where ${\omega }_{ij}$ is the asymmetric rotation tensor. Considering the symmetry of ${\sigma }_{ij}$, Eq. (23) can be written as

$2W={\iiint }_{\Omega }^{}\frac{\partial {\sigma }_{ij}}{\partial x_j}{u}_{i}d\Omega +2U\left(\mathit{u}\right)$ (24)

As a result, the equivalence between the total work done by external forces and the total strain energy is observed only if $\partial {\sigma }_{ij}/\partial x_j=0$. However, this is valid only for the zeroth-order finite-volume theory. Consequently, for the higher-order versions of the generalized finite-volume theory, it is verified the need to investigate the different aspects that involve mechanical energy evaluation.

4.1 Total strain energy for the generalized finite-volume theory

For linearly elastic materials, the total strain energy can be defined as

$U=\underset{\Omega }{\overset{}{\iiint }}\frac{1}{2}{\mathit{\sigma }}^T\mathit{\epsilon }d\Omega ={\iiint }_{\Omega }^{}\frac{1}{2}{\mathit{\epsilon }}^T\mathit{C}\mathit{\epsilon }d\Omega$ (25)

where $\mathit{\sigma }$ is the stress tensor, $\mathit{\epsilon }$ is the strain tensor and $\mathit{C}$ is the material stiffness tensor. Considering the displacement approximation presented in Eq. (1), the strain tensor of a subvolume can be described as

${\mathit{\epsilon }}^{\left(q\right)}\left(x_1^{\left(q\right)},x_2^{\left(q\right)}\right)={\mathit{E}}_0^{\left(q\right)}\left(x_1^{\left(q\right)},x_2^{\left(q\right)}\right){\mathit{W}}^{\left(q\right)}+{\mathit{E}}_1^{\left(q\right)}\left(x_1^{\left(q\right)},x_2^{\left(q\right)}\right){\mathit{W}}_{\nabla }^{\left(q\right)}+{\mathit{E}}_2^{\left(q\right)}\left(x_1^{\left(q\right)},x_2^{\left(q\right)}\right){\mathit{W}}_{{\nabla }^2}^{\left(q\right)}$ (26)

where ${\mathit{E}}_m^{\left(q\right)}\left(x_1^{\left(q\right)},x_2^{\left(q\right)}\right)$ are matrices written as functions of the local coordinates and relate the strain tensor components to the unknown coefficients of the displacement field.

The strain energy at the subvolume level can be evaluated as

$U^{\left(q\right)}={\int }_{-\frac{h_q}{2}}^{\frac{h_q}{2}}{\int }_{-\frac{l_q}{2}}^{\frac{l_q}{2}}\frac{1}{2}{{\mathit{\epsilon }}^{\left(q\right)}}^T{\mathit{C}}^{\left(q\right)}{\mathit{\epsilon }}^{\left(q\right)}d{x}_1^{\left(q\right)}d{x}_2^{\left(q\right)}$ (27)

Substituting Eq. (26) in Eq. (27) and integrating in terms of the local coordinates, the local strain energy, considering the second-order version of the generalized finite-volume theory, can be defined as

$U^{\left(q\right)}=\frac{1}{2}\left({\mathit{W}}^{\left(q\right)}{\mathit{D}}_{00}^{\left(q\right)}{\mathit{W}}^{\left(q\right)}+{\mathit{W}}^{\left(q\right)}{\mathit{D}}_{01}^{\left(q\right)}{\mathit{W}}_{\nabla }^{\left(q\right)}+{\mathit{W}}_{\nabla }^{\left(q\right)}{\mathit{D}}_{10}^{\left(q\right)}{\mathit{W}}^{\left(q\right)}+{\mathit{W}}_{\nabla }^{\left(q\right)}{\mathit{D}}_{11}^{\left(q\right)}{\mathit{W}}_{\nabla }^{\left(q\right)}+{\mathit{W}}_{{\nabla }^2}^{\left(q\right)}{\mathit{D}}_{22}^{\left(q\right)}{\mathit{W}}_{{\nabla }^2}^{\left(q\right)}\right)$ (28)

where

${\mathit{D}}_{mn}^{\left(q\right)}={\int }_{-\frac{h_q}{2}}^{\frac{h_q}{2}}{\int }_{-\frac{l_q}{2}}^{\frac{l_q}{2}}\frac{1}{2}{{\mathit{E}}_m^{\left(q\right)}}^T{\mathit{C}}^{\left(q\right)}{\mathit{E}}_n^{\left(q\right)}d{x}_1^{\left(q\right)}d{x}_2^{\left(q\right)}$ (29)

Besides, ${\mathit{D}}_{02}^{\left(q\right)}$, ${\mathit{D}}_{12}^{\left(q\right)}$, ${\mathit{D}}_{20}^{\left(q\right)}$ and ${\mathit{D}}_{21}^{\left(q\right)}$ are found to be null matrices. The total strain energy is obtained by considering the individual contribution of each subvolume; thus, it can be written as

$U=\underset{\Omega }{\overset{}{\iiint }}\stackrel{-}{U}d\Omega =\sum _{q=1}^{N_q}{U}^{\left(q\right)}$ (30)

The strain energy for the lower order versions of the generalized finite-volume theory can be obtained by uncoupling the vector ${\mathit{W}}_{{\nabla }^2}^{\left(q\right)}$, in the case of the first-order version, and ${\mathit{W}}_{{\nabla }^2}^{\left(q\right)}$ and ${\mathit{W}}_{\nabla }^{\left(q\right)}$, in the case of the zeroth-order version.

4.2 Total work done by external loadings based on the generalized finite-volume theory

The horizontal displacement at the subvolume vertical faces can be expressed by three Legendre polynomials as follows

$u_1^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)={\mu }_{1\left(0\right)}^{\left(\mathrm{2,4}\right)}+x_2^{\left(q\right)}{\mu }_{1\left(1\right)}^{\left(\mathrm{2,4}\right)}+\frac{1}{2}\left(3{x_2^{\left(q\right)}}^2-\frac{h_q^2}{4}\right){\mu }_{1\left(2\right)}^{\left(\mathrm{2,4}\right)}$ (31)

where ${\mu }_{1\left(m\right)}^{\left(\mathrm{2,4}\right)}$ are unknown coefficients of the horizontal displacement at the subvolume vertical faces. From Eq. (31), the surface-averaged kinematic quantities can be evaluated as

$\begin{array}{c}{\stackrel{-}{u}}_1^{\left(\mathrm{2,4}\right)}=\frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}{u}_1^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}={\mu }_{1\left(0\right)}^{\left(\mathrm{2,4}\right)}\\ {\stackrel{-}{\theta }}_{12}^{\left(\mathrm{2,4}\right)}=-\frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}\frac{d{u}_1^{\left(\mathrm{2,4}\right)}}{d{x}_2^{\left(q\right)}}\left(x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}=-{\mu }_{1\left(1\right)}^{\left(\mathrm{2,4}\right)}\\ {\stackrel{-}{\kappa }}_{12}^{\left(\mathrm{2,4}\right)}=-\frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}\frac{d^2{u}_1^{\left(\mathrm{2,4}\right)}}{d{x}_2^{2\left(q\right)}}\left(x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}=-3{\mu }_{1\left(2\right)}^{\left(\mathrm{2,4}\right)}\end{array}$ (32)

Similarly, for linearly elastic materials, the normal traction acting on the subvolume vertical faces can be expressed by three Legendre polynomials as follows

$t_1^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)={\tau }_{1\left(0\right)}^{\left(\mathrm{2,4}\right)}+x_2^{\left(q\right)}{\tau }_{1\left(1\right)}^{\left(\mathrm{2,4}\right)}+\frac{1}{2}\left(3{x_2^{\left(q\right)}}^2-\frac{h_q^2}{4}\right){\tau }_{1\left(2\right)}^{\left(\mathrm{2,4}\right)}$ (33)

where ${\tau }_{1\left(m\right)}^{\left(\mathrm{2,4}\right)}$ are unknown coefficients of the normal tractions acting on the vertical faces. From Eq. (33), the surface-averaged static quantities can be evaluated by

$\begin{array}{c}{\stackrel{-}{t}}_1^{\left(\mathrm{2,4}\right)}=\frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}{t}_1^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}={\tau }_{1\left(0\right)}^{\left(\mathrm{2,4}\right)}\\ {\stackrel{-}{t}}_{1/2}^{\left(\mathrm{2,4}\right)}=-\frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}\frac{d{t}_1^{\left(\mathrm{2,4}\right)}}{d{x}_2^{\left(q\right)}}\left(x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}=-{\tau }_{1\left(1\right)}^{\left(\mathrm{2,4}\right)}\\ {\stackrel{-}{t}}_{1/22}^{\left(\mathrm{2,4}\right)}=-\frac{1}{h_q}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}\frac{d^2{t}_1^{\left(\mathrm{2,4}\right)}}{d{x}_2^{2\left(q\right)}}\left(x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}=-3{\tau }_{1\left(2\right)}^{\left(\mathrm{2,4}\right)}\end{array}$ (34)

Using Eqs. (33) and (34), the resultant static quantities associated to the subvolume vertical faces are determined as follows

$\begin{array}{c}{R}_1^{\left(\mathrm{2,4}\right)}=\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}{t}_1^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}={\stackrel{-}{t}}_1^{\left(\mathrm{2,4}\right)}{h}_q\\ R_2^{\left(\mathrm{2,4}\right)}=\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}{t}_2^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}={\stackrel{-}{t}}_2^{\left(\mathrm{2,4}\right)}{h}_q\\ M_R^{\left(\mathrm{2,4}\right)}=-\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}{t}_1^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)x_2^{\left(q\right)}d{x}_2^{\left(q\right)}={\stackrel{-}{t}}_{1/2}^{\left(\mathrm{2,4}\right)}\frac{h_q^3}{12}\\ S_R^{\left(\mathrm{2,4}\right)}=\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}{t}_1^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)\frac{1}{6}\left(3{x_2^{\left(q\right)}}^2-\frac{h_q^2}{4}\right)d{x}_2^{\left(q\right)}={\stackrel{-}{t}}_{1/22}^{\left(\mathrm{2,4}\right)}\frac{h_q^5}{720}\end{array}$ (35)

where $R_i^{\left(\mathrm{2,4}\right)}$ are the resultant forces, $M_R^{\left(\mathrm{2,4}\right)}$ is the resultant bending moment and $S_R^{\left(\mathrm{2,4}\right)}$ is the second-order bending moment acting on the subvolume vertical faces. Using Eqs. (31), (32), (33) and (34), the total work done on a vertical face can be evaluated by

$W^{\left(\mathrm{2,4}\right)}=\frac{1}{2}\underset{-\frac{h_q}{2}}{\overset{\frac{h_q}{2}}{\int }}{t}_i^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)u_i^{\left(\mathrm{2,4}\right)}\left(x_2^{\left(q\right)}\right)d{x}_2^{\left(q\right)}=\frac{1}{2}\left({\stackrel{-}{u}}_1^{\left(\mathrm{2,4}\right)}{\stackrel{-}{t}}_1^{\left(\mathrm{2,4}\right)}{h}_q+{\stackrel{-}{u}}_2^{\left(\mathrm{2,4}\right)}{\stackrel{-}{t}}_2^{\left(\mathrm{2,4}\right)}{h}_q+{\stackrel{-}{\theta }}_{12}^{\left(\mathrm{2,4}\right)}{\stackrel{-}{t}}_{1/2}^{\left(\mathrm{2,4}\right)}\frac{h_q^3}{12}+{\stackrel{-}{\kappa }}_{12}^{\left(\mathrm{2,4}\right)}{\stackrel{-}{t}}_{1/22}^{\left(\mathrm{2,4}\right)}\frac{h_q^5}{720}\right)$ (36)

Extending the results presented on Eq. (36) to all subvolume faces and using Eq. (35), the work done by external loadings for the second-order version of the generalized finite-volume theory can be determined as follows

$W^{\left(q\right)}=\frac{1}{2}\left({\mathit{R}}^{\left(q\right)}{\stackrel{-}{\mathit{u}}}^{\left(q\right)}+{\mathit{M}}_{\mathit{R}}^{\left(q\right)}{\stackrel{-}{\mathit{\theta }}}^{\left(q\right)}+{\mathit{S}}_{\mathit{R}}^{\left(q\right)}{\stackrel{-}{\mathit{\kappa }}}^{\left(q\right)}\right)$ (37)

where ${\mathit{R}}^{\left(q\right)}$ is the local resultant force vector, ${\mathit{M}}_{\mathit{R}}^{\left(q\right)}$ is the local resultant bending moment vector and ${\mathit{S}}_{\mathit{R}}^{\left(q\right)}$ is the local resultant second-order bending moment vector acting on the subvolume faces. The total work done by external loadings is given by the sum of each subvolume contribution, which can be expressed as

$W=\sum _{q=1}^{N_q}{W}^{\left(q\right)}$ (38)

For the lower-order versions of the generalized finite-volume theory, the expressions for the work done by external loadings can be obtained by uncoupling the surface-averaged curvatures, in the case of the first-order version, and rotations and curvatures, in the case of the zeroth-order version.

5 NUMERICAL RESULTS

On this contribution, two different examples are analyzed, employing the three versions of the generalized finite-volume theory and Q4 and Q8 elements of the finite element method. The studied examples are a cantilever beam subjected to a concentrated load and a Messerschmitt-Bölkow-Blohm (MBB) beam. To confirm the checkerboard-free property of the approaches based on the finite-volume theory, the examples are firstly analyzed without employing filtering or image processing techniques. After that, the same examples are analyzed employing a filter that modifies the elements or subvolumes sensitivities, as presented in Eq. (19), for mesh-independency.

Some numerical aspects are investigated during the analysis, such as the number of iterations, processing time, convergence, and relative compliance. The continued penalization scheme is adopted, where the penalty factor increases gradually ($∆p=0.5$) from 1 to 4, as suggested by Talischi et al. (2012). As a convergence criterion, the tolerance for the maximum change between relative densities of successive steps is assumed to be 1%. In the absence of filtering techniques, each approach's damping factor is adjusted to avoid any divergence during the optimization process. The damping factor is set up as close as possible to 1/2, since no oscillation in the displacement field is verified when the algorithm is performed. The adopted damping factor for each simulation is shown in the following Tables and was obtained by varying increments of 0.1 as follows: 1/2, 1/2.1, 1/2.2, …, until its findings convergence in the optimization process.

5.1 Cantilever beam

A classical problem in the topology optimization of bidimensional structures is the cantilever beam, whose analysis domain and boundary conditions are presented in Figure 3. The proposed optimization problem consists of minimizing the structural compliance function, defined as twice the total strain energy, with a volume constraint of 40% of the total volume. The computational environment, in terms of programming language and machine, can be described as MatLab R2016a (64-bits)/Intel® CoreTM i7 CPU 2.93 GHz/16.0 GB RAM/64-bits. Consistent units are employed for the physical and geometrical parameters.

Figure 3:
Cantilever beam.

Figure 4 shows the optimum topologies obtained for each studied mesh size and employing the zeroth, first and second-order finite-volume theory (FVT0^th, FVT^1st, and FVT^2nd, respectively) and the Q4 and Q8 elements. Additionally, Table 1 presents the investigated numerical aspects. From Figure 4, the optimum topologies obtained employing the finite-volume theory approaches have shown to be checkerboard-free. However, the approaches based on the finite element method have generated optimum topologies with the checkerboard pattern issue. The checkerboard pattern problem in optimum topologies is directly related to the displacement assumptions of the finite element method, leading to structures artificially rigid (Díaz and Sigmund, 1995). On the other hand, the satisfaction of equilibrium equations and continuity conditions through the faces of adjacent subvolumes guarantees the checkerboard-free property for the different versions of the finite-volume theory, even when no filtering technique is employed.

Table 1 presents the total number of iterations, the processing time, the number of degrees of freedom (NDOF), and the adopted damping factor, set up to avoid divergence in the optimization process. In general, the number of iterations has varied from one approach to another, presenting higher values when the first-order finite-volume theory and Q4 approaches are employed, and the lowest value was obtained for the second-order finite-volume theory followed by the Q8 and FVT0th approaches. The zeroth-order finite-volume theory has been approximately 1.08 times slower than the Q4 approach for the finest mesh in terms of computational cost. The Q8 approach has presented the highest computational cost: 1.20 times slower than the first-order finite-volume theory and 1.04 times slower than the second-order finite-volume theory, for the finest mesh. The number of degrees of freedom explains these differences in the computational cost partially since it defines the size of the global system of equations.

Figure 4:
Optimum topologies for the cantilever beam analysis evaluating the compliance using the strain energy (No filtering).

Table 1:
Convergence analysis for the cantilever beam evaluating the compliance using the strain energy (No filtering).

Analysis	Mesh	NDOF	Number of Iterations	Processing Time (s)	Damping factor
Q4	42x21	1892	270	36.95	1/2.5
	82x41	6972	360	635.26	1/2.5
	162x81	26732	710	19495.87	1/2.6
Q8	42x21	5546	266	157.18	1/2.3
	82x41	20666	397	3138.15	1/2.4
	162x81	79706	485	54375.84	1/2.6
FVT0th	42x21	3654	252	51.10	1/2.6
	82x41	13694	306	764.33	1/2.6
	162x81	52974	488	21133.43	1/2.6
FVT1st	42x21	5481	394	141.96	1/2.8
	82x41	20541	376	1500.75	1/2.8
	162x81	79461	599	45207.53	1/3.6
FVT2nd	42x21	7308	287	85.42	1/3.0
	82x41	27388	326	2933.66	1/3.0
	162x81	105948	453	52327.55	1/3.2

Figure 5:
Optimum topologies for the cantilever beam analysis evaluating the compliance using the work done by external loadings (No filtering).

Similarly, the proposed optimization problem can also be solved by defining the structural compliance as twice the work done by external loadings. The obtained optimum topologies are shown in Figure 5 for the first and second-order finite-volume theory since the total strain energy and the work done by external forces are equivalent to the approaches based on the finite element method and the zeroth-order finite-volume theory. The numerical aspects investigated for convergence analysis can be found in Table 2.

Table 2:
Convergence analysis for the cantilever beam evaluating the compliance using the work done by external loadings (No filtering).

Analysis	Mesh	NDOF	Number of Iterations	Processing Time (s)	Damping factor
FVT1st	42x21	5481	374	64.51	1/6.2
	82x41	20541	401	1297.19	1/8.0
	162x81	79461	696	49691.47	1/10.0
FVT2nd	42x21	7308	287	84.94	1/3.0
	82x41	27388	514	3414.41	1/8.3
	162x81	105948	641	72673.70	1/10.2

When the compliance function is estimated using the external work done by external loadings, the optimum topologies tend to show more bars and length scale issues, as illustrated on the optimum topologies presented in Figure 5. The damping factors for these approaches have shown to be much lower when compared to the same approaches employing the strain energy, which turns the convergence process slower and increases the computational cost. The number of iterations tends to be higher, making the approaches employing the external work done by external loadings more computational costly, as shown in Table 2. In fact, for the current example, the objective function is better estimated when the compliance function is defined as twice the total strain energy.

As shown in Figures 4 and 5, although the checkerboard pattern issue can be overcome by the topology optimization approach based on the finite-volume theory, the mesh dependence between successive meshes persists. Therefore, in this contribution, the mesh-independency filter, presented in section 3.1, is employed to avoid mesh dependence, in the approaches based on the finite-volume theory, and checkerboard pattern and mesh dependence, in the approaches based on the finite element method. The optimum topologies for the same problem employing the sensitivity filtering are presented in Figure 6. In this case, the compliance function is evaluated as twice the structural strain energy, and the damping factor is adjusted as 1/2 for all investigated approaches.

The optimum topologies presented in Figure 6 are checkerboard-free, and the mesh dependence is better controlled in this scenario. There are some differences between the optimum topologies obtained by the finite-volume theory and the approaches based on the finite element method. The optimum topologies obtained by the zeroth-order finite-volume theory have dramatically reduced the mesh dependence between successive meshes, providing topologies with fewer bars and reducing the length scale issue, which are desirable features for manufacturing.

Figure 6:
Optimum topologies for the cantilever beam analysis evaluating the compliance using the strain energy (filtering).

The most critical topologies are obtained for the Q4 element approach; in this case, the optimum topologies present slender bars with length scale issues, undesirable features for manufacturing. The approaches based on the Q8 element, first and second-order finite-volume theories have presented similar optimum topologies with more bars when compared to the zeroth-order finite-volume theory approach, and fewer bars and length scale issues when compared to the Q4 element approach. In general, the optimum topologies obtained by the zeroth-order finite-volume theory approach are well behaved and more indicated for the design of optimum structures.

Table 3 presents the results obtained for the overall convergence of the different topology optimization approaches, considering the application of a mesh independent filter that regularizes the element or subvolume sensitivities. In general, the number of iterations has changed from one approach to another, where the minimum values are observed for the second-order finite-volume theory and the Q4 approaches. In terms of computational cost, the Q8 approach has presented the highest processing time, while the Q4 approach has presented the lowest computational cost. The approach based on the zeroth-order finite-volume theory is 1.8 times slower than the same approach based on the Q4 element, for the finest mesh. The number of degrees of freedom explains the computational efficiency of the Q4 approach partially since it defines the size of the global system of equations.

In Table 3, it is presented a numerical parameter, denoted as relative compliance, permitting the comparison between the different approaches in terms of the lowest compliant structure. This value is obtained by recalculating each optimum topology's structural compliance employing a Q8 element; after that, this result is divided by the compliance obtained for the Q8 optimum topologies for the same mesh sizes. The relative compliance values show that the stiffest structures are obtained when the zeroth-order finite-volume theory is employed, especially for the finest mesh. The approaches based on the Q4 element and the second-order finite-volume theory have shown to be more flexible and less indicated to optimize the analyzed structure.

Table 3:
Convergence analysis for the cantilever beam evaluating the compliance using the strain energy (filtering).

Analysis	Mesh	NDOF	Number of Iterations	Processing Time (s)	Relative compliance
Q4	70x35	5112	230	306.37	1.00053
	150x75	22952	293	5270.43	1.00056
	230x115	53592	485	52254.96	1.00364
Q8	70x35	15122	153	645.97	1.00000
	150x75	68402	320	25743.90	1.00000
	230x115	160082	543	248635.27	1.00000
FVT0th	70x35	10010	156	195.15	1.01019
	150x75	45450	329	8341.55	0.99068
	230x115	106490	630	94014.33	0.99536
FVT1st	70x35	15015	156	403.17	0.99833
	150x75	68175	285	15400.58	0.99952
	230x115	159735	526	156290.92	0.99731
FVT2nd	70x35	20020	152	667.21	1.00309
	150x75	90900	254	20218.72	1.00419
	230x115	212980	441	199348.83	1.00246

5.2 Messerschmitt-Bölkow-Blohm (MBB) beam

Another classical problem in topology optimization is known as Messerschmitt-Bölkow-Blohm (MBB) beam, whose analysis domain and boundary conditions are shown in Figure 7. The optimization problem consists of finding the stiffest structure with a given volume fraction of 50%. Taking advantage of the structure symmetry, only half of the structure is analyzed, employing boundary conditions that reflect this symmetry. Additionally, in the model conception, consistent units for the physical and geometric parameters are employed. The computational environment for this example, in terms of programming language and machine, can be described as MatLab R2018a (64-bits)/Intel® CoreTM i7-8550U CPU @ 1.80 GHz 1.99 GHz/16.0 GB RAM/64-bits.

Figure 7:
Messerschmitt-Bölkow-Blohm (MBB) beam.

Figure 8:
Optimum topologies for the MBB beam analysis evaluating the compliance using the strain energy (No filtering).

Figure 8 shows the optimum topologies for the approaches based on the finite-volume theory and Q4 and Q8 elements in the absence of filtering or image processing techniques, where the compliance function is evaluated as twice the total strain energy. As presented in the previous example, the checkerboard pattern is an issue for the finite element-based approaches, mainly when the Q4 element is employed. Even for the finest mesh of the Q8 element, the checkerboard pattern appears very locally, which is an issue for manufacturing. On the other hand, for the three versions of the generalized finite-volume theory, it is not verified the presence of any checkerboard regions. The damping factor for each simulation was adjusted to avoid any divergence during the optimization process. The adopted damping factors are shown in Table 4.

Table 4:
Convergence analysis for the MBB beam evaluating the compliance using the strain energy (No filtering).

Analysis	Mesh	NDOF	Number of Iterations	Processing Time (s)	Damping factor
Q4	45x15	1472	248	9.96	1/2.5
	90x30	5642	717	429.20	1/2.5
	180x60	22082	763	12306.13	1/2.6
Q8	45x15	4292	322	42.81	1/2.5
	90x30	16682	422	1889.06	1/2.5
	180x60	65762	545	41312.91	1/2.5
FVT0th	45x15	2820	253	13.44	1/2.6
	90x30	11040	422	315.81	1/2.6
	180x60	43680	584	14110.88	1/2.6
FVT1st	45x15	4230	401	47.84	1/3.2
	90x30	16560	458	1116.78	1/3.0
	180x60	65520	665	32551.37	1/2.7
FVT2nd	45x15	5640	285	39.70	1/3.0
	90x30	22080	362	1220.22	1/3.0
	180x60	87360	606	43860.14	1/3.4

Table 4 presents the obtained results for the overall convergence analysis employing the different versions of the generalized finite-volume theory, and Q4 and Q8 elements of the finite element method, where the structural compliance is defined as twice the total strain energy. The number of iterations has varied from one approach to another, mainly when the Q4 element is employed. The adopted damping factor explains the high number of iterations partially for the finite-volume theory approaches since it provides a slow convergence for the optimization process. In terms of computational cost, the second-order finite-volume theory has presented the highest processing time, while the Q4 approach has presented the lowest computational cost. As a result, the zeroth-order finite-volume theory is 1.15 slower than the Q4 approach and 2.93 times faster than the Q8 approach. The second-order finite-volume theory is 1.06 times slower than the Q8 approach and 3.56 slower than the Q4 approach. The first-order finite-volume theory is 2.65 times slower than the Q4 approach and 1.27 times faster than the Q8 approach.

Figure 9:
Optimum topologies for the MBB beam analysis evaluating the compliance using the work done by external loadings (No filtering).

Table 5:
Convergence analysis for the MBB beam evaluating the compliance using the work done by external loadings (No filtering).

Analysis	Mesh	NDOF	Number of Iterations	Processing Time (s)	Damping factor
FVT1st	45x15	4230	368	54.40	1/9.8
	90x30	16560	643	1393.07	1/8.9
	180x60	65520	795	41768.21	1/11.8
FVT2nd	45x15	5640	404	102.53	1/9.6
	90x30	22080	571	2747.20	1/9.8
	180x60	87360	1077	87893.58	1/10.9

As mentioned before, the objective function can also be defined as twice the work done by external loadings, since the external work is equivalent to the strain energy for a conservative mechanical system. However, for the first and second-order versions of the finite-volume theory, there is a residual difference when a not sufficiently refined mesh is employed. Figure 9 shows the optimum topologies obtained when the structural compliance function is defined as twice the work done by external loadings. Table 5 presents the investigated numerical aspects during the optimization process.

The obtained optimum topologies are similar from one approach to another, although the obtained numerical aspects are worst when the external work is employed, presenting an increase in the number of iterations occasioned by a reduction in the damping factor. Consequently, it is also registered an increase in the computational cost. Thus, the objective function is better estimated when the structural compliance function is defined as twice the total strain energy.

Figure 10:
Optimum topologies for the MBB beam analysis by evaluating the compliance using the strain energy (filtering).

The topology optimization process is also performed, employing a sensitivity filter for mesh-independency. Figure 10 shows the obtained optimum topologies when the filtering technique is employed, which controls the mesh dependence, in the case of the finite-volume theory approaches, and the checkerboard effect and mesh dependence, in the case of the finite element method. Also, the optimum topologies presented in Figure 10 are practically the same for the different approaches.

Table 6 presents the numerical results obtained for the convergence analysis of the different employed approaches, considering the application of the mesh-independency filter. In general, the numbers of iterations are similar for the different approaches, with the Q4 approach showing more substantial differences in comparison to the other ones. In terms of computational cost, the Q8 approach has presented the highest processing time, followed by the second-order finite-volume theory. The Q4 approach has presented the lowest computational cost, followed by the zeroth-order finite-volume theory. For the current example, the stiffest structure was obtained for the Q8 approach, presenting the smallest compliance, where the values shown in Table 6 are relative to the optimum topology obtained by the Q8 element approach. In the relative compliance calculation, the compliance is evaluated employing the Q8 finite element for the optimum topologies obtained by the different approaches, for a fair comparison.

Table 6:
Convergence analysis for the MBB beam by evaluating the compliance using the strain energy (filtering).

Analysis	Mesh	NDOF	Number of Iterations	Processing Time (s)	Relative compliance
Q4	90x30	5642	603	363.85	1.00082
	180x60	22082	367	5109.44	1.00068
	270x90	49322	590	55524.63	1.00057
Q8	90x30	16682	447	2145.40	1.00000
	180x60	65762	650	51972.51	1.00000
	270x90	147242	745	324878.83	1.00000
FVT0th	90x30	11040	482	746.57	1.00333
	180x60	43680	664	19848.28	1.00089
	270x90	97920	761	110059.62	1.00063
FVT1st	90x30	16560	472	1278.69	1.00252
	180x60	65520	684	31081.14	1.00051
	270x90	146880	754	202667.73	1.00033
FVT2nd	90x30	22080	428	2464.42	1.00459
	180x60	87360	672	56573.36	1.00262
	270x90	195840	711	286714.04	1.00175

6 CONCLUSIONS

The topology optimization for compliance minimization algorithms based on the three versions of the generalized finite-volume theory has shown to be efficient, especially in the absence of filtering technique, where the algorithms have demonstrated the checkerboard-free property. This efficiency has its origins in the satisfaction of continuity conditions in a surface-averaged sense between adjacent subvolumes, which provides interfacial connections among the subvolumes. On the other hand, the checkerboard pattern effect is definitively a problem for the Q4 and Q8 elements. In those cases, this issue appears due to the finite element method's assumptions, such as the satisfaction of equilibrium and continuity conditions through the nodes, which provides nodal connections, resulting in checkerboard regions. In the case of the higher-order versions of the finite-volume theory, the evaluation of the compliance function using the strain energy shows to be more efficient than using the work done by external loadings.

The continued penalization scheme is adopted during the optimization, guaranteeing a gradual convergence for the overall process. In the absence of filtering techniques, the OC method's damping factor is adjusted to avoid divergence during the optimization process, since a non-maximum number of iterations is established. The damping factor was set up to be as close as possible to the value of 1/2 and avoid the oscillatory phenomenon during the optimization process. For the approaches that employ the mesh-independency filter, the damping factor was set up as 1/2, providing a faster convergence.

The sensitivity filter is employed to solve the mesh dependence and length scale problems. In the case of the finite element method, this filtering technique is employed to avoid the formation of checkerboard regions additionally. In terms of processing time, the approach based on the Q4 element is the most efficient, while the approach based on the Q8 element is usually the less efficient, with the finite-volume theory exhibiting the intermediates values, with higher processing times for the higher-order versions. For the cantilever beam example, the optimum topologies obtained by the standard (or zeroth-order) finite-volume theory present fewer bars and most with lower slenderness when compared with the topologies obtained by the approaches based on the Q4 and Q8 elements, which are desired features for manufacturing. Even though the Q4 approach shows the fastest convergence for this example, the obtained optimum topologies are less desirable, presenting more bars and most with higher slenderness when compared with the topologies obtained by the other approaches.

It is adopted a unique expression to evaluate the filter radius for all analyzes, considering only the neighbor elements/subvolumes (with shared nodes). Different values for the filter radius can affect the obtained topologies, but this investigation can be conducted in future works.

Based on the obtained results, the continuation of this investigation is justified by exploring the different aspects that evolve the finite-volume theory, especially in the case of heterogeneous materials with periodic microstructure, where the finite-volume theory has shown to be also efficient.

Acknowledgments

The authors acknowledge the financial support provided by the Brazilian federal agency CNPq.

References

• Aboudi, A., Pindera, M.-J., Arnold, S. M. (1999). Higher-order theory for functionally graded materials. Composites Part B: Engineering 30(8): 777-832.
• Araujo, M.V.O., Lages, E.N., Cavalcante, M.A.A. (2020). Checkerboard free topology optimization for compliance minimization applying the finite-volume theory. Mechanics Research Communications 108: 103581.
• Bansal, Y., Pindera, M.-J. (2003). Efficient reformulation of the thermoelastic higher-order theory for functionally graded materials. Journal of Thermal Stresses 26(11-12): 1055-1092.
• Bansal, Y., Pindera, M.-J. (2005). A second look at the higher-order theory for periodic multiphase materials. Journal of Applied Mechanics 72(2): 177-195.
• Bendsøe, M. P., Kikuchi, N. (1988). Generating optimal topologies in structural design using a homogenization method. Computer Methods in Applied Mechanics and Engineering 71(2): 197-224.
• Bendsøe, M. P. and Sigmund, O. (2003). Topology optimization: theory, methods and applications, Springer-Verlag (Berlin).
• Cavalcante, M. A. A., Marques, S. P. C., Pindera, M.-J. (2007a). Parametric formulation of the finite-volume theory for functionally graded materials - Part I: analysis. Journal of Applied Mechanics 74(5): 935-945.
• Cavalcante, M. A. A., Marques, S. P. C., Pindera, M.-J. (2007b). Parametric formulation of the finite-volume theory for functionally graded materials - Part II: numerical results. Journal of Applied Mechanics 74(5): 946-957.
• Cavalcante, M. A. A., Marques, S. P. C., Pindera, M.-J. (2008). Computational aspects of the parametric finite-volume theory for functionally graded materials. Computational Materials Science 44(2): 422-438.
• Cavalcante, M. A. A. and Pindera, M.-J. (2012a). Generalized finite-volume theory for elastic stress analysis in solid mechanics - part I: framework. Journal of Applied Mechanics 79(5): 051006.
• Cavalcante, M. A. A. and Pindera, M.-J. (2012b). Generalized finite-volume theory for elastic stress analysis in solid mechanics - part II: results. Journal of Applied Mechanics 79(5): 051007.
• Cavalcante, M. A. A. and Pindera, M.-J. (2014a). Generalized FVDAM theory for periodic materials undergoing finite deformations - part I: framework. Journal of Applied Mechanics 81(2): 021005.
• Cavalcante, M. A. A. and Pindera, M.-J. (2014b). Generalized FVDAM theory for periodic materials undergoing finite deformations - part II: results. Journal of Applied Mechanics 81(2): 021006.
• Chen, Q., Wang, G., Chen, X. (2018). Three-dimensional parametric finite-volume homogenization of periodic materials with multi-scale structural applications. International Journal of Applied Mechanics 10(4): 1850045.
• Collet, M., Bruggi, M., Duysinx, P. (2017). Topology optimization for minimum weight with compliance and simplified nominal stress constraints for fatigue resistance. Structural and Multidisciplinary Optimization 55(3): 839-855.
• Díaz, A., Sigmund, O. (1995). Checkerboard patterns in layout optimization. Structural and Multidisciplinary Optimization 10(1): 40-45.
• Ferretti, E. (2020). DECM: A Discrete Element for Multiscale Modeling of Composite Materials Using the Cell Method. Materials 13(4): 1-38.
• Gattu, M., Khatam, H., Drago, A. S., Pindera, M.-J. (2008). Parametric finite-volume micromechanics of uniaxial, continuously reinforced periodic materials with elastic phases. Journal of Engineering Materials and Technology 130(3): 031015.
• Haber, R. B, Jog, C. S., Bendsøe, M. P. (1996). A new approach to variable-topology shape design using a constraint on perimeter. Structural and Multidisciplinary Optimization 11(1-2): 1-12.
• Khatam, H., Pindera, M.-J. (2009). Parametric finite-volume micromechanics of periodic materials with elastoplastic phases. International Journal of Plasticity 25(7): 1386-1411.
• Khatam, H., Pindera, M.-J. (2010). Plasticity-triggered architectural effects in periodic multilayers with array wavy microstructures. International Journal of Plasticity 26(2): 273-287.
• Lopes, C. G., Santos, R. B., Novotny, A. A. (2015). Topological derivative-based topology optimization of structures subject to multiple load-cases. Latin American Journal of Solids and Structures 12(5): 834-860.
• Ma, Z.-D., Kikuchi, N., Hagiwara, I. (1993). Structural topology and shape optimization for a frequency response problem. Computational Mechanics 13(3): 157-174.
• Michell, A. G. M. (1904). The limits of economy of material in frame structures. The London, Edinburgh, and Dublin Philosophical Magazine and Journal Science 8(47): 589-597.
• Montes, M. A. (2016). Topology optimization algorithms for the solution of compliance and volume problems in 2D, Master Thesis (in English), Investigation Center in Mathematics, Mexico.
• Pomezanski, V., Querin, O. M., Rozvany, G. I. N. (2005). CO-SIMP: extended SIMP algorithm with direct corner contact control. Structural and Multidisciplinary Optimization 30(2): 164-168.
• Poulsen, T. A. (2002). A simple scheme to prevent checkerboard patterns and one-node connected hinges in topology optimization. Structural and Multidisciplinary Optimization 24(5): 396-399.
• Rozvany, G. I. N. (2009). A critical review of established methods of structural topology optimization. Structural and Multidisciplinary Optimization 37(3): 217-237.
• Rozvany, G. I. N., Querin, O. M., Gaspar, Z., Pomezanski, V. (2003). Weight-increasing effect of topology simplification. Structural and Multidisciplinary Optimization 25(5-6): 459-465.
• Shobeiri, V. (2016). Structural topology optimization based on the smoothed finite element method. Latin American Journal of Solids and Structures 13(2): 378-390.
• Sigmund, O. (2001). A 99 line topology optimization code written in MatLab. Structural and Multidisciplinary Optimization 21(2): 120-127.
• Sigmund, O. (2007). Morphology-based black and white filters for topology optimization. Structural and Multidisciplinary Optimization 33(4-5): 401-424.
• Sigmund, O., Petersson, J. (1998). Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima. Structural and Multidisciplinary Optimization 16(1): 68-75.
• Talischi, C., Paulino, G. H., Pereira, A., Menezes, I. F. M. (2012). PolyTop: a Matlab implementation of a general topology optimization framework using unstructured polygonal finite element meshes. Structural and Multidisciplinary Optimization 45(3): 329-357.
• Vieira, C. S., Marques, S. P. C. (2019). A new three-dimensional finite-volume model for evaluation of thermal conductivity of periodic multiphase composites. International Journal of Heat and Mass Transfer 139: 412-424.
• Versteeg, H. K. and Malalasekera, W. (2007). An introduction to computational fluid dynamics: the finite volume method, Prentice-Hall (New York).
• Wang, H., Liu, J., Qian, X., Fan, X., Wen, G. (2017). Continuum structural layout in consideration of the balance of the safety and the properties of structures. Latin American Journal of Solids and Structures 14(6): 1143-1169.
• Wang, S., Wang, M. Y. (2006). Radial basis functions and level set method for structural topology optimization. International Journal for Numerical Methods in Engineering 65(12): 2060-2090.
• Zhong, Y., Bansal, Y., Pindera, M.-J. (2004). Efficient reformulation of the thermal higher-order theory for FGMs with locally variable conductivity. International Journal of Computational Science and Engineering 5(4): 795-831.

Edited by

Guest Editors:

Volnei Tita and Nicholas Fantuzzi.

Publication Dates

• Publication in this collection
  09 Nov 2020
• Date of issue
  2020

History

• Received
  07 Apr 2020
• Reviewed
  16 July 2020
• Accepted
  29 Aug 2020
• Published
  01 Sept 2020