A computational framework for topology optimisation of flexoelectricity at finite strains considering a multi-field micromorphic approach

This paper presents a novel in-silico framework for the design of flexoelectric energy harvesters at finite strains using topology optimisation. The main ingredients of this work can be summarised as follows: (i) a micromorphic continuum approach is exploited to account for size dependent effects in the context of finite strains, thus permitting the modelling and simulation of flexoelectric effects in highly deformable materials such as dielectric elastomers. A key feature of the multi-field (mixed) formulation pursued is its flexibility as it permits, upon suitable selection of material parameters, to degenerate into other families of high order gradient theories such as flexoelectric gradient elasticity. (ii) A novel energy interpolation scheme is put forward, whereby different interpolation strategies are proposed for the various contributions that the free energy density function is decomposed into. This has enabled to circumvent numerical artifacts associated with fictitious high flexoelectric effects observed in the vicinity of low and intermediate density regions, where extremely high strain gradients tend to develop. (iii) A weighted combination of efficiency-based measures and aggregation functions of the stress is proposed to remedy the shortcomings of state-of-the-art efficiency-based functionals, which promotes the development of hinges with unpractical highly localised large strain gradients. Finally, a series of numerical examples are analysed, studying the development of direct flexoelectricity induced by bending, compression and torsional deformations


Introduction
Flexoelectricity [1,2] is a two-way coupling mechanism between electric polarisation and inhomogeneous deformations (i.e., bending, torsion), only relevant at small length scales.To be more precise, direct flexoelectricity is responsible for the generation of polarisation in the material due to strain gradients.On the other hand, converse flexoelectricity accounts for the development of stresses within the material due to gradients of the electric field.Recent advancements in fabrication and design of miniaturised components have indeed motivated the interest of the scientific and engineering community for flexoelectric materials, leading to a vast number of publications in the experimental characterisation of this phenomenon [3][4][5][6].
Some of these works study the flexoelectric effect (at small length scales) inherent in piezoelectric crystals, such as Barium Titanate, devoid from a centrosymmetric internal structure.In these materials, it is interesting to quantify the electric energy originated either from piezoelectricity or flexoelectricity, when subjected to a given mechanical loading [5,7].However, flexoelectricity is not exclusive from non-centrosymmetric materials, as shown in Refs.[8], demonstrating flexoelectricity in perovskite ferroelectrics [6], graphene nano-shells and biological membranes [3,4], all of them centrosymmetric.This is achievable when the inversion symmetry of these materials is broken through the application of non-uniform strains.Even dielectric elastomers, which are isotropic, and hence centrosymmetric, can exhibit flexoelectricity at small length scales [8].The latter, namely dielectric elastomers, open up for new possibilities in the field of flexoelectric-based nano-actuators, sensors and energy harvesters, taking advantage of their capability of exhibiting large deformations [9][10][11][12].At the macroscopic level, their ability to undergo large electrically induced deformations, in conjunction with their fast response, led to their identification as ideal candidates in the field of soft robotics.For instance, electrically induced area expansions of over 380% on dielectric elastomer thin films placed on the verge of snap-through configurations have been reported by Li et al. [13].Using a theoretical framework, Ref. [14] demonstrates also the capability of dielectric elastomers balloons to undergo snap-through electrically induced actuations.Other applications for dielectric elastomers include Braille displays, deformable lenses, haptic devices and energy generators, to name but a few.
Flexoelectricity represents a generalisation of the simpler electro-elasticity continuum theory [15][16][17][18].The variational formulation of the latter is extremely well established.In the most standard formulation, displacements and the scalar electric potential [15,19,20] are modelled as the unknown fields.In this formulation, the constitutive information is encapsulated in the Helmholtz energy functional via its invariant-based representation depending upon kinematic strain measures and the electric field [21][22][23].However, for more complex constitutive models than that of an ideal dielectric elastomer, the saddle point nature of the Helmholtz functional (convex with respect to the deformation gradient tensor and concave with respect to the electric field in the small strain/small electric field regime), makes in general impossible to define a priori constitutive models which satisfy the ellipticity condition [24,25], which ensures the well-posedness of the problem.Motivated by the possible loss of ellipticity of the Helmholtz functional, Gil and Ortigosa [26][27][28][29] extended of the concept of polyconvexity [25,[30][31][32][33][34][35][36][37][38] to the field of electro-mechanics by postulating a definition of the internal energy functional convex with respect to an extended set of arguments, proving that this definition satisfied the ellipticity condition unconditionally.
The previous classical continuum theory of electro-mechanics cannot be applied for the modelling of flexoelectricity.Instead, high order gradient theories and generalised continua [39] need to be used in order to account for size dependent effects inherent to flexoelectricity.For the particular case of dielectric elastomers, in the context of finite strains, there exist two approaches relying on two different types of high order continuum theories.On one hand, the works in [40][41][42][43] rely on continuum frameworks that lie within the classification of (flexoelectric) gradient elasticity.From the numerical point of view, the authors in [40,41] resort to different types of C 1 finite element approximations of the displacement field in order to handle the higher regularity entailed by the fourth order derivatives featuring in the formulation.However, the requirement of C 1 continuity in the displacement field approximation can be eliminated when mixed-type finite element formulations are employed.In this regard, see the work in [42,43], the latter using a couple stress-based approach for flexoelectricity.On the other hand, a more general approach for the modelling of flexoelectric materials is that pursued in [44], where a (flexoelectric) micromorphic approach in conjunction with a mixed-type finite element formulation is advocated for.The great advantage of the latter resides on its flexibility, in the sense that a suitable tuning of the material parameters featuring in the constitutive model permits to degenerate the formulation into that of flexoelectric gradient elasticity succinctly described previously.
Mathematical tools such as Topology Optimisation (TO) represent a powerful instrument with the potential to assist experienced researchers in the design of smart or stimuli responsive materials with complex underlying physical behaviours.The design of flexoelectric materials can indeed benefit from the use of such tools.TO of smart materials has been widely explored mainly in the context of linearised elasticity, making use of the wide spectrum of robust approaches available for TO, ranging from density-based methods, with the Solid Isotropic Material with Penalisation (SIMP) method as their maximum representative [45], level-set methods [46,47], phasefield methods [48], topological derivative methods [49] and evolutionary methods [50].Kang and Wang [51] used the SIMP method for the TO of piezoelectric ceramics, an electro-active material which unlike DEs, is restricted to scenarios characterised by very small deformations and/or displacements.Zhang et al. [52] also investigated the TO of piezoelectric sensors with active vibration control purposes.Other works [53][54][55][56] have investigated the simultaneous optimisation of polarisation and layout of piezoelectric ceramics over a fixed host (passive) material, and when the host structure is also included in the optimisation process [57][58][59].Furthermore, the work in [60] carries out the TO of piezoelectric ceramic-based micro-grippers considering large displacements but small deformations.Relevant works in the field of TO of piezoelectric composites include [61][62][63].Some recent work includes the design by means of TO of thermoelectric coolers and generators [64][65][66].More recently, some works have ventured in the TO of dielectric elastomers.For instance, the work in Reference [67] delves into the TO of these materials with the aim of maximising the electrically induced rotation of a rotary device.Extremely interesting is the work in [68], where the TO is applied with the objective of conceiving intelligent microstructures for the design of wide tunable band gaps.Although not using a TO framework, the tunability of band gaps in the context of dielectric elastomers has also been investigated in a vast number of works (see for instance [69][70][71]) Furthermore, the recent works by Ortigosa et al. [72][73][74] presented a thorough numerical study on density-based and phase field-based TO in the context of nonlinear electro-mechanics [15][16][17]75,76], considering both DEs and piezoelectric polymers.That work proposes a novel energy interpolation scheme, where suggestions regarding the mechanical and electrical exponents featuring in the proposed interpolation scheme are made in order to yield designs devoid from undesired intermediate density regions.
Topology optimisation can help designing efficient energy harvesters taking advantage of the flexoelectric effect at small length scales.Few works exist on the effective application of TO techniques for the design of energy harvesters based on the principle of flexoelectricity [77,78].The latter works are however restricted to the case of small strains.As shown in [77,78], embedding TO algorithms with the aim of maximising an efficiency-based objective function, typically measuring the ratio between the global (integral) mechanical energy and the global electric energy generated through flexoelectricity, leads to designs characterised by a poor structural performance, with hinges where high stress concentration occurs.This mechanical deficiency is indeed promoted by the inherent flexoelectric nature of the material, which promotes the development of hinges with highly localised large strain gradients.The latter leads to an extremely rapid growth of the electric energy generated when compared with the associated mechanical energy decay, yielding extremely efficient designs according to the metric chosen, but inappropriate from a structural/mechanical standpoint.In order to remedy the above shortcoming, this paper presents a new TO framework for the robust (safe) design of flexoelectric actuators with the following key novelties: (a) consideration of dielectric elastomers, hence, of finite strain scenarios; (b) definition of objective functions considering a weighted combination of efficiency-based measures and p-norm/KS aggregation/regularisation [79] functions of the Von Mises stress distribution within the material.
The layout of this paper is as follows: Section 2 introduces the flexoelectric micromorphic elasticity formulation presented in Ref. [44] in a self-contained manner.Section 3 describes all the ingredients of the optimisation problem by means of the SIMP method.Finally, Section 4 shows some numerical examples demonstrating the applicability of the proposed formulation.Finally, Section 5 provides some concluding remarks.

