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

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.


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 higherorder 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 volumeaverage 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 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.
Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 3/21 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.

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 Pindera (2003 and 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 Pindera (2009 and2010) 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 surfaceaveraged 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 surfaceaveraged sense. Thus, the additional coefficients of the displacement field can be expressed in terms of the new surfaceaveraged 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 vertexcentered 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-Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 4/21 scale modeling of heterogeneous materials with periodic microstructure (Ferretti, 2020). These approaches can also be explored to solve the checkerboard pattern of topology optimization.

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 1 − 2 plane with 0 ≤ 1 ≤ and 0 ≤ 2 ≤ , which is discretized in horizontal subvolumes and vertical subvolumes. The subvolume dimensions are and ℎ for = 1, … , , where = • is the total number of subvolumes. In this cartesian formulation of the generalized finitevolume 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): where = 1,2 and ( ) ( ) 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).  Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 5/21 Figure 2a illustrates the kinematic quantities associated with each face of a generic subvolume . Therefore, these quantities in terms of surface-averaged displacements, rotations, and curvatures can be defined, respectively, as follows 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 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 where � ( ) is the local surface-averaged displacement vector, � ( ) is the local surface-averaged rotation vector, � ( ) is the local surface-averaged curvature vector, ( ) is the vector containing the unknown coefficients related to the zerothorder finite-volume theory, ( ) is the vector formed by the unknown coefficients related to the first-order finite-volume theory, 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 Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 6/21 In the static analysis, the continuity conditions between the third and first faces of adjacent subvolumes are imposed as follows 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 where ̅ ( ) is the local surface-averaged traction vector, ̅ ( ) is the local surface-averaged normal traction first derivative vector and ̅ 2 ( ) is the local surface-averaged normal traction second derivative vector 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 where 1 ( ) = , 2 ( ) = ℎ , 3 ( ) = and 4 ( ) = ℎ are the subvolume edges lengths and ̅ ( ) ( ) is taken from Eq. (7) and can be expressed as where (2×16) ( , ) are submatrices of selected components of the matrix (16×16) ( ) related to the surface-averaged tractions acting on a face of the subvolume . Replacing Eqs. (4) and (7) in Eqs. (9) and (8) The vector (00) ( ) can be obtained from Eq. (10), which is given by Replacing Eq. (11) in Eq. (4), the following expression can be obtained: . Thus, the local system of equations for a generic subvolume can be obtained by replacing Eq. (12) in Eq. (7), which gives 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 where is the total number of degrees of freedom, ( ×1) and ( ×1) are the global surface-averaged static and kinematic vectors, respectively, and ( × ) is the global stiffness matrix evaluated by 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 finitevolume 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.

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 Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 8/21 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 where ( ) and � are the material and the reference domain volumes, respectively, is the local displacement vector, 0 is the local stiffness matrix for a unitary relative density, is the relative density vector, is the penalty factor, is the prescribed volume fraction, is the minimum relative density, is the relative density associated with each element and 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 where is the current iteration, is the move-limit, is the damping factor and is given by where 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 is directly related to the method performance, once this affects the speed variation of (Montes, 2016). A high value for 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 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 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.

Mesh-independency filter
To avoid the occurrence of mesh dependency, it is suggested the modification of the elements' sensitivities by the following expression: where � is the convolution operator (weighting function) given as Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 9/21 where dist( , ) is the distance between the element center of and the element center of (Sigmund, 2001).
To consider the contribution only of the neighbor elements (with shared nodes), it is adopted a filter radius of = 1.01�( ) 2 + (ℎ ) 2 .

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 ; thus, this function can be expressed as where ( , ) is the total strain energy, � ( , ) is the specific strain energy, ( , ) is the stress tensor, ( ) is the strain tensor, ( ) is the stiffness tensor and is the reference domain. In the absence of body forces and considering = 0 in , the work done by external loadings can be defined as where are the traction components, are the displacement components, is the external surface subjected to external loadings, is the external surface with predicted displacements and = ∪ .
Applying Cauchy's law and the divergence theorem to Eq. (22), it follows where is the asymmetric rotation tensor. Considering the symmetry of , Eq. (23) can be written as As a result, the equivalence between the total work done by external forces and the total strain energy is observed only if ⁄ = 0. However, this is valid only for the zeroth-order finite-volume theory. Consequently, for the higherorder versions of the generalized finite-volume theory, it is verified the need to investigate the different aspects that involve mechanical energy evaluation.

Total strain energy for the generalized finite-volume theory
For linearly elastic materials, the total strain energy can be defined as where is the stress tensor, is the strain tensor and is the material stiffness tensor. Considering the displacement approximation presented in Eq. (1), the strain tensor of a subvolume can be described as Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 10/21 where ( ) � 1 ( ) , 2 ( ) � 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 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 where Besides, 02 ( ) , 12 ( ) , 20 ( ) and 21 ( ) 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 The strain energy for the lower order versions of the generalized finite-volume theory can be obtained by uncoupling the vector 2 ( ) , in the case of the first-order version, and 2 ( ) and ( ) , in the case of the zeroth-order version.

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 1 (2,4) � 2 ( ) � = 1(0) (2,4) + 2 ( ) 1(1) (2,4) + where 1( ) (2,4) are unknown coefficients of the horizontal displacement at the subvolume vertical faces. From Eq. (31), the surface-averaged kinematic quantities can be evaluated as Similarly, for linearly elastic materials, the normal traction acting on the subvolume vertical faces can be expressed by three Legendre polynomials as follows Using Eqs. (33) and (34) 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 where ( ) is the local resultant force vector, ( ) is the local resultant bending moment vector and ( ) 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 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. 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 (∆ = 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.

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 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 secondorder 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.   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. 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 Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 15/21 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 Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 16/21 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.

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 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 elementbased 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 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.  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. 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 finitevolume 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 Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e308 19/21 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.

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.