Mixed variational formulations for structural topology optimization based on the phase-field approach

We propose a variational principle combining a phase-field functional for structural topology optimization with a mixed (three-field) Hu-Washizu functional, then including directly in the formulation equilibrium, constitutive, and compatibility equations. The resulting mixed variational functional is then specialized to derive a classical topology optimization formulation (where the amount of material to be distributed is an \emph{a priori} assigned quantity acting as a global constraint for the problem) as well as a novel topology optimization formulation (where the amount of material to be distributed is minimized, hence with no pre imposed constraint for the problem). Both formulations are numerically solved by implementing a mixed finite element scheme, with the second approach avoiding the introduction of a global constraint, hence respecting the convenient local nature of the finite element discretization. Furthermore, within the proposed approach it is possible to obtain guidelines for settings proper values of phase-field-related simulation parameters and, thanks to the combined phase-field and Hu-Washizu rationale, a monolithic algorithm solution scheme can be easily adopted. An insightful and extensive numerical investigation results in a detailed convergence study and a discussion on the obtained final designs. The numerical results clearly highlight differences between the two formulations as well as advantages related to the monolithic solution strategy; numerical investigations address both two-dimensional and three-dimensional applications.


Introduction
Having assigned a design region, a load distribution, and suitable boundary conditions, the goal of structural topology optimization is to identify an optimal distribution of material within the design region. Classically, optimality is reached when the obtained material distribution minimizes a measure of structural compliance and satisfies mechanical equilibrium. The minimization of structural compliance is in general expressed in terms of the work done by the assigned (assumed to be constant) load distribution by the corresponding displacements.
Several methods have been developed in order to solve topology optimization problems. Originally, a discrete formulation was introduced where areas of dense material and voids are alternated without any transition region (Bendsøe, 1983). Known as 0-1 topology optimization problem, this first approach also leads to many difficulties both from an analytical and a numerical point of view (Sigmund and Petersson, 1998).
Possible alternative approaches are based on homogenization methods, where the macroscopic material properties are obtained from microscopic porous material characteristics (Allaire et al., 2004;Suzuki and Kikuchi, 1991), or on the Solid Isotropic Material Penalization (SIMP) method (Zhou and Rozvany, 1991), which consists of penalizing the density region, different from the void or bulk material, by choosing a suitable interpolation scheme for material properties at the macroscopic scale (Bendsøe and Sigmund, 1999;Bendsøe, 1983). Another approach used by many authors (cf., e.g., Burger, 2003;Osher and Santosa, 2001) is the level set method, that allows for topology changes but presents difficulties in creating holes. An alternative to SIMP and level set is based on the phase-field method, for the first time introduced by Bourdin and Chambolle (2003) and successively employed by Burger and Stainko (2006) for stress constrained problems as well as by Takezawa et al. (2010) for minimum compliance and eigenfrequency maximization problems. More recently, Penzler et al. (2012) have solved nonlinear elastic problems by means of the phase-field approach, while Dedè et al. (2012) have applied this method in the context of isogeometric topology optimization.
The phase-field approach delivers an efficient method to solve the free boundary problem stemming from the fact that the boundary of the region filled by the material is unknown. Phase-field topology optimization penalizes an approximation of the interface perimeter in such a way that, by choosing a very small positive penalty term, one can obtain a sharp interface region separating solid materials and voids (Blank et al., 2014a). The phase-field approach to structural optimization problems has been recently used by different authors (cf., e.g., Auricchio et al., 2019;Carraturo et al., 2019;Blank et al., 2014b), the main advantage being the fact that it allows to handle topology changes as well as nucleation of new holes. Moreover, especially for three-dimensional applications, another advantage is that the phase-field approach produces regular patterns in the design domain that require little post-processing effort to interpret results (Deaton and Grandhi, 2014). On the other hand, density-based topology optimization approaches (e.g., based on SIMP or power-laws) require post-processing filtering for translating results in manufacturable data, which may result case-and user-dependent and might lead to undesirable effects, like an artificial increase of the final volume fraction (Zegard and Paulino, 2016).
Although promising results have been obtained on phase-field topology optimization, there are still some key aspects which merit further investigation. Firstly, the adopted values of phase-fieldrelated parameters affect both the obtained topology distribution and the convergence behavior of the numerical solution; hence, selecting proper values of these parameters for the specific problem at hand requires a time-consuming preliminary set of investigations.
Furthermore, traditional schemes for topology optimization problems generally fix the amount of material to be distributed in the design region, amount which acts as a global constraint in the optimization process. This gold-standard strategy is easy to implement, since it adds only a single global equation in the computation; however, when the evolution of the topology is itself the solution of a finite element computation, the presence of such a global constraint requires to introduce global variables, altering the intrinsic element-based subdivision of the spatial discretization and introducing a coupling term between all elements. As an obvious consequence, the usual banded structure of the tangent stiffness matrix is altered and such an increase of the band-width of the tangent stiffness matrix might have consequences on the convergence rate of topology optimization procedures. Currently, the literature is missing a formulation which involves only local constraints for enforcing/minimizing the amount of material to be distributed.
An additional issue of high interest is related to classically adopted numerical solution schemes. These are generally based on staggered algorithms: the evolution of the phase-field variable is obtained without the respect of the equilibrium, which is restored at the successive incremental step wherein the phase-field variable is instead frozen. In the framework of topology optimization, such an approach is referred to as Nested Analysis and Design (NAND). Allowing for the use of existing FE packages, NAND is easy and convenient to implement and it is often preferred also as a consequence of the presence of the global constraint. Nevertheless, NAND can be computationally demanding, especially in the presence of non-linearities and for large-scale problems. On the other hand, a Simultaneous Analysis and Design (SAND) approach treats the optimization design (e.g., phase field) and equilibrium (e.g., displacement) unknowns simultaneously. In other words, the combined optimization and equilibrium problem is solved in a monolithic way by considering the evolution of topology and equilibrium conditions in a single optimization routine. Despite the fact that SAND may be competitive with respect to conventional NAND schemes, to the best of authors' knowledge, SAND (monolithic) approaches are still largely to be explored, as compared with the more classical NAND (staggered) solution strategies.
Finally, the phase-field variable should be physically bounded in the range [0, 1]. To enforce this constraint, the phase-field variable obtained from the solution of the topology evolution equation (which possibly lies outside the admissible bound) is generally projected on the admissible set at each step of the iterative solution procedure. This strategy obviously alters stationary. A variational formulation allowing for the direct enforcement of the physical bounds, without any ad hoc algorithmic solution, would be more consistent. Accordingly, in the described context of literature and open problems, the goals of the present paper are the following ones.
-To propose a mixed variational principle combining a phase-field functional for topology optimization with a (three-field) Hu-Washizu functional, then including directly in the formulation equilibrium, constitutive, and compatibility equations. The variational formulations includes also a penalization term to directly enforce that the phase-field variable lies within the admissible physical bound [0, 1]; -To specialize the mixed variational functional to a classical formulation for the topology optimization problem, referred to as formulation with volume constraint. Here, the amount of material to be distributed within the design domain is imposed a priori, acting as a global constraint and operating on the global stiffness matrix of the corresponding finite element formulation; -To propose a novel formulation for the topology optimization problem, leading at the same time to a volume minimization. The amount of material is not imposed a priori by means of a global constraint, but rather related to a cost of the material, which is minimized together with the structural compliance. Such an approach avoids the introduction of a global constraint, respecting the convenient local nature of the finite element discretization; -To implement the proposed variational formulations in the context of mixed finite elements.
Thanks to the combined phase-field and Hu-Washizu rationale, a SAND monolithic algorithm is proposed, allowing to compute the evolution law of the phase-field variable under the respect of equilibrium; -To trace general guidelines for fixing the values of phase-field-related simulation parameters.
In addition, a theoretical estimate for setting the cost of the material (which plays the role of a mass penalty parameter) in the volume minimization formulation as function of a target desired solution is proposed and verified; -To perform a comparative analysis in terms of both computational efficiency and obtained final designs for monolithic and staggered solution strategies as well as for the two investigated formulations.
-To show the effectiveness of the proposed volume minimization formulation and its finite element mixed implementation both in two-dimensional and three-dimensional applications.