Flexoelectric micromorphic elasticity
An introduction into nonlinear continuum electromechanics considering a micromorphic approach will be covered in this section, as well as the relevant governing equations.

Kinematics in the macro-scale
Let us consider the motion of a continuum with reference configuration B 0 ⊂ R n (n = {2, 3}).After the motion, the solid occupies a deformed configuration B ⊂ R n (refer to Fig. 1).We assume the existence of a deformation mapping φ (X) linking material particles from the reference configuration X ∈ B 0 to the deformed configuration x ∈ B according to x = φ (X).Associated with φ (X), the deformation gradient tensor F [80][81][82][83] is defined as which permits to define its co-factor and its Jacobian as

The micromorphic problem
In this paper, a micromorphic approach is considered in order to mathematically describe the continuum B 0 .Hence, attached to every macroscopic material point X ∈ B 0 , there exists a microelement B X 0 , centred in X. Microscopic points in B X 0 are denoted as Ξ ∈ B X 0 .Similarly, attached to every macroscopic spatial point x ∈ B, there exists a microelement B x , centred in x.Microscopic points in B x are denoted as ξ ∈ B x .The latter are related with their material counterparts Ξ ∈ B X 0 through the following affine transformation [39] where G ∈ R n×n represents the microscopic deformation gradient tensor.Notice that in the micromorphic continuum, both microscopic and macroscopic deformation gradient tensors are independent fields, i.e.F ̸ = G.Associated with G, it is possible to define its gradient with respect to the macroscale material coordinates X ∈ B 0 as

Governing equations
In the absence of inertia effects, the local form of the conservation of linear momentum [83] in the solid B 0 can be written as where f 0 represents a body force per unit undeformed volume B 0 and t 0 , the traction force per unit undeformed area on ∂ t 0 B 0 .Furthermore, φ represents the value of the Dirichlet boundary condition on ∂ φ B 0 , and with N the outward unit vector on the boundary, and Finally, P represents the macroscopic first Piola-Kirchhoff stress tensor.With regard to the equations for electrostatics in the macroscale, these reduce to Gauss's and Faraday's laws on B 0 .The local form of the Gauss' law [76,84] can be written in a Lagrangian formalism as where D 0 is the Lagrangian electric displacement vector, ρ 0 represents an electric volume charge per unit of undeformed volume B 0 and ω 0 , an electric surface charge per unit of undeformed area ∂ ω 0 B 0 ⊂ ∂B 0 .Furthermore, under the same hypothesis, the local form of the Faraday's law can be written in a Lagrangian setting as where E 0 is the Lagrangian electric field vector and ϕ, the scalar electric potential.In (7), ∂ ϕ B 0 represents the part of the boundary ∂B 0 where essential electric potential boundary conditions are applied such that ∂ ω 0 B 0 ∪∂ ϕ B 0 = ∂B 0 and ∂ ω 0 B 0 ∩ ∂ ϕ B 0 = ∅.The system of coupled PDEs is completed with the governing equations associated with the micromorphic continuum, which are expressed in a Lagrangian setting as [44] DIVΣ where Σ G ∈ R 3×3×3 and Σ G ∈ R 3×3 represent the micromorphic double stress and the micromorphic Piola stress, respectively.Furthermore, T 0 ∈ R 3×3 represents the micromorphic Piola traction on ∂ T 0 B 0 and Ḡ, the prescribed value of G on Remark 1.As indicated in Steinman et al. [44], the physical meaning of Dirichlet boundary conditions on the field G is not completely understood.Nonetheless, this type of boundary conditions have been included in (8) for the sake of completeness.However, in the numerical examples in Section 4, only natural (or Neumann) boundary conditions are considered for this field.