Introductory settings
We indicate the design region with Ω ⊂ R n (with n = 2 or n = 3), characterized by a volume V , defined as: The material distribution within Ω is described by a scalar phase variable φ, ideally a binary field (i.e., either 0 or 1), with φ = 0 corresponding to the absence of material (i.e., void or no material), and φ = 1 corresponding to the presence of material. After topology optimization, the material distributed inside the design region Ω occupies a (dimensionless) volume fraction v of the total volume V , defined as: However, since a binary (i.e., either 0 or 1) phase parameter φ would call for sharp interfaces (between the region with material, i.e., with φ = 1, and the region with no material, i.e., with φ = 0) and since problems with infinitely sharp interfaces are in general difficult to numerically treat and solve, as classically done, we consider the phase parameter φ as a continuous real value field, with values in the interval [0, 1]. Ideally, the optimal phase variable φ should be close to an ideal binary field.

Elastic problem
Limiting the discussion to the case of linear elastic materials, given a material distribution φ : Ω → [0, 1], the elastic problem (equilibrium, constitutive, and compatibility field equations, in combination with boundary conditions) can be expressed as: with σ the stress tensor field, b the applied body-force density vector field per unit volume, C(φ) a positive-definite material elastic tensor function of the phase variable φ, ε the strain tensor field, u the displacement vector field, div and ∇ s the divergence and the symmetric gradient operators, Γ D the region of ∂Ω where we apply Dirichlet boundary conditions (assumed to be homogeneous), Γ N the region of ∂Ω where we apply Neumann boundary conditions with t the surface load vector and n the outward-pointing versor normal to Γ N . The elastic problem can be classically condensed in just one field equation in terms of displacements (referred to as Navier's equation), as follows: (2.4) Assuming voids to be modeled as a very soft material, in the present work we adopt the following expression for C(φ): where C A is a constant positive-definite material elastic tensor (assumed to describe an isotropic behaviour with Young's modulus E A and Poisson's ratio ν A ) corresponding to the solid dense material, p can be any positive value, and 0 < δ 1 governs the low (but non-null) stiffness of the voids. Equation (2.5) is such that C(1) = C A and C(0) → δC A for p → +∞; as an example, p = 10 already corresponds to C(0) − δC A < 10 −4 .
It is worth highlighting that the particular choice of C(φ) from Equation (2.5) is different from state-of-the-art expressions, which generally might read as one of the following: (2.6b) As shown in Fig. 1 for their scalar counterparts, these lasts expressions either bring to nonmonotonic expressions (i.e., C 1 (φ)) or fast diverge as soon as φ exceeds the admissible range [0, 1] (i.e., C 2 (φ)). On the one hand, a monotonically increasing behaviour would be beneficial in terms of stability of the ensuing numerical scheme and it would be more consistent from the physical point of view. On the other hand, a smooth behaviour outside the admissible range would be beneficial because the property φ ∈ [0, 1] is sometimes violated in numerical iterative solution schemes, often requiring an ad hoc truncation of the stiffness matrix or of the phase-field variable, and hence introducing non-smoothness effects. The expression proposed in Eq. (2.5) and adopted in this work allows to overcome drawbacks inherited by definitions in Eq. (2.6) (see Fig. 1). In fact, C(φ) turns out to be monotonically increasing, everywhere-defined, and uniformly positive-definite for all φ, being at the same time smooth with all its derivatives. In addition, our choice for C(φ) entails that the elastic energy density (φ, ε) → C(φ)ε : ε/2 is convex, contributing to the stability of the approximation. (2.6a); f 2 C (φ), scalar counterpart of Eq. (2.6b); and f C (φ), scalar counterpart of the proposed expression (2.5). Parameters: C A = 10, δ = 10 −3 and p = 10.

Topology optimization
The goal of topology optimization is then to design an optimal structure, i.e., to identify the optimal material distribution φ, such that structural compliance is minimized when solving the elastic problem. We introduce the structural compliance C, with the unit of measure of a work, as: where b and t are the actual body and surface loads (assumed to be constant), while u is the displacement, satisfying the elastic problem (2.4). Accordingly, the topology optimization problem can be written as: It should be remarked that in (2.4) the boundary conditions are defined on a fixed subset of ∂Ω. Once a solution is numerically obtained, one has to check whether φ = 0 at least on some portion of Γ D and Γ N . This generally follows from minimality. In particular, it is fulfilled by all our computations.

Phase-field topology optimization
As mentioned, the phase parameter φ is allowed to take values in the interval [0, 1] and, accordingly, the region where 0 < φ < 1 corresponds to a smooth material interface, representing a transition from void to solid material. To introduce a measure of the extension or the thickness of such a transition interface and to favor fields taking values as close as possible to 0 or 1, we introduce a new functional P to be minimized, defined as: where γ φ > 0 is the thickness penalization parameter (with the unit of measure of a length) and ψ 0 is the double-well potential function, defined as: (2.10) For small values of γ φ , the minimization of P(φ) penalizes oscillatory material distributions and favours continuous fields φ taking values close to 0 or 1. For γ φ → 0, the quantity P(φ) converges in the sense of Γ -limits to the perimeter functional (cf. Modica (1987)), namely a measure of the perimeter of the interfaces between regions with material (φ = 1) and regions with no material (φ = 0). In the following we augment the compliance functional C by a multiple of P, introducing the functional J , defined as: where the perimeter stiffness parameter κ φ > 0 (force divided by length) measures the mutual relevance of structural compliance and perimeter in the minimization. Accordingly, the phase-field topology optimization problem can be written as: 2.5 Phase-field constraint We start recalling that ideally the phase-field variable φ should be as much as possible a binary field but also that any value outside the interval, i.e., φ / ∈ [0, 1], is physically meaningless and it can be the source of numerical instabilities. The minimisation of J in Eq. (2.12) already favors φ to be close to 0 or 1 through the presence of ψ 0 , but φ ∈ [0, 1] just in the limit as γ φ → 0; for all γ φ > 0 (which is the case of any numerically-obtained solution) the constraint φ ∈ [0, 1] does not necessarily hold through the minimisation of J .
With these concepts in mind, it is clear why in most iterative schemes available in the topology optimization literature the variable φ is projected onto the admissible range [0, 1] at each step, however violating the stationary conditions of the corresponding adopted functional.
As an alternative, we propose to further augment the functional J by a contribution B, defined as: and modulated by a bounding stiffness parameter κ b > 0 (force per unit area).

Formulation with volume constraint
As classically found in the literature, we may now constrain the topology optimization problem so that the portion v of material to be distributed inside the design domain Ω is equal to an a priori (2.14) The corresponding topology optimization problem reads as: Before moving on, we wish to remark that problem (2.15) can be equivalently rewritten as Our specific choice (2.5) for C makes the functional under minimization to be convex, up to the lower-order ψ 0 -term in P(φ). We propose to solve the above minimization problems by investigating the stationarity of the Lagrangian L vc , defined as: where E el is the elastic energy, defined as: The superscript vc in L vc recalls that we are here considering a volume-constrained formulation.

Formulation with volume minimization
Instead of prescribing a priori a fraction of material v to be distributed in the design domain Ω, we may aim at exploring what is the minimum amount of material v we could distribute. In terms of minimization, a corresponding topology optimization problem would read where κ v > 0 is a volume penalty parameter (force per unit area), representing for example a measure of the "cost" of the material per unit volume. It is interesting to emphasize that in (2.19) we are not minimizing the amount of material v to be distributed, but rather Ω φ 2 dΩ; this choice makes the problem more stable and it can be proved to be equivalent to minimize v, as long as φ exclusively takes values 0 or 1. As for the volume-constrained case, problem (2.19) can be equivalently rewritten as Again, the choice (2.5) for C makes the functional under minimization to be convex, up to the lower-order ψ 0 -term in P(φ). We propose to solve the topology optimization problem (2.19) by looking at the stationarity of the Lagrangian L vm , defined as: The superscript vm in L vm recalls that we are looking at a volume-minimization formulation.

Formulation with volume constraint
To compute a stationary point of the Lagrangian L vc (2.17), we set its variations to zero, i.e.: (3.23) We recall that the last three relations in Eq. (3.22) are, respectively, the equilibrium, the compatibility, and the constitutive equations, which are then automatically enforced through the stationarity of L vc . Equations (3.22) can be rewritten in a residual form as: To solve such a residual system using an iterative scheme, a good starting point could be to set φ(X) = c ≤ 1 ∀X ∈ Ω, i.e., to consider a design domain homogeneously covered by a partially dense material. Natural choices are c = 0.5 or c = v.
In addition, one can follow a classical Allen-Cahn-like approach by introducing a gradient-flow dynamics. To do so, the phase-field variable φ is assumed to be dependent on a pseudo-time variable t, i.e., φ = φ(X, t), and the new Allen-Cahn volume-constraint Lagrangian L vc AC is introduced: where the parameterτ φ > 0 (unit of measure of force times squared time per unit area) represents viscosity with respect to the pseudo-time t. Exploring the stationary condition of L vc AC with respect to φ (assuming δφ to be independent from t) returns the condition: At the time-discrete level, denoting with a subscript n quantities evaluated at the previous time t n and with no subscript quantities evaluated at the current time t n+1 , we obtain: with ∆t = t n+1 −t n and the discrete viscosity parameter τ φ =τ φ /∆t (unit of measure of force times time per unit area) introduced such to have the discretized temporal derivative in the evolution equation (3.27). In a more explicit format, (3.28) We may note that the Allen-Cahn term contributes to the monotonicity of the map: delivering enhanced stability to the numerical scheme, although the map φ → D δφ L vc AC (φ). Indeed, such map turns out to be monotone, apart from the polynomial ψ 0 -term and the bilinear term φb · u in the compliance. By choosing τ φ /∆t large with respect to 1/γ φ , the Allen-Cahn term dominates the ψ 0 -term. In particular, the function In the literature it is often chosen τ φ = γ φ (modulo fixing dimensions). This choice is of course compatible with Eq. (3.29), as long as ∆t is taken to be small enough. Nevertheless, motivated by Eq. (3.29), the following relationship between τ φ and other parameters in the functional will be introduced: where T φ is a characteristic time constant. Substituting Equation (3.22) 1 with (3.28) we obtain the new residual equation corresponding to the formulation with volume constraint in an Allen-Cahn approach

Formulation with volume minimization
To compute a stationarity point of the Lagrangian L vm (2.21), we set its variations to zero, i.e.: We recall again that the last three relations in Eq. (3.32) are, respectively, the equilibrium, the compatibility, and the constitutive equations, which are then automatically enforced through the stationarity of L vm . The term E el ,φ (φ, ε) in Eq. (3.32) 1 is given in Eq. (3.23). Equations (3.32) can also be rewritten in a residual form as: (3.33) To solve such a residual system using an iterative scheme, as for the formulation with volume constraint, we can follow an Allen-Cahn-like approach, i.e., introducing a gradient-flow dynamics.
To do so, the phase-field variable φ is again assumed to be dependent on a pseudo-time variable t, i.e., φ = φ(X, t), and a new Lagrangian L vm AC is introduced: The stationarity of L vm AC with respect to φ (assuming δφ to be independent from t) corresponds to (3.36) The Allen-Cahn term again contributes to the monotonicity of φ → D δφ L vm AC (φ). Compared with the volume-constraint case, the situation is here more favorable, since the κ v term is also contributing to convexity. In particular, the function (3.37) In particular, comparing (3.37) with (3.29), one has that from (3.37) convexity still holds for τ φ = 0 (no Allen-Cahn regularization), as long as κ v is large enough. Also for this formulation, τ φ will be defined as function of other parameters as in Eq. (3.30). Substituting Equation (3.32) 1 with (3.36) we obtain the new residual system corresponding to the formulation with volume constraint in an Allen-Cahn approach 3.3 Solution algorithm, convergence criteria, and output quantities In a time-discrete setting, at each time instant t n+1 we need to solve the non-linear residual equation (3.31) for the formulation with volume constraint and the non-linear residual equation (3.38) for the formulation with volume minimization. Furthermore, assuming to properly solve at each time instant t n+1 the corresponding non-linear residual problem, we need to follow the time-discrete dynamics to convergence to a stationary solution. In the following, stationarity is measured by controlling changes in time of the material distribution φ and of the displacement field u. Accordingly, we adopt an algorithm with a double iteration loop, i.e., an external iteration loop on time (controlling a norm on material distribution φ and displacement field u changes) and an internal iteration loop (controlling a norm on the satisfaction of the residual problem).
The stopping criterion for the external iteration loop on time is defined on the basis of the following error relative to the Allen-Cahn procedure, or briefly Allen-Cahn error: where: with T φ the time constant introduced in Eq. (3.30) and A ∆u a normalization constant, computed as: Here,ū sol is the reference displacement field obtained by solving the equilibrium problem with a fixed φ = 1 everywhere in the design domain. The iterations in time are stopped when E AC < c conv AC , with c conv AC an Allen-Cahn convergence parameter to be set. At each time instant t n+1 , the internal iteration loop consists in the solution of the residual systems (3.31) and (3.38) by means of the Newton-Raphson method which is implicit since time discretization has been performed by means of a backward Euler strategy. Starting from the functionals L, derivatives required for computing the residual vector R and the tangent system matrix D (i.e., linearization of R) are performed by means of the Mathematica package AceGen (Korelc and Wriggers, 2016), allowing for combined symbolic-numeric programming. The solution of the finite element resulting system is obtained by means of the Mathematica package AceFEM (Korelc and Wriggers, 2016). The stopping criterion for the internal iteration loop is defined on the basis of the norm of the residual vector R corresponding to the governing equations, i.e.: The iterations on the residual equations are stopped when E R < c res R , with c res R a residual convergence parameter to be set, returning the updated solution of primary variables at time instant t n+1 . Divergence is obtained when a maximum number N max iter = 15 of Newton-Raphson iterations is reached. In this case, the time-step size ∆t is decreased and the internal loop at t n+1 is started again. The tuning of the time-step size is performed throughout the solution by means of a path-following procedure with an adaptive time stepping which returns the optimal ∆t around a reference value ∆t 0 , here chosen in the range ∆t ∈ (10 −5 , 10 3 )∆t 0 , (Korelc and Wriggers, 2016).
An overview of the solution algorithm is reported in Algorithm 1. It is noteworthy that, in our approach, we do not solve residual equations in a staggered form (as sometimes proposed in the literature) but we numerically iterate in an implicit form on all residual equations, which means that our approach is fully implicit and monolithic. In other words, the implemented solution strategy corresponds to a Simultaneous Analysis and Design (SAND) approach. For the sake of comparison, results will be compared also by adopting a Nested Analysis and Design (NAND) approach, which corresponds to a staggered implementation between the phase-field evolution problem and the mechanical equilibrium conditions. This is achieved by assuming that: -The elastic energy rate term in Eqs. (3.22) 1 and (3.32) 1 is computed at the previous time step, that is E el ,φ (φ n , ε n ) (see Eq. (3.23)); -The constitutive response considered in the mechanical equilibrium balance (i.e., Eqs. (3.22) 5 and (3.32) 4 ) is computed on the basis of the phase-field computed at the previous time step, that is C(φ n ). For the staggered NAND approach, the phase-field variable is projected within the admissible range [0, 1] at each step, in agreement with existing literature (e.g. Carraturo et al., 2019). This projection has been verified to be necessary for convergence issues, since otherwise C(φ n ) takes values leading to instability/divergence of the numerical iteration scheme. On the other hand, for the monolithic SAND approach, the phase-field variable is not ad hoc projected at each step. This choice is more sound from the theoretical point of view since the bounding constraint is already included in the variational functional, and obtained results will show that it does not lead to divergence issues.
Algorithm 1: Implicit and monolithic solution algorithm for the phase-field topology optimization problem with an Allen-Cahn-like strategy (SAND approach). Once a converged final solution is obtained (and in particular the corresponding phase-field φ sol , displacement u sol , and stress σ sol ), the solution quality is expressed by computing and evaluating the following three quantities: -A scalar measure of the structural displacement U sol , obtained by dividing the final value of the compliance C sol = C(φ sol , u sol ) with the total resultant F tot of the applied loads: -The diffused-perimeter measure P sol = P(φ sol ).
To compare different solutions (in particular corresponding to different material distributions) from an engineering point of view, we also compute the Von Mises stress distribution σ vm (X), the maximum Von Mises stress σ max vm , and the average Von Mises stress σ avg vm , defined respectively as: where s sol is the deviatoric part of the stress tensor σ sol at solution, · is the standard Euclidean norm, V sol is the effective total volume occupied by the material at solution, i.e.

Finite-element approximation and parameters settings
To solve the problem under investigation, we introduce a space discretization, considering a standard decomposition of the design domain Ω into a set of non-overlapping finite elements, and for each single element we introduce the approximations described in the following. Both twodimensional (2D) and three-dimensional (3D) examples will be considered.
-The displacement field u is approximated in a C 0 -continuous element-by-element form as: where N i (ξ) are Lagrange-type shape functions defined on the parent coordinate system ξ andû i are the nodal degrees of freedom, with n e n being the total number of nodes for each element. A standard isoparametric mapping from the parent coordinate system ξ to the spatial (physical) coordinate system X is implemented.
-The phase field φ is approximated in a C 0 -continuous element-by-element form as: whereφ i denotes nodal degrees of freedom and N i , n e n are defined above. -The stress field σ is approximated in a discontinuous element-by-element form, introducing element unknowns collected in vectorσ, (Djoko et al., 2006). In the parent coordinate system ξ, the second-order stress tensor σ ξ is interpolated as: where vec( · ) denotes the vector Voigt representation of stresses, and N σ represents the shape functions for the stresses. In the spatial (physical) coordinate system X, the second-order stress tensor σ is computed from the one in the parent domain, σ ξ , via: where T is a second-order transformation tensor which allows to map stresses from the parent ξ to the spatial X cartesian space. The transformation matrix must: 1. Produce stresses in spatial cartesian space which satisfy the patch test (i.e., can produce constant stresses and be stable); 2. Be independent of the orientation of the initially chosen element coordinate system and numbering of element nodes (invariance requirement). For the second-order transformation tensor T , Pian and Sumihara (1984) proposed a constant array (to preserve constant stresses) deduced from the Jacobian J(ξ) associated to the geometrical mapping between the physical coordinates X and the parent coordinates ξ, computed at the centre of the element, i.e., T = J| ξ=0 = J 0 . In the present paper, we employ the normalized counterpartJ 0 of it: (4.47) -The strain field ε is approximated in a discontinuous element-by-element form, introducing element unknowns collected in vectorε, (Djoko et al., 2006). In the parent coordinate system ξ, the second-order strain tensor ε ξ is interpolated as: 48) where N ε represent the shape functions for the strains. In the spatial coordinate system X, the second-order strain tensor ε is computed from the one ε ξ through the transformation tensor T as 1 : where J(ξ) = Det(J(ξ)) and J o = Det(J o ). -For the case of the formulation with volume constraint (i.e., L vc ), the Lagrange multiplier λ is assumed as a constant global field, hence: λ =λ . (4.50) In numerical applications, four node quadrilateral elements (QUAD-4 or Q1 elements) are implemented for 2D applications, while eight node hexahedral elements (HEX-8 or H1 elements) for 3D applications. Therefore, Lagrange-type shape functions N i correspond to bilinear and trilinear polynomials, respectively for 2D and 3D applications. More information on the interpolation of the stress and strain fields, together with the specific forms of N σ and N ε employed in this work (Weisman, 1996;Cao et al., 2002;Djoko et al., 2006), are given in Appendix A. Static condensation of stress and strain variables is performed at element level to reduce the dimensions of the numerical problem, which then correspond to the one of a pure-displacement phase-field formulation.