Variational formulation
In this section we will present the variational formulation associated with the nonlinear PDEs in Eqs. ( 5), ( 6), ( 7) and (8), which stems from a three-field variational principle with unknown fields U = {φ, ϕ, G} belonging to suitable functional spaces where H 1 (B 0 ) stands for the Sobolev space of functions belonging to L 2 (B 0 ) (square summable functions) and with gradient also in L 2 .The three-field variational principle, denoted as Π (U ), can be defined as where Ψ (F ) represents the free energy density of the material.In above equation, it is implicitly assumed that {F, E 0 , G} depend upon the gradients of the unknown fields U = {φ, ϕ, G} through Eqs. ( 1), ( 7) and (4), respectively.In addition, the external contribution Π ext (φ, ϕ, G) is decomposed into volumetric and surface contributions, i.e. with The directional derivatives of Π (U ) in (10) with respect to the weak forms of ( 5), ( 6) and (8), respectively, namely Integration by parts an application of the Gauss' theorem on Eqs.(13), permits to conclude that these correspond with the weak forms of the boundary value problem associated with Eqs. ( 5), ( 6), (7) and (8), provided that { P, D 0 , Σ G , Σ G } are related with the derivatives of the free energy density Ψ (F, E 0 , G, G) featuring in (13) through the following relationships The stationary conditions in (A.3), namely DΠ [ p φ ] = 0, DΠ [ p ϕ ] = 0 and DΠ [ p G ] = 0, are coupled through the nonlinearity of the free energy density Ψ (F, E 0 , G, G).It is customary to make use of a Newton-Raphson scheme where the three weak forms can be linearised with respect to incremental fields {∆φ, ∆ϕ, ∆G} ∈ which permit the update of {φ, ϕ, G} ∈ V φ h × V ϕ h × V G h at a given Newton-Raphson iteration k + 1 as In (15), the three terms D 2 Π [ p φ ; ∆φ], D 2 Π [ p φ ; ∆ϕ] and D 2 Π [ p φ ; ∆G], stemming from the linearisation of DΠ [ p φ ] with respect to ∆U , are obtained as Furthermore, in (15), the three terms stemming from the linearisation of DΠ [ p ϕ ] with respect to ∆U , are obtained as Finally, linearisation of DΠ [ p G ] with respect to ∆U yields The details concerning the Finite Element discretisation of the nonlinear stationary conditions in (13) and their linearisations in ( 17)-( 19) can be found in Appendix A.

Constitutive equations
The free energy density of the material Ψ (F ) (with F = {F, E 0 , G, G}) was introduced in the three-field variational principle Π (U ) in Eq. ( 10), with As shown in Eq. ( 14), the derivatives of Ψ (F ) with respect to F permit to obtain their respective work conjugates, namely { P, D 0 , Σ G , Σ G }.The objective of this section is to specify the constitutive model followed in this paper, namely, the explicit expression for the free energy density Ψ (F ).In order to do so, we anticipate that the have followed the same model as that considered in Ref. [44], where the free energy density is additively decomposed into four contributions as where Ψ mac (F, E 0 ) stands for the macroscopic contribution, Ψ mic (G) refers to the micromorphic contribution, Ψ scale (F, G) is the scale-bridging contribution and finally, Ψ flexo (E 0 , G, G) is the flexoelectric contribution, being responsible for the flexoelectric effect.In the subsequent sections, each of the contributions of the free energy density will be specified.

The macroscopic free energy density
The macroscopic energetic contribution, namely Ψ mac (F, E 0 ), corresponds with the free energy density function of a dielectric elastomer, capable of exhibiting large electrically induced deformations.As customary in nonlinear electro-mechanics, the energy functional Ψ mac (F, E 0 ) is additively decomposed into a purely mechanical contribution and a coupled contribution, denoted as Ψ macro m and Ψ macro em , respectively, as For the mechanical contribution Ψ mac m (F), polyconvex [25] free energy density function will be considered, namely any function that can be re-written as a convex multivariable function in terms of the extended set of kinematic entities V = {F, H, J }, i.e.
. Specifically, in this paper we consider the Mooney-Rivlin free energy function to describe the macroscopic contribution of the free energy Ψ (F ), which is defined as with is clearly convex with respect to V, provided that the material parameters {µ 1 , µ 2 , λ} are positive.Notice that in (24), || • || represents the Frobenius norm of the argument •.These can be related with the first and second Lamé parameters in the origin (i.e.F = I, with I ∈ R 3×3 the second order identity matrix), denoted as {λ 0 , µ 0 }, respectively, according to the following expressions For the electro-mechanical contribution Ψ mac em (F, E 0 ), we consider the free energy function of an ideal dielectric elastomer, which is defined as where ε 0 represents the electric permittivity of vacuum, being ε 0 = 8.854 × 10 −12 F m −1 , and ε r the relative electric permittivity.
Remark 2. Notice that Ψ mac em (F, E 0 ) is the dual of the following internal energy density function, where the relationship between both functions is given by the following Legendre transformation where ( 27) can be equivalently written as The authors (see Refs. [26][27][28]) have shown that the definition of W mac em (J, d) is indeed convex with respect to both of its arguments, provided that J > 0 (which is a physical condition preventing inter-penetrability of matter), namely the Hessian ] if positive definite.In Refs.[26][27][28], the authors postulated the extension of polyconvexity to the field of electromechanics by permitting the internal energy density e(F, D 0 ) to be convex with respect to any of the elements within the set {F, H, J, D 0 , d}.The use of the free energy Ψ mac em (F, E 0 ) in ( 26) may give the impression that it can jeopardise the polyconvexity condition, through the use of the minus sign.However, re-expressing Ψ mac em (F, E 0 ) in terms of its dual e mac em (F, D 0 ) dissipates any doubt.Therefore, Ψ mac em (F, E 0 ) is indeed polyconvex in the context of electromechanics in the sense of [26][27][28].

Micromorphic free energy density
Following the work in [44], the expression adopted in this work for the micromorphic contribution Ψ mic (G) is the following where l ≥ 0 is the length-scale parameter and with µ 0 in Eq. (25).Notice that ||G|| 2 in (28) represents the Frobenius norm of the third order tensor G, and its expression has been shown in (28), where the use of repeated indices implies summation.

Scale-bridging free energy density
For the scale-bridging contribution Ψ scale (F, G), we use the same expression as that followed in Ref. [44], namely where the material parameter κ in ( 29) is a penalty-like parameter coupling both macro-and micro-deformation gradients, i.e.F and G, respectively.

Flexoelectric free energy density
Finally, the flexoelectric contribution Ψ flexo (E 0 , G, G), according to Ref. [44] adopts the following expression where the flexoelectric coefficient ν f is defined as where ν f represents the dimensionless flexoelectric coefficient.
2.5.5.Specific form of the work conjugates and second order derivatives of the free energy density The specification, through Sections 2.5.1-2.5.4, of each of the energetic contributions that the free energy density Ψ (F ) is decomposed into, permits to give the explicit definition of the work conjugates { P, D 0 , Σ G , Σ G } according to (14), namely where the cross product operation between second order tensors is defined as with A ⊗ B denoting the dyadic product for any nth order tensor A and B (n ∈ N) and with E the third order alternating tensor.In addition, we show in this section the various second partial derivatives of the free energy, namely ∂ 2 XY Ψ (F ) (for any X, Y ∈ F), featuring in the different contributions of the linearised expressions in ( 17)- (19).In particular, for the specific choice of the constitutive model described in Section 2.5, these are as follows +κI; with ∀ A, B ∈ R 3×3 and A ∈ R 3×3×3×3 .In addition, the derivatives Furthermore, the derivatives ∂ 2 GG Ψ (F ) and ∂ 2 GG Ψ (F ) can be obtained as with Finally, the last derivative featuring in (A.6) is 2.6.The flexoelectric micromorphic approach in the context of microcontinuum theories Fig. 2 showcases the convenience and flexibility of the micromorphic approach adopted in this manuscript.As discussed in Ref. [44], one of the great advantages of the micromorphic approach considered, denoted as FM-Elasticity (which stands for coupled Flexoelectricity and Micromorphic Elasticity), resides in the fact that it can be particularised, by tuning some of the material parameters featuring in the different energy contributions described in the preceding sections, into a wide range of less generic modelling approaches.For instance, for the limiting case when the scale bridging coefficient κ in ( 29) is κ → ∞, the formulation reduces to Flexoelectric Gradient Elasticity (FG-Elasticity) [41], where the micro-and macroscopic deformation gradient tensors G and F coincide, namely F = G, and hence, G = ∇ 0 F. FG-Elasticity permits the use of mixed-type variational formulations, but it can also be modelled through a two field variational principle (with only {φ, ϕ} as unknown fields), forcing a higher regularity for the discretisation of the field φ, requiring C 1 (B 0 )-continuous finite element basis functions.The latter is precisely the formulation pursued in [40,41].
Furthermore, setting the flexoelectric coefficient ν f = 0 in ( 30)-( 31), FM-Elasticity reduces to a family of continuum formulations devoid from flexoelectricity.The most generic of them is EM-Elasticity (i.e.Electro Micromorphic Elasticity).The latter can reduce into other less generic continuum approaches such as Electromechanics Gradient Elasticity (EG-Elasticity) and Electromechanics Elasticity (E-Elasticity), and into formulations excluding at all any source of electromechanical coupling, such as Micromorphic Elasticity (M-Elasticity), Gradient Elasticity (G-Elasticity), and finally, Elasticity.

Introduction
As customary in density-based topology optimisation (see [45]), we seek the optimised design, in this specific case, of a flexoelectric energy harvester, through the introduction of the design density field ρ(X) ∈ [0, 1], which represents a continuous scalar field introduced in order to differentiate solid from void regions in B 0 .In order to prevent spurious results, we make use of standard filtering techniques [85,86] such as the conic filter, which can be stated as the convolution product where ∆ : R + → R + is the so-called convolution kernel, ρ(X) is referred to as the filtered density field, and with ∆(r ) = max{0, 1 − r/R}, being R the filter radius.Alternatively, it is possible to use the Helmholtz' filter, resulting from the solution of the following PDE where R represents the filter radius and where ∆ 0 (•) represents the material or Lagrangian Laplacian operator.In addition, we use the classical smoothed Heaviside projection function proposed in [87] ρ(X) = where ρ(X) is known as the physical density field and β and η are parameters carefully selected and updated throughout the optimisation [87].The field ρ : B 0 → [0, 1] is used in order to interpolate the free energy density of the body Ψ (X, F ), which depends upon each macroscopic point X ∈ B 0 through ρ(X), namely The additive decomposition of the free energy considered, shown in Eqs. ( 21) and ( 22), permits to establish an independent energy interpolation scheme for each of the five energetic contributions considered, namely The functional dependence of each of the five contributions in Eq. ( 43) with respect to the field ρ(X), is dictated by the choice of the energy interpolation scheme considered.Irrespectively of that choice, the five energy interpolation schemes in (43) must verify that, for the limiting case ρ = 1, the interpolated energy must coincide with that of the solid material considered.On the other hand, for the limiting case ρ = 0, the interpolated energetic contributions are forced to coincide with a small fraction (controlled by a parameter α ≪ 1) of that of the solid material (a zero contribution is avoided in order to prevent numerical stability issues).For the macroscopic electromechanical contribution Ψ mac em , however, we impose that for the limiting case ρ = 0, the energetic contribution must coincide with that of vacuum, denoted as Ψ vac em , and defined as All of that can be mathematically stated below as Unfortunately, the limiting values of Ψ mac m ( ρ(X), F) and Ψ mac em ( ρ(X), F, E 0 ) in ( 45) for ρ = 0 can induce artificial numerical instabilities.In order to remedy this shortcoming, it is customary [72,73] to replace them with the following more stable definitions where it can be seen that the electromechanical contribution in the void region Ψ void em (F, E 0 ) is indeed a particularisation of that of the vacuum Ψ vac em (F, E 0 ) for the case H = I.Inspired by the SIMP method [45], a possible energy interpolation scheme for each of the five energy contributions that complies with the extreme or limiting cases established in (45) and ( 46) is the following where the five parameters { p m , p e , p mic , p scale , p flexo } featuring in the five interpolated energy contributions in (47), , play the role of the penalising exponent in the SIMP method in the context of elasticity, where only one exponent is needed.