Settings of simulation parameters
The performance of phase field formulations highly depend on the values of the chosen parameters. These can be divided in physical and simulation parameters.
On one hand, physical parameters are the material properties contained in the tangent stiffness matrix C and, for the volume constraint formulation L vc , the final amount of material distributed in the design domain.
1 This choice satisfies with Eq. (4.46) the principle of energy equivalence in the parent and spatial coordinate systems: On the other hand, simulation parameters are the thickness penalization parameter γ φ (with the unit of length), the perimeter stiffness parameter κ φ (with the unit of force per unit length), and the bounding stiffness parameter κ b (with the unit of force per unit area). Moreover, for the volume minimization formulation L vm , the final amount of material distributed in the design domain (i.e., v sol ) depends on another simulation parameter represented by the volume penalty parameter κ v (with the unit of force per unit area).
The following considerations are traced for properly settings the values of simulation parameters. As regards the volume constraint formulation L vc : -The thickness penalization parameter γ φ is chosen comparable to a typical element dimension h e in the mesh discretization, that is γ φ ≈ h e ; -The perimeter stiffness parameter κ φ is chosen such that the term κ φ /γ φ (governing the driving force which allows to obtain a binary-like solution for the phase field) is comparable with the strain-energy density rate e el ,φ in Eq. (3.23). The latter quantity is clearly not constant during the solution. A proper estimate can be anyway obtained as κ φ ≈ γ φē el ,φ , withē el ,φ = e el ,φ (1,ε) being the strain-energy density rate obtained on the "full material case". In detail,ε is the reference strain field obtained by solving the equilibrium problem with a fixed φ = 1 everywhere in the design domain; -The bounding stiffness parameter κ b should be at least 3 order of magnitudes greater than κ φ , that is κ b > 10 3 κ φ .
As regards the volume minimization formulation L vm , the same rules-of-thumb can be introduced for γ φ and κ b . Nevertheless, contrarily to the volume constraint formulation, additional care is needed in order to obtain an a priori estimate of the final volume fraction v sol . In fact, although the volume minimization principle inherits advantages from the design viewpoint (as the following results will show), a preliminary estimate of v sol would be useful for obtaining a reference design solution. Therefore, the following is proposed: -The volume penalty parameter κ v is determined on the basis of a target estimate of the final volume fraction v tar sol . To this aim, a function v tar sol (κ v ) is introduced such that the following three conditions are met: (4.51) The first two conditions enforce the limit conditions that a material with "null cost" (i.e., κ v = 0) will fill the entire domain design, while for a material with an extremely "high cost" (i.e., κ v → +∞) the domain design will tend to be empty. The last condition prescribes that, starting from "null material cost", the variation of v tar sol associated with an increase of κ v is inversely proportional toē el ,φ . Conditions (4.51) are met by the following definition of the which allows to fix κ v on the basis of the target desired value of v tar sol ; -Once κ v is fixed, the perimeter stiffness parameter κ φ is now chosen such that the term κ φ /γ φ ≈ κ v , or κ φ = κ v γ φ . This choice follows the same rationale of the volume constraint formulation, but it allows to tune κ φ as a function of the target volume fraction.

Numerical results
The following Section presents an extensive numerical investigation resulting in a detailed convergence study and a discussion on the obtained final designs. The numerical results clearly highlight differences between the two investigated formulations as well as advantages related to the monolithic solution strategy, with numerical simulations addressing both two-dimensional and three-dimensional applications. In particular, Section 5.1 presents the results of a 2D cantilever and performs a comparison between different formulations and solution approaches, while Section 5.2 extends the discussion to two 3D examples, i.e., a 3D cantilever and a 3D bridge (see Figure 2).
In all the simulations, null body forces are considered, e.g., b = 0. Moreover, the characteristic time constant is chosen as T φ = 1 s, the discrete viscosity parameter τ φ as in Eq. (3.30), the reference time step increment ∆t 0 = 10 −2 s, and the convergence parameter for the residual internal iteration loop fixed as c res R = 10 −8 J.

2D simulations: parametric analyses
The 2D cantilever is solved starting from a two-dimensional design domain of rectangular shape with dimensions a × b (see Fig. 2), assuming plane strain conditions, with out-of-plane thickness c. The domain is clamped on the left vertical side and loaded by a downward vertically-oriented surface traction g on a segment of length ∆ at the right bottom end of the domain. Table 1 reports geometrical and material parameters. If not differently specified, the initial conditions for the phase variable φ in the entire domain Ω are set as follows: φ 0 =v for all simulations relative to the formulation with volume constraint; φ 0 = 1 for all simulations relative to the formulation with volume minimization. Colour maps of the phase-field variable φ are represented in the region of Ω where the phase variable φ ≥ 0.1. Moreover, the simulation time is normalized with respect to the time constant T φ and, when not specified, a monolithic (SAND) solution strategy is adopted.  The design domain is discretized with (120 × 60) square elements along the two coordinate axes. As a result of a preliminary convergence analysis, the employed mesh density allows to obtain mesh-independent results for all cases under investigation. Table 2 reports the values of simulation parameters which are employed, if not object of a parametric investigation. It is noteworthy that adopted values respect the general rules-of-thumb traced in Section 4.1. In fact, the mesh size is h e = 0.016 m andē el ,φ = 80 MPa (with g = 500 MPa). The target final volume fraction with the volume minimization functional hence results v tar sol ≈ 0.4 with κ v = 100 MPa, that is comparable with the volume constraint valuev = 0.4 assigned for the volume constraint functional.

Staggered (NAND) versus monolithic (SAND) solution strategy: convergence behaviour
The first set of analyses compares the convergence behaviour of a staggered (NAND) versus a monolithic (SAND) solution strategy. This comparison is made only with reference to the classical volume constraint approach, i.e., functional L vc AC . Very similar results are however to be found for the novel volume minimization approach, i.e., functional L vm AC , as well. The evolution of the compliance C and of the Allen-Cahn error E AC versus the simulation time is shown in Fig. 3a. It can be observed that the error of the staggered approach decreases at the beginning of the simulation but starts to oscillate when reaching E AC ≈ 10 −3 , while the error of the monolithic strategy is substantially monotonic, down to the imposed value of 10 −6 .
The obtained final solution is then investigated for different values of the convergence parameter c conv AC , value with respect to which convergence of the Allen-Cahn procedure is checked. Figure 3b reports the volume fraction v sol and the structural displacement U sol at the converged final solution, together with the final interface perimeter P sol . As previously noted, the staggered approach does not return a converged solution when enforcing c conv AC ≤ 10 −3 . The obtained solutions obviously coincide in terms of final volume fraction (since enforced by means of a global constraint), but not in terms of structural displacement and interface perimeter. A converged final value is not obtained with the staggered approach, while a converged solution is obtained with the monolithic one for c conv AC ≤ 10 −3 , although c conv AC = 10 −2 would be practically sufficient. This outcome is confirmed by the maps on the distribution of the phase field variable at different values of c conv AC , reported in Fig. 4. Finally, Figure 3b shows also the comparison on the total number of Newton-Raphson iterations, confirming the significant out-performance of the monolithic solution strategy (SAND) versus the staggered one (NAND). Therefore, all following analyses will adopt a monolithic (SAND) approach.

Volume constraint versus volume minimization: convergence behaviour
This section presents a comparison between the convergence behaviour of the volume constraint versus the volume minimization formulation. For both formulations, two different initial conditions are investigated, that is φ 0 = 1 and φ 0 =v for L vc AC , and φ 0 = 1 and φ 0 = 0.5 for L vm AC . Figure 5a shows the evolution of the volume fraction v and of the compliance C along the simulation time. Contrarily to L vc AC , the volume fraction v with L vm AC evolves during the simulation, starting from the assigned initial condition. Structure compliance follows coherently. It is noteworthy that, for the chosen set of parameters, the final obtained solution between the two formulations is practically identical in terms of volume fraction and compliance, making the comparison consistent. Figure 5a reports also the Allen-Cahn convergence error E AC along the simulation time. In all cases, results show an acceptable convergence behaviour, with a similar decrease rate of E AC between the two formulations. However, the evolution of E AC for the volume constraint functional L vc AC is significantly more oscillatory than the one obtained for the volume minimization one L vm AC . The final solution is then investigated for different values of the convergence parameter c conv AC (with respect to which convergence of the Allen-Cahn procedure is checked). The final values of v sol and U sol are reported in Fig. 5b, together with P sol . Below a given threshold (ca. 10 −2 −10 −4 ), both formulations show a robust and converged solution. This outcome is confirmed by the maps on the distribution of the phase field variable at different values of c conv AC , reported in Fig. 6. Moreover, Fig. 5b reports the number of Newton-Raphson iterations required for the solution, clearly proving the computational out-performance of the volume minimization approach L vm AC with respect to the volume constraint one L vc AC . The effect of the initial condition for the volume constraint functional L vc AC is negligible both in terms of final solution and numerical performances. This is due to the fact that, at the very first step, the global constraint changes the value of φ from φ 0 tov. On the other hand, for the volume minimization functional L vm AC , the choice on φ 0 affects the solution. In fact, only minor differences are obtained in terms of final obtained solution (see also Fig. 6), but the convergence behaviour is significantly different, φ 0 = 1 requiring a significantly lower number of Newton-Raphson iterations than φ 0 = 0.5. This is also observable in Fig. 5a by noting that the adaptive solution scheme with φ 0 = 0.5 significantly reduces the reference time-step increment δt 0 = 10 −2 , while φ 0 = 1 not. Accordingly, in what follows, the initial condition φ 0 = 1 will be adopted for the volume minimization functional L vm AC , while φ 0 =v for L vc AC .

Volume constraint versus volume minimization: amount of distributed material
This section analyzes the two formulations by varying the amount of distributed material. This is achieved by varyingv for the volume constraint functional L vc AC , and κ v for the volume minimization functional L vm AC . For the latter, the value of the perimeter stiffness parameter κ φ varies with κ v according to the relationship κ φ = γ φ κ v , as described in Section 4.1.
As shown in Fig. 7, for L vc AC , the assigned valuev clearly corresponds to the obtained final volume fraction v sol for the entire range of the parametric analysis, verifying the correctness of the implementation of the global constraint. On the other hand, for L vc AC , the relationship between κ v and v sol is non-linear. Remarkably, the obtained trend is very-well captured by the theoretical target estimate v tar sol (κ v ) provided in Eq. (4.52). Figure 8 shows that, for a given value of final volume fraction v sol , the solutions obtained from the two formulations are practically identical in terms of structural displacement U sol , as well as maximum σ max vm and average σ avg vm Von-Mises stresses. The effective distribution of the phase-field variable is also very similar, although differences are observable in Fig. 9 and again in Fig. 8 as regards the interface perimeter P sol .
Finally, the analysis of the performance of the two formulations (in terms of Newton-Raphson iterations, see Fig. 8) show that the volume minimization functional L vm AC is significantly more efficient (> 50% iteration saving) than L vc AC for the more challenging cases cases where the final volume fraction is low, that is v sol < 0.5, cases which are also more interesting from the engineering viewpoint.

Volume constraint versus volume minimization: effect of load variations
The last comparison between the two formulations addresses the effect of the applied load magnitude g on the final solution. In this case, the volume constraintv (for L vc AC ) and the volume penalty κ v (for L vm AC ) are held constant, as given in Table 2. Figure 10 clearly show the difference on the rationale upon which the two functionals are built. In fact, the load magnitude g highly affects the final obtained structural displacement U sol when employing a volume constrained principle (i.e., L vc AC ) because the final volume fraction v sol is assigned a priori. This is accompanied by high variations in stresses, while the final design is practically independent from the load (see Fig. 11). On the other hand, the load magnitude g highly affects the final volume fraction v sol obtained when employing a volume minimization principle (i.e., L vm AC ). At different load levels, final designs associated with similar values of structural displacement U sol and stresses σ max vm and σ avg vm are obtained. As shown in Fig. 11, the obtained design with L vm AC is now highly affected by the applied load value g.   (a) Volume fraction v(φ) (left), compliance C(φ, u) (center) and Allen-Cahn error E AC (right).   3.0

3D simulations
In this section, the novel volume minimization functional L vm AC is tested on two 3D problems. In particular, the addressed examples are (see Fig. 2): -A cantilever beam, representing the 3D counterpart of the case study presented in Section 5.1.
The domain is now discretized by means of 80 × 40 × 40 hexahedral elements, resulting in an average element size h e = 0.025 m. The geometrical and material features, as well as the applied load magnitude, correspond to the values listed in Table 1, while values of simulation parameters to the ones in Table 2 apart from γ φ = 0.02 m (due to a slight increase in the dimension of the elements in the mesh), and c conv AC = 10 −2 (fixed on the basis of the convergence analysis in Section 5.1.2); -A bridge, as the one investigated in Zegard and Paulino (2016). The design domain is represented by a rectangular cuboid x × y × z . The domain is fixed on the bottom plane at Z = 0 along strips of length l bc along X. Moreover, a passive solid slab on the top surface at Z = z of height h s m is considered, on top of which a distributed load g is applied in the Z-direction.
Material properties correspond to the one of the previous example, being reported in Table 1.
Values of geometrical and simulation parameters, as well as load magnitude, are given in  In particular, the 3D cantilever has been addressed to show that the proposed formulation can be straightforwardly implemented in a 3D computational framework, analyzing differences obtained from a 2D to a 3D setting for the same case study. The bridge is instead defined such to have significant differences in the physical dimension of the problem, in order to prove the robustness of the parameter settings guidelines traced in Section 4.1, with particular reference to κ v , γ φ , and κ φ .
Since the volume minimization functional is adopted, the setting of the volume penalty parameter plays a fundamental role. For the cantilever, we haveē el φ = 84 MPa and hence v tar sol ≈ 0.4 with κ v = 100 MPa; for the bridge, we haveē el φ = 170 MPa and hence v tar sol ≈ 0.2 with κ v = 1000 MPa. Figure 12 shows the comparison between the estimated and final obtained volume fractions, confirming the effectiveness of the proposed procedure for obtaining the target solution.
The final obtained designs are shown in Fig. 13 for the cantilever, and in Fig. 14 for the bridge. Three-dimensional patterns clearly arise, with regular external and internal surfaces obtained without the need of post-processing filtering techniques thanks to the adopted phase-field rationale.

Comparison between investigated formulations and conclusions
The present work addresses a rigorous study on basic properties of structural topology optimization problems with phase-field. In particular, we performed the following five tasks: 1. Proposed mixed Hu-Washizu variational formulations for the topology optimization problem with phase-field to directly impose equilibrium, constitutive, and compatibility equations in the formulation; 2. Investigated two different topology optimization principles, e.g.,: i) a formulation imposing a priori the amount of material to be distributed within the design domain (formulation with volume constraint); ii) a formulation based on a minimization of material to be distributed, given that a cost (i.e., a penalty parameter) is assigned to the material (formulation with volume minimization); 3. Introduced a Simultaneous Analysis and Design (SAND) monolithic solution strategy, thanks to the Hu-Washizu functional rationale and based on an Allen-Cahn scheme, where the phase-field variable evolves under the respect of mechanical equilibrium at each computational incremental step; 4. Analyzed both numerical convergence behaviours and obtained final designs, based on comparative analyses between simulation strategies (monolithic SAND vs. staggered NAND) and between topology optimization principles (volume constraint vs. minimization); 5. Analyzed the performance of the proposed variational formulation based on volume minimization also on three-dimensional case studies.
On the basis of theoretical considerations, general guidelines are traced for properly setting the value of simulation parameters characterizing the phase-field evolution (i.e., τ φ , κ φ and γ φ ). Following this careful setting, the implemented finite element formulations and solution algorithm are robust, allowing to conduct a wide campaign of parametric simulations without encountering any numerical convergence issue. In addition, the volume penalty parameter κ v in the volume minimization functional has been related to target values of volume fraction. From the obtained results, we may conclude the following: -A monolithic SAND solution strategy significantly outperforms a staggered NAND solution strategy (see Fig. 3); -The volume minimization formulation L vm shows an Allen-Cahn convergence behavior with better properties than the volume constraint formulation L vc , since the decrease of the Allen-Cahn error measure E AC is highly oscillatory for L vc but not for L vm (see Fig. 5); -Comparing designs with the same final volume fraction v sol , the two formulations lead to practically coincident solutions in terms of compliance C sol (or structural displacement U sol ) and interface perimeter P sol (see Fig. 8), although final designs are slightly different (see Fig.  9); -The final design is practically independent from the applied load with the volume constraint functional L vc , while is highly affected with the volume minimization one L vm (see Figs. 10 and 11), which is seen as a significative advantage of this latter formulation; -The number of Newton-Raphson iterations required for solving the problem with the the volume minimization functional L vm is in most cases significantly lower (less than 50%) than the ones employed with the volume constraint functional L vc ; -No convergence issues are encountered in 3D applications (see Figs. 13 and 14), allowing to obtain (thanks to the phase-field rationale) complex but regular patterns without the need of any post-processing filtering technique (as required by density-based approaches) which might lead to undesirable (artificial) effects (Zegard and Paulino, 2016).
In conclusion, the obtained results highlight fundamental aspects of structural topology optimization problems with phase-field, allowing to optimize the solution strategy and to trace guidelines for settings simulation parameters beforehand. Referring in particular to the proposed volume minimization principle, we believe that these results can be a starting point for more advanced developments of phase-field topology optimization that considers loading uncertainties (Dunning et al., 2011) or multi-target strategies, e.g., controlling both geometry and compliance (Strömberg, 2010), both geometry and stresses (Burger and Stainko, 2006), the structure life-cycle cost (Sarma and Adeli, 2002), or manufacturing costs (Liu et al., 2019).

A Interpolation of the stress and strain fields
The interpolation matrices of the stress and strain fields employed for 2D and 3D applications are detailed in what follows. Choices follow arguments presented in Weisman (1996), Cao et al. (2002) and Djoko et al. (2006). In detail, minimum distributions required for stability are included as a common basis between the assumed stress and strain fields. In addition, the assumed strain field is enriched with respect to the stress one with strain modes not already contained in the minimum one. A minimal strain enrichment, a priori orthogonal to the assumed stress field, is considered. Interesting extensions of the present work would be to investigate more refined enrichment strategies, accompanied by a generalization of the introduced mixed formulations, in order to deal with material constrained behaviours (e.g., incompressibility or inextensibility).