The choice of the objective function
In this paper, we seek to obtain the optimised design of a flexoelectro energy harvester B 0 .Specifically, we aim at optimising an objective function embedding the energy harvesting properties of the flexoelectric material.With that objective in mind, previous references in the field of topology optimisation of flexoelectric materials [77,78], identify the objective function as the efficiency, denoted as J η ( ρ, φ, ϕ), and defined as the (negative) ratio between the mechanical and electromechanical contributions of the macroscopic energy contribution, namely However, as it will be shown in the numerical examples section, identification of the efficiency in (48) as the objective function of the problem leads to designs characterised by a poor structural performance, with hinges where high stress concentration occurs.This mechanical deficiency is indeed promoted by the inherent flexoelectric nature of the material, which promotes the development of hinges with highly localised large strain gradients.The latter leads to an extremely rapid growth of the electric energy generated when compared with the associated mechanical energy decay, yielding extremely efficient designs according to the metric chosen, but inappropriate from a structural/mechanical standpoint.In order to remedy the above shortcoming, it is possible to propose an objective function defined as a combination of the efficiency in (48) and a p-norm aggregation/regularisation function of the Von Mises stress distribution within the material. 1The latter, denoted herein as J τ ( ρ, φ), can be defined as where the quality of the approximation of the maximum operator is clearly affected by the choice of the exponent γ ∈ R + .Typically, whilst high values of γ yield a better approximation, they might compromise the stability of the optimisation algorithm due to the high nonlinearity of the regularised functional J τ ( ρ, φ).Furthermore, the local stress measure τ ( ρ(X), φ) in ( 49) is defined as where τ dev ( ρ(X), φ) is the deviatoric part of the interpolated Kirchhoff stress tensor τ ( ρ(X), φ), namely where the fourth order tensor P dev is defined as where δ i j refers to the −i j component of the Kronecker delta tensor.Notice in above Eq.( 51) that the interpolated energy Ψ mac m ( ρ(X), F) is defined similarly to the mechanical contribution of the macroscopic interpolated energy In fact, the only difference between Ψ mac m ( ρ(X), F) in ( 53) and Ψ mac m ( ρ(X), F) in (47) resides in the choice of the penalising exponent, which is either p τ or p m , respectively.This difference is motivated by the possible development of classical shortcomings, i.e. stress singularity, associated with stress minimisation or stress constrained topology optimisation [79].In order to remedy these, it is customary to consider a value of the penalising exponent of p τ = 0.5 [79].Having introduced both functionals J η ( ρ, φ, ϕ) in ( 48) and J τ ( ρ, φ) in (49), it is now possible to define the objective function J ( ρ, φ, ϕ) as a weighted combination of the two according to where J η 0 and J τ 0 refer to initial values of J η and J τ , respectively, i.e. their values at the initial topology optimisation iteration.In Eq. ( 54), ω represents the weight associated with the stress-based objective functional J τ ( ρ, φ).Alternative definitions for the objective function J ( ρ, φ, ϕ) can be considered, as the one presented below Remark 3. The Von Mises-based functional J τ in (49) has been introduced with the aim of reducing high stress values which, in turn, prevents the development of hinges and hence, of extremely localised deformations.This, therefore, can help limiting the value of the electric field attained through the inherent flexoelectric behaviour of the material.As a result, a reduction of high electric field concentrations is achieved in an indirect manner.Alternative functionals to J τ could have been defined incorporating as a constraint, or in a multi-objective fashion, the maximum value of the norm of the electric field attained through flexoelectricity.Through this approach, a direct control on the maximum electric field can also be obtained.Either of above two approaches would permit to incorporate important limiting physical phenomena inherent to dielectric elastomers, such as electrical breakdown [88], into the optimisation problem.
Clearly, in order to obtain the descend direction DL X (ρ, U , p U )[∆ρ], it is necessary first to obtain the values of the unknown fields U and those of the adjoint states p U .The first, namely the unknown fields U , are obtained from the optimality condition of the Lagrangian L X (ρ, U , p U ) with respect to δ p U (alternatively called as the forward problem), namely where the expression for DΠ (ρ, U )[δ p U ] are the same as those for DΠ (ρ, U )[ p U ] in (61), exchanging δ p U with p U .Clearly, the optimality conditions for both Lagrangians are exactly identical, namely and therefore, the forward problem needs to be solved only once.The adjoint fields p U can be obtained from the optimality condition of the Lagrangian L X (ρ, U , p U ) with respect to δU (alternatively called as the adjoint problem), namely where we have made use of the Schwartz theorem that permits to use ; δU ], due to their equivalence.Notice in above equation that the expressions for D 2 Π (ρ, U )[δU ; p U ] in (63) are in fact the same as the linearised expressions in Eqs. ( 17)- (19).Finally, the explicit expressions for DJ X ( ρ, φ, ϕ)[δφ] for X = {η, τ } will be shown in Sections 3.3.1 and 3.3.2.The different expressions for both cases where X = {η, τ } entails that the adjoints states p U corresponding to both Lagrangian functionals L X ( ρ, φ, ϕ) are indeed different, and therefore, Eq. ( 63) needs to be solved twice, one for each case, namely for X = η and for X = τ .Finally, the Finite Element implementation of both adjoint problems is shown Appendix B.

Algorithmic implementation
Finally, the reader interested in the algorithmic details of the optimisation algorithm implemented in our inhouse platform, can check the main ingredients necessary for computational implementation purposes in Algorithm 1, presented in a pseudo-code format.

Numerical examples
The objective of this section is to reveal the performance of the density-based topology optimisation approach described throughout Section 3. Specifically, the goal is to verify the capability of the proposed approach in the context of the Flexoelectric Micromorphic Elasticity (FM-Elasticity) formulation described through Section 2.
Prior to that, we describe the common features of all the numerical examples contained in the forthcoming subsections.First, we comment on the choice of the various penalising exponents { p m , p e , p mic , p scale , p flexo } featuring in the proposed interpolation scheme in Eq. ( 47).The reader can observe the values for these in Table 1.Notice from that table that p m = p e = 3.This choice is motivated by previous experience from the authors in the field of topology optimisation in nonlinear continuum electromechanics [72].The choice of the other two parameters, namely { p mic , p scale }, is motivated by our empirical observations.With regard to the penalising scheme in the flexoelectric energy contribution, namely p flexo , this deserves a careful discussion.After a detailed empirical analysis (i.e.numerical simulations), we observed that a considerably larger value was needed in comparison to the other four penalising exponents.The underlying reason for this difference in values lies in the fact that high strain gradients tend to occur in areas characterised by intermediate and low density values, generating a fictitious flexoelectric effect.In order to avoid this artificial effect, we advocate for higher values of the flexo-penalising exponent, hence, capturing flexoelectricity only in the purely solid regions and neglecting those associated with low density regions.
With regard to the Finite Element discretisation chosen, tri or (bi-)quadratic hexahedral Finite Element interpolations of the unknown fields {φ, ϕ} and tri or bi-linear for the field G have been considered in all the examples Set design variables ρ complying with volume constraint.
while user defined criterion do Nonlinear forward problem in (62) and Appendix A. Get solution fields U = {φ, ϕ, G}.
Evaluate objective function J in terms of individual contributions J η and J τ according to (55).(see Fig. 3), for 3-dimensional or 2-dimensional cases, respectively.Although a theoretical or numerical proof of the stability or inf-sup condition [89] associated with the aforementioned discretisation is not provided, our choice of Finite Element spaces is inspired by that in Ref. [77] in the linear regime.The authors in [77] did prove the stability of their discretisation by performing a numerical test following the work of Chapelle and Bathe [90].Our discretisation is similar to that in [77], with the difference that we advocate for FM-Elasticity, whilst the formulation in [77] lies within FG-Elasticity.Therefore, the latter requires an extra degree of freedom in order to enforce (in a weak sense) the equivalence between the micro-and macro-deformation gradient tensors (i.e.G = F).
As with regard to the material parameters of the constitutive model, these are very similar to those considered in Ref. [44].The various material parameters featuring in the different energetic contributions of the free energy density, described in Sections 2.5.1-2.5.4 can be observed in Table 2. Table 2 Numerical examples.Material parameters featuring in the various energy contributions described through Sections 2.5.1-2.5.4.

Material parameters
Fig. 4. Numerical Example 1. Geometry and boundary conditions.
With regard to the length scale parameter l featuring in the micromorphic contribution of the energy density in (28) and on the scaled flexoelectric coefficient in (31), this has been taken corresponding to three times the size of the representative size of the smallest Finite Element on each of the examples.Finally, the level of discretisation of the underlying Finite Element meshes used on each example is merely motivated by the degree of resolution sought for the density field.Relatively fine Finite Element discretisations, which permit the development of fine structural elements through the topology optimisation algorithm, are in general advocated for throughout the subsequent numerical examples.

Numerical example 1: bending induced flexoelectricity
The objective of this example is to demonstrate the applicability of the proposed formulation in a bending dominated problem.The geometry and boundary conditions of the problem can be seen in Fig. 4. The Finite Element discretisation comprises {161202, 80601, 81204} degrees of freedom for the unknown fields {φ, ϕ, G}.The volume fraction for this example (volume constraint in Eq. ( 56)) was set to f = 0.4.
With regard to the electrodes depicted in 4, notice that a prescribed value of ϕ = 0 is used for the grounded electrode (bottom electrode in Fig. 4).However, the other electrode does not have a prescribed value.Instead, the condition that needs to be prescribed is that the electric potential needs to be constant on the region (line) representing the electrode.This can be done in a variety of ways.In this work, we enforce it following a penaltytype formulation, where, for all the nodes i = {1, 2, . . ., α, . . ., N } associated with the discretisation of ϕ and contained on the electrode of interest, we enforce that all the values of ϕ at those nodes, namely ϕ i , should be as close as possible to that of node α contained in the electrode.This can be done mathematically by modifying the Lagrangian in (57) as where κ ∈ R should have a sufficiently large value as to achieve the desired goal of enforcing that all the nodes have almost the same value of electric potential.The optimal designs obtained can be observed in Fig. 5, for various values of the weighting coefficient ω in Eq. (55).Clearly, when ω = 0, only the efficiency functional J η is active, which entails that the optimisation algorithm promotes the development of hinges (extremely thin structural elements), characterised by high strain gradients which obviously, yield a pronounced flexoelectric effect.This can be clearly reflected in Fig. 6.In the first row of this figure, we can observe that the highest extreme values for both Von Mises stress τ in (50) and electric field E 0 correspond with the case ω = 0.In addition, the deformed configuration corresponding with this case is considerably larger when compared with that associated with other values of ω.The designs corresponding with increasing values of ω can also be observed in Fig. 5. Clearly, for the lowest value of ω ̸ = 0, namely ω = 0, this design is very similar from a topological standpoint, to its counterpart with ω = 0.As the value ω increases, the dissimilarity in the designs with respect to the case ω = 0 become much more apparent.This increasing topological dissimilarity manifest in Fig. 6 in the descend in the extreme values for both Von Mises stress and electric field, as expected.
As observed from Fig. 6, increasing values of ω entail a reduction of localised effects.This is reflected also in the values of efficiency J η (48) and J τ (49) for the various values of ω, which are shown in Fig. 5.
It is worth mentioning that the performance of the optimised flexoelectric design might be subjected to several sources of uncertainty including aleatory uncertainties, e.g.manufacturing imperfections, unknown loading conditions and variations of the material properties, and epistemic uncertainties due to the lack of knowledge.As a result, the designed performance predicted under deterministic assumption might be degraded.This issue becomes especially important in the manufacturing of micro-and nano-structures, where spatially varying manufacturing errors may cause parts of the structure to become thinner or thicker than intended.In the same vein, spatial variations of the electromechanical properties might lead to deviation from the predicted electric field.Although the field of topology optimisation under uncertainty is out of the scope of the current research, it is worth mentioning here due to its importance from the manufacturing point of view.Several probabilistic [91] and non-probabilistic [92] formulations have been presented in the topology optimisation literature addressing aforementioned issues.However, to the best of the authors' knowledge, the problem of topology optimisation under uncertainty has not been addressed so far in the context of flexoelectricity.

Numerical example 2: compressed energy harvester
In this example, we pursue a similar analysis to that in the example in Section 4.1, particularised to the geometry and boundary conditions of the problem can be seen in Fig. 7.The Finite Element discretisation comprises {161202, 80601, 81204} degrees of freedom for the unknown fields {φ, ϕ, G}.The volume fraction for this example (volume constraint in Eq. ( 56)) was set to f = 0.4.
Unlike in the preceding one in Section 4.1, where the flexoelectric effect was induced in a bending dominated problem, this example studies the generation of flexoelectricity through compression (see Fig. 7).The optimal designs obtained can be observed in Fig. 8, for two values of the weighting coefficient ω in Eq. (55).Clearly,  when ω = 0, only the efficiency functional J η is active, which entails that the optimisation algorithm promotes the development of hinges (extremely thin structural elements), characterised by high strain gradients which obviously, yield a pronounced flexoelectric effect.In fact, an extreme rotation yielding interpenetration can be observed in the deformed configuration associated with the case ω = 0 on the left column in Fig. 9.
Once the value of ω is increased, the extreme deformation observed in 9 is considerably decreased in Fig. 9 (right column), the latter being also characterised by the presence of smaller values of the electric field and Von Mises stress distribution.

Numerical example 3: torsion-induced flexoelectricity
The geometry and boundary conditions of the problem considered in this example can be seen in Fig. 10.The Finite Element discretisation comprises {321602, 160801, 161604} degrees of freedom for the unknown fields {φ, ϕ, G}.The volume fraction for this example (volume constraint in Eq. ( 56)) was set to f = 0.5.This particular example investigates the generation of flexoelectricity through torsional induced deformations.
In addition, for a fixed value of the weighting parameter ω in Eq. ( 55), the influence in the design of the length scale parameter R featuring in the Helmholtz or PDE filter in Eq. ( 40) is investigated.Specifically, two values of R have been considered, namely R = 0.08R 1 and R = 0.027R 1 .The designs corresponding with the two values of R can be seen in Fig. 11, together with their respective values for functionals J η and J τ .Both designs show that, as it was expected, the case with smaller length scale parameter R exhibits smaller design details, specially in the four diagonal elements and their unions with the fixed regions in the centre of the design.Finally, the Von Mises stress distribution and the norm of the electric field for both designs in Fig. 11 are shown in Fig. 12.

Numerical example 4: 3D compressed flexoelectric cube
The final example explores the Flexo Micromorphic Finite Element formulation embedded into our in-house platform into a three-dimensional application.The geometry and boundary conditions are depicted and described in Fig. 13.The Finite Element discretisation comprises {985527, 85750, 385875} degrees of freedom for the unknown fields {φ, ϕ, G}, corresponding with one quarter of the domain (symmetry boundary conditions have been applied).40) and for a fixed value of weighting coefficient ω in Eq. ( 55).
For the configuration described, Figs.14a and b show the designs for a volume fraction of f = 0.45 for different values of the weighting coefficient ω in Eq. (55).The difference between both designs becomes more apparent in Fig. 15, where their respective contour plot distribution of Von Mises stress and electric field is shown.As expected, the design with ω = 0.3 yields smaller extreme values as its counterpart with ω = 0.1.In addition, the design and Von Mises stress and electric field distribution for the case of smaller volume fraction (i.e.f = 0.3) with ω = 0.3 is shown in Figs. 14 and 15 too.

Conclusions
This paper presents the first work where topology optimisation is applied for the design of dielectric elastomerbased flexoelectric energy harvester, hence, in the context of finite strains.The micromorphic continuum formulation described in Ref. [44] has been adopted for the modelling of flexoelectricity, in conjunction with a mixed-type Finite Element formulation, permitting the use of C 0 approximations for the different fields involved, namely, displacements, electric potential and the micro deformation gradient tensor.
From the topology optimisation standpoint, a novel energy interpolation scheme has been presented, whereby different penalising exponents are considered for the various contributions that the free energy function is additively A multi-objective type formulation has been proposed whereby a weighted combination of efficiency-based measures and aggregation functions of the Von Mises stress distribution within the material are considered.This was motivated by the necessity to remedy some shortcomings of the efficiency-based functional, promoting the development of hinges with highly localised large strain gradients, hence, yielding designs which are inappropriate from a structural standpoint.The effect of the consideration of the stress-based functional has been evidenced in the various numerical examples shown in Section 4. These have shown that an increasing contribution of the latter function alleviates and even eliminates the presence of hinges, which in turn, manifest in an obvious decrease of the electrical efficiency, but with the benefit of yielding designs with a better structural performance.
An important physical phenomenon inherent in dielectric-based flexoelectric devices is the possible development of electrical breakdown [88], occurring when the electric field exceeds the dielectric strength of the dielectric material.Consideration of electrical breakdown in the design of dielectric devices and flexoelectric devices, to the best of the authors' knowledge, has never been incorporated in the design problem, despite its importance.In addition, the value of the external mechanical forces in all the numerical examples on this paper has led to moderate deformations, except in the cases where a vanishing weighting factor associated with the aggregation function of the Von Mises stress distribution was considered.Although the latter parameter, up to a certain extent, helps preventing the development of mechanically induced instabilities, consideration of buckling-type constraints in the optimisation problem can also be investigated as an alternative form of preventing localised deformations leading to poor structural performance.
In the numerical examples section, the consideration of a vanishing factor ω (i.e.negligible influence of the Von Mises type objective function) yields extremely thin structural elements or hinges, which are promoted by the optimisation process as an efficient way to create strain gradients and, hence, electric fields.The deformed configurations associated with these figures show very large deformations, yielding rotations around the hinges of the design.For the cases when ω is different from zero, hinges are precluded.Nonetheless, the optimisation process promotes an alternative mechanism which favours the development of strain gradients through the creation of structural curved elements.For the values of the external mechanical loads considered in this paper, relatively moderate deformations have been obtained in these cases.Therefore, mechanically induced instabilities have not been observed.For larger values of applied forces, the global stability of the proposed design might be indeed compromised.For that case, consideration of buckling-type constraints could be incorporated with the aim of increasing the stability of the design.in (13), are discretised using the functional spaces where, for any field U ∈ {φ, ϕ, G}, n U node denotes the number of nodes per element of the discretisation associated with the field U ; and N U a : B e 0 → R, the ath shape function used for the interpolation of U .In addition, U a represents the value of the field U at the ath node of a given finite element.Consideration of the functional spaces in (A.2) enables the weak forms DΠ [ p φ ], DΠ [ p ϕ ] and DΠ [ p G ] in (13) to be written in terms of their associated elemental residual contributions, namely where N denotes the number of elements for the underlying discretisation, and where each of the residual contributions R φ a,e , R ϕ a,e and R G a,e can be expressed as2 Discretisation of the test functions and incremental fields permits Eq. ( 17) to be written in terms of their associated elemental stiffness contributions, namely Assembly over all the elements of the underlying discretisations leads to the following algebraic linear system of equations where K XY refers to the assembled counterpart of the element contribution K XY ab,e , being X and Y any member of the set {φ, ϕ, G}.In addition, in Eq. (A.8), R X represents the assembled residual vector associated with the element contribution R X a,e .Finally, ∆φ, ∆ϕ and ∆G represent the vector of nodal values associated with the unknown fields φ, ϕ and G, respectively.

Appendix B. Finite element implementation of adjoint problem
We recast Eq. (B.1), where the adjoint problems (recall that there is one adjoint problem associated with the efficiency-based functional J η (48) and another associated with the stress-based functional J τ (49)) can be written in the following compact manner with X referring to the type of functional considered, i.e.X = η, τ .Taking {δφ, δϕ, δG} ∈ V where each of the stiffness contributions is exactly the same as the stiffness matrices arising from the linearisation of the forward problem (see Eq. (A.6)).Assembly over all the elements of the underlying discretisations leads to the following algebraic linear system of equations where K XY refers to the assembled counterpart of the element contribution K XY ab,e , being X and Y any member of the set {φ, ϕ, G}.In addition, in Eq. (B.6), R X represents the assembled residual vector associated with the element contribution R X X a,e .Finally, pφ , pϕ and pG represent the vector of nodal values associated with the adjoint fields p φ , p ϕ and p G , respectively.

Fig. 1 .
Fig. 1.Deformation mapping φ (X, t).The DE in its material (B 0 ) and deformed (B) configurations.Microelement around point P ( p in spatial configuration) both in the material and spatial configurations, namely B X P 0 and B X p , respectively.

Fig. 3 .
Fig.3.Finite Element used in all numerical simulations for the two-dimensional case.

Fig. 6 .
Fig. 6.Numerical Example 1. Contour plot distribution of Von Mises stress τ (50) and the norm of the Lagrangian electric field E 0 for optimal designs in Fig. 5, corresponding with various values of the weighting coefficient ω in Eq. (55).

Fig. 8 .
Fig. 8. Numerical Example 2. Optimal designs corresponding with various values of the weighting coefficient ω in Eq. (55).Values of J η and J τ scaled with those of the case for ω = 0.

Fig. 9 .
Fig. 9. Numerical Example 2. Contour plot distribution of Von Mises stress τ (50) and the norm of the Lagrangian electric field E 0 for optimal designs in Fig. 8, corresponding with the two values of the weighting coefficient ω in Eq. (55) considered.

Fig. 10 .
Fig. 10.Numerical Example 3. Geometry and boundary conditions.The green area is fully clamped (fixed displacements).The value of the clockwise tangential force per unit length applied on the external circumference is t = 6.36 × 10 −6 N mm −1 .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 12 .
Fig. 12. Numerical Example 3. Contour plot distribution of Von Mises stress τ (50) and the norm of the Lagrangian electric field E 0 for optimal designs in Fig. 11, corresponding with the two values of the length scale parameter R in the Helmholtz or PDE filter in Eq. (40) and for a fixed value of weighting coefficient ω in Eq. (55).

Fig. 13 .
Fig. 13.Numerical Example 4. Geometry and boundary conditions for cube of side 1 mm.The two internal planes, in addition to the six lateral faces are connected to ground (ϕ = 0).Fixed displacements are imposed in the red region in the bottom lateral face.The value of the compressive force per unit area in the top face of the cube is t = 10 −6 N mm −1 .(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 15 .
Fig. 15.Numerical Example 4. Contour plot distribution of Von Mises stress τ (50) and the norm of the Lagrangian electric field E 0 for optimal designs in Fig. 14.