FEM–CM as a hybrid approach for multiscale modeling and simulation of ferroelectric boundary value problems

Constitutive modeling of ferroelectrics is a challenging task, spanning physical processes on different scales from unit cell switching and domain wall motion to polycrystalline behavior. The condensed method (CM) is a semi-analytical approach, which has been efficiently applied to various problems in this context, ranging from self-heating and damage evolution to energy harvesting. Engineering applications, however, inevitably require the solution of arbitrary boundary value problems, including the complex multiphysical constitutive behavior, in order to analyze multifunctional devices with integrated ferroelectric components. The well-established finite element method (FEM) is commonly used for this purpose, allowing sufficient flexibility in model design to successfully handle most tasks. A restricting aspect, especially if many calculations are required within, e.g., an optimization process, is the computational cost which can be considerable if two or even more scales are involved. The FEM–CM approach, where a numerical discretization scheme for the macroscale is merged with a semi-analytical methodology targeting at material-related scales, proves to be very efficient in this respect.


Introduction
The predictive investigation of electromechanical behaviors of piezoelectric, and in particular polycrystalline ferroelectric/elastic smart structures, based on computational analyses, is a demanding challenge.Processes and features on different scales have to be taken into account, from switching of crystal unit cells on the lowest scale, domain wall motion and intergranular interactions in a mesoscopic range, to field gradation due to stress concentrators on a structural macroscopic level.Mutual interactions among these issues finally require a comprehensive scale bridging modeling approach.A representative volume element (RVE), comprising a sufficient number of poly-domain grains, can be the basis of constitutive modeling, whereby sophisticated approaches nowadays reproduce experimental data in terms of commonly presented hysteresis loops very well.These are either based on a phenomenological framework [8,22,25,30,33] or rely on microphysical considerations of domain switching [7,20,21,28,57].The latter, in general, comprise a multitude of internal variables, on the other hand, however, mostly get by with relatively few parameters, whereupon physical interpretation and identification are rather straightforward.It should be noted that due to an almost inexhaustible amount of liter-ature on some of the topics taken up in this introduction, the cited references are intended only as representative examples.
A discretization scheme is required for the solution of a macroscopic boundary value problem (BVP) with regard to engineering applications.A piezoelectric stack actuator with electric and mechanical field concentrations at electrode edges [46,52] is a contemporary example, just as ferroelectric energy harvesting devices, exhibiting interfaces of smart and non-smart components [26,27], or multifunctional composites [42,56].The finite element method (FEM) is probably still the most common numerical tool for structural analyses and has been extensively exploited towards linear piezoelectric [1,13,14] and ferroelectric [9,24,39,57] material behavior.In the latter case, mesostructural evolution of polycrystals has, e.g., been considered within a micromechanical framework in [32], where domains of an RVE are allowed to switch individually without mutual interactions and driven by local macroscopic fields.In [4] each integration point (IP) of an element constitutes a grain, thus microstructural aspects are strongly correlated to the numerical discretization.Tan and Kochmann [57] recently presented a micromechanically motivated constitutive model of ferroelectric ceramics and implemented it into a 3D FEM.The sophisticated approach takes into account rate effects and thermal activation on the single grain level, while polycrystalline behavior is obtained by simple averaging over all grains, dispensing with intergranular interactions.
Scale bridging methods and homogenization approaches, respectively, providing effective quantities on the macroscale, have been the subject of research for approximately one century.Pioneering work traces back to Voigt [60] and Reuss [47], assuming constant strain and stress, respectively, in an RVE.Further classical mean field approaches such as Mori-Tanaka [40], differential scheme [44], the self-consistent method [19,31] or Hashin-Shtrikman type formulations [17,18], have long been established in solid mechanics and provide the basis of (semi-)analytical, numerical or hybrid calculations of heterogeneous structures on different scales.Contemporary contributions focusing on ferroelectrics are found, e.g., in [22,51,55], while [29] presents a general survey on modern mean field theories.
Comprehensive multiscale analysis without geometric restrictions is always based on numerical methods, inter alia incorporating concepts of the classical mean field approaches.A different idea is adopted in FE-FFT-based methods [16,41,54], including a finite element (FE) approach for the macroscale simulation and an evaluation of the behavior at the microscale based on fast Fourier transforms.Another powerful contemporary tool for numerical homogenization and two-scale simulation is the FE 2 [11,53], where nested FE calculations of macro-and meso-or microscales are carried out in separate models with a con-tinuous exchange of information.The local macroscopic consistent tangent is thus derived from the simulations on the lower scale.More recently, this approach has been elaborated for multiphysical problems [49,59].Despite of efforts to reduce the computational cost, e.g., in terms of a monolithic solution scheme in [34], the FE 2 method is still among the computationally most expensive two-scale approaches.
Semi-analytical or hybrid methods aim to cope with this drawback, still enabling the solution of arbitrary macroscopic BVPs with microstructural interactions on one or several scales.In [62] a Hashin-Shtrikman type formulation is taken as basis for the microscale.Other methodologies are, e.g., the uniform or non-uniform transformation field analysis [10,12] and more recently the self-consistent clustering analysis [38,63].
In the work at hand, the so-called condensed method (CM) [35,37,48] is exploited within a FE framework, combining the efficiency of the CM in calculating micro-and mesostructural evolutions from unit cell to polycrystalline RVE level with the flexibility of the FEM in solving macroscopic BVPs.The CM constitutes a semi-analytical approach for the simulation of nonlinear multiphysical constitutive behavior, accounting for interactions of the scales involved, in particular of grains in an RVE or domains in a grain.Resulting residual fields, in return, contribute to evolution equations and energy dissipation.Since 2015 the CM has been successfully applied to various problems, including damage and life span of ferroelectrics exposed to electromechanical cyclic loading [36], ferromagnetic modeling and multiferroic composites [48], tetragonal-rhombohedral ferroelectrics and loading-induced phase transition [48,58], ferroelectric heating [61] and optimization of ferroelectric energy harvesting [6].
Spanning four scales of a ferroelectric device from unit cell to grain to polycrystalline RVE to macroscopic structure, the presented approach, denoted as FEM-CM, is highly efficient with regard to computational cost, and has previously been validated on the constitutive level.In contrast to the FE 2 or FE-FFT methods, on the other hand, it does not allow for a spatially resolved micro-/mesostructure simulation, and is thus not suitable if geometrical shapes, deterministic arrangements etc. are of predominant significance.For statistically arranged micro-and mesostructures, e.g., associated with grains or ferroic domains, despite of their deterministic orientations, the FEM-CM turns out to be a very powerful approach.
In the following sections, the theoretical basis of the CM is outlined in brief.The interfaces of CM and FEM are further depicted, just as essential specifics of the FEM-CM implementation.Sampling points (SP) are introduced on an independent grid, in order to decouple mesostructural information, in particular grain size distribution, from the FE discretization.Being elaborated for the general 3D case so far, examples are restricted to the 1D case of a rod element, which is appropriate to demonstrate numerical issues, and finally to prove the appropriateness of the approach.

Some phenomenological basics of ferroelectric continua
Nonpolar solid continua under quasistatic conditions are considered in the following.Employing analytical notation as a mathematical basis, whereat summation is indicated by repeated indices, the balance equations read: Equation ( 1) describes the balance of momentum, containing Cauchy stresses σ i j and volume forces b i , whereas Eq. ( 2) implies electrostatic equilibrium with electric displacement D i , disregarding volume charges.Here, as well as in the following, a comma in a subscript denotes partial spatial derivation.Moreover, infinitesimal mechanical deformations are assumed, leading to a linear relation between the strain ε i j and the displacement gradient: On the other hand, the electric field E i is related to the gradient of the electric potential φ: Due to dissipation as a result of domain wall motions, thermodynamic consistency is a major issue in ferroelectric constitutive laws.In e.g.[37], the electric enthalpy density h as a thermodynamic potential of ferroelectric continua reads: with the electric field and the total strain as independent variables.Internal variables, accounting for domain wall motion, are involved in the elastic tensor C i jkl , the dielectric tensor κ i j and the piezoelectric coupling tensor e ikl as well as the irreversible strain ε irr i j and polarization P irr i .Differentiating Eq. ( 5) with respect to the independent variables, the associated variables stress and electric displacement are obtained as representing the ferroelectric constitutive law.Equations ( 6) and ( 7) are valid within sufficiently small changes of state, at which material tensors act as linear tangent moduli.Inelastic nonlinear ferroelectric behavior is achieved adapting material tangents and irreversible quantities iteratively, see Sect.3.1.
Considering thermodynamic consistency, a mathematical description of ferroelectricity has to satisfy the balance of entropy.Starting at the Clausius inequality according to where q i , ρ, r and denote specific heat flux, mass density, volume heat source and temperature, respectively, the generalized Clausius-Duhem inequality is obtained by inserting the local energy balance into Eq.( 8).
Here, s denotes the specific entropy, u the specific internal energy and the strain and electric displacement are decomposed into reversible (ε rev i j , D rev i ) and irreversible (ε irr i j , P irr i ) parts.While the statement −q i ,i / ≥ 0 is always fulfilled due the oppositional directions of temperature gradients and heat fluxes q i , the remaining terms hold the inequality in case of reversible processes, i.e., εirr i j = 0 and Ṗirr i = 0 [45].Consequently, an irreversible change of state due to domain wall motion has to satisfy for the sake of thermodynamic consistency.

Multiscale modeling of a ferroelectric continuum
In order to obtain realistic macroscopic hystereses as a result of processes in grains, domains and crystal lattices and to be able to investigate arbitrary structures, multiscale modeling is indispensable.For this purpose, four different scales are defined and illustrated in Fig. 1.
Here, as in the following, a short notation is introduced representing all material coefficients, i.e.

Transition of scales micro and meso-1
A transverse isotropy is assumed, at which the local polarization vector is perpendicular to the isotropic plane.Accordingly, material coefficients of a domain in a global coordinate system depend on the direction of spontaneous polarization of the associated unit cells.Due to the heterogeneity of the domain structure, properties of a grain, on the other hand, have to be described as effective quantities resulting from volume averaging, according to an approach by Huber et al. [20] introducing volume fractions as internal variables.
A polycrystalline MVE with a representative excerpt of a grain m, indicating the 3D tetragonal domain structure, is illustrated in Fig. 2. In the model each of the N = 6 species n is defined by its direction of spontaneous polarization along with their corresponding volume fractions ν (n) , satisfying the conditions It should be noted that the microphysical approach works for arbitrary types of crystals, e.g., rhombohedral or mixed types [48].The effective properties of a grain are obtained from the material tensors of its allocated domain species by weighted summation: Fig. 2 3D domain model of a grain consisting of tetragonal unit cells with vectors of spontaneous and average polarization and internal variables ν (n) representing volume fractions The average polarization P (m) i is derived accordingly.The internal variables ν (n) are controlled by a switching criterion [21], resulting in an evolution law [37]: Accordingly, the dissipative work w diss of unit cells in n necessarily has to reach an energy barrier w crit in order to initiate switching.While this requirement is stipulated by the first Heaviside function H(...), the second one satisfies the principle of potential energy minimization in terms of ensuring the maximum energy of all possible options to dissipate.Here, the superscript k denotes the adopting species for the providing domains n, consequently satisfying Eq. ( 12).Furthermore, dν 0 is a model parameter, expressing the magnitude of incremental volume change within an iteration step.It has to be chosen sufficiently small for the sake of numerical stability and large enough with regard to computational cost.The energy barriers for 90 • and 180 • switching of tetragonal unit cells are given as [20] w involving the magnitudes of spontaneous polarization P 0 and the coercive field E c .The dissipative energy, on the other hand, reads consisting of strain as well as dielectric energies of switching unit cells and always being positive, see Eq. ( 10).The evolution of material tensors according to Eq. ( 13) is disregarded in Eq. ( 17), making a significantly smaller contribution to dissipation.Furthermore, the evolution law, in contrast to e.g.[57], does not account for rate dependence of switching or aspects of thermal activation.Spontaneous strain ε sp i j and the change of spontaneous polarization P sp i are determined from lattice parameters, see Table 2.The mechanical stress and the electric field are considered as driving forces in Eq. (17), which are assumed constant during switching.It will be illuminated in the following section from which scale both quantities have to be adopted.Moreover, domain wall motion causes irreversible strain and polarization of a grain m, depending on the spontaneous strain ε sp(n) i j and polarization P sp(n) i , respectively, associated with each domain species n [37]: Unlike material tensors, the irreversible quantities are recursively described by the change of the internal variables, depending on the loading history.The Heaviside functions guarantee ε irr(m) i j and P irr(m) i only to be modified by donating domain species.Alternatively, ε sp(n) i j(n→k) could have been introduced, relating the n tensors of spontaneous strain to the cubic phase constituting a reference or interim configuration in terms of switching [43].In this case, Eq. ( 18) basically takes the form of Eq. ( 13).

Transition of scales meso-1 and meso-2
The meso-1-meso-2-transition is based on the so-called condensed method [35,48], which is a semi-analytical approach bridging scales in multiphysical systems accounting for interactions of heterogeneities.In the following, may be an arbitrary quantity, including independent and associated fields, material tensors according to Eq. ( 11) as well as ε irr i j and P irr i .Assuming local fluctuations due to the grain structure on the meso-2-scale, effective macroscopic quantities are considered on the MVE level, whereupon angled brackets denote volume averaging: As outlined in [35], simplified polycrystals are considered which are characterized by homogeneous grains of uniform size, i.e., = f(x i ) in each grain.Thus, Eq. ( 19) evaluates to whereat M denotes the total number of grains of a MVE and (m) is given, e.g., by Eq. ( 13) in case of material tensors.
Accordingly, in Fig. 1 the model of the polycrystal (meso-2) is illustrated by concentric crosses, representing the prevailing polarizations of each grain, distinguished by different colors.In [37] a stochastic grain size distribution was implemented in connection with Eq. ( 19).Compared to a uniform grain size it was found that its impact is negligible, unless a relative standard deviation of 15-20% is exceeded.
In order to derive the macroscopic constitutive law from the meso-2 scale, Eqs. ( 6) and (7) are likewise applied to a grain m: While material tensors and irreversible quantities are continually updated based on the evolution law, the independent and associated variables have to be determined subsequently.Therefore, a generalized Voigt approximation is considered assuming homogeneous total strain and electric field inside each MVE, i.e. ε i j , E i = f (m), which is marked by an overbar.Consequently, this leads to With these relations the macroscopic constitutive law results from inserting Eqs. ( 21) and ( 22) into Eq.( 20), yielding: with Selecting, e.g., σ i j = σ ext i j and E i = E ext i as loading quantities, the strain of a MVE results from the constitutive law Eq.( 25) as and the electric displacement is obtained from Eqs. ( 26) and (29) according to Interactions between grains, inducing residual stresses and electric displacements, are involved in Eqs. ( 29) and ( 30) by imposing constraints upon strains and electric fields of grains in order to satisfy Eqs. ( 23) and (24).Butterfly and dielectric hysteresis loops thus exhibit smooth shapes as known from experiments, see Sect. 5. Individual stresses of each grain, with external and residual contributions, are derived as follows: Consequently, the dissipative work of Eq. ( 17) has to be adapted in the evolution law according to since individual stresses of grains act as driving forces of domain wall motions within a generalized Voigt approximation of the CM.

Transition of scales meso-2 and macro
Macroscopic problems can be handled on the basis of Eqs. ( 29) and (30), as long as a MVE is representative for the whole structure.As soon as, e.g., notches or electrode edges are involved, a discretization scheme is required.In order to solve arbitrary BVPs, the meso-2-macro transition is realized by incorporating the multiscale constitutive law of Eqs. ( 25)- (28) in conjuction with the evolution law Eq.( 14) and the dissipative work according to Eq. ( 32) into a FE approach.
Starting at the principle of virtual displacements with δW a denoting the virtual work of external loads, the internal potential i is related to the electric enthalpy density introduced in Eq. ( 5) where V Cont represents the volume of the continuum on the macroscale and R denotes the number of MVEs in a structure.
Again, the angled bracket indicates the volume average, here of the electric enthalpy density in a MVE to be considered at this point, constituting one essential scale-bridging interface.Consistently, Eqs. ( 25) and ( 26) are obtained by differentiating h with respect to the macroscopic independent variables ε i j and E i .The weak formulation of the coupled nonlinear electromechanical BVP follows from Eq. ( 33), applying the first variation to the internal potential of Eqs. ( 34) and ( 35), at which the latter three integrals collectively describe the virtual work of external volume and boundary loads.Equation (36) provides a basis for the multiscale FE analysis of a polycrystalline, polydomain ferroelectric structure.Whereas ε i j and E i have been introduced at transition from meso-1 to meso-2 levels, their interpretation within a FE context will be enlightened in the following section.

Finite element formulation
In this section, fundamentals of the FE implementation of the multiscale computational approach are briefly outlined first.Subsequently, the problem of decoupling mesostructural properties from geometrical discretization is addressed.
Finally, the focus is on a quadratic rod element, which will be the basis of first investigations to demonstrate the feasibility of the approach.

General approach
In the following, an algebraic system of linear equations (SLE) of a FE is obtained from the weak formulation according to Eq. ( 36), assuming the structure to consist of a finite number of elements with individual volumes V El , generally not coinciding with the number of MVEs R. A generalized displacement method is considered, at which mechanical displacement and the electric potential act as primary nodal variables.Furthermore, isoparametric FEs ensure numerical efficiency [5], interpolating nodal quantities in an element according to where the index k marks the K nodes of a FE and the interpolation functions depend on nondimensional local coordinates (ξ, η, ζ ).Curly brackets indicate column matrices with dimensions 3K for { x}, { û} and K for { φ}, respectively, in a 3D problem and the circumflex implies nodal variables.Furthermore, the approximated independent variables at each MVE are related to the primary nodal variables as whereat [D u ] and D φ denote differential operators regarding the relations in Eqs. ( 3) and ( 4) and the Voigt index p ranges from 1 to 6. Bridging the scales of meso-2 and macro, Eq. ( 38) constitutes a second interface of CM and FEM, in that {ε} and {E} appear in Eqs. ( 25), ( 26) and (32).In this context, σ (m) i j in Eqs. ( 31) and ( 32) is substituted by Moreover, it is obvious that a generalized Voigt approximation on the meso-2 scale requires the generalized displacement method, as the interpolation functions in Eq. ( 38) deliver total strain and electric field for each MVE.
Consequently, an alternative approximation, e.g., a mixed Voigt-Reuss based on strain and electric displacement, comes along with another interpolation approach in terms of displacement and charge as primary nodal variables.
Considering the fundamental lemma of calculus of variations, the algebraic SLE of a FE results from Eq. ( 36): at which the generalized stiffness matrix includes elastic, piezoelectric and dielectric contributions as defined by corresponding effective material properties: On the other hand, the nodal loads read containing surface loads and volume forces as well as contributions due to domain wall motions in { a irr,σ } and { a irr,D }.The former will vanish, if ∂ V El is not part of the structure's boundary.
The flowchart in Fig. 3 illustrates how to solve a BVP based on the outlined multiscale FEM-CM formulation.In addition, a procedure restricting to the investigation of polycrystalline constitutive behavior based on the CM is depicted in Fig. 3b for comparison.Introducing individual material tensors as a result of different grain orientations, associated random angles are initially prescribed in each MVE.Afterwards, material coordinates as well as ε sp(n) i j and P sp(n) i of each domain species are transformed into a global coordinate system.Subsequently, the simulation begins imposing external loads in an outer loop.An inner loop, on the other hand, focuses on the evolution of the internal variables, thus

Decoupling of grain structure and FE discretization
In most of the FE implementations of micromechanical ferroelectric constitutive models, FEs [2,3] or IPs of each FE [4,15,50] represent grains or multiple-grain MVEs, see Fig. 4a.In a physically inconsistent manner, this approach thus implies a correlation between the FE mesh and the mesoscopic structure.Therefore, a more suitable approach is realized at this point, introducing a so-called MVE grid, see Fig. 4b.It implies an arrangement of SPs, being independent from the FE mesh.Each SP is assigned a MVE, thus the domain evolution is calculated according to Fig. 3 at these points.At any other positions the macroscopic material tangents and irreversibilities are obtained by interpolation.For the calculation of the integrals of Eqs. ( 41) and ( 42) values are required at the IP of the FE mesh.Besides being a fundamental issue, this concept is particularly reasonable if stress gradients, e.g., at notches or electrode edges, require local mesh refinement.Furthermore, a variation of IPs for the sake of accurancy of calculations must not be linked with mesostructural length scales.Whereas the SPs in Fig. 4 are simply located on a square grid for the sake of illustration, their arrangement can generally be chosen arbitrarily.Since the distribution of SPs represents the local grain density as mesostructural feature, a refinement of the MVE grid, in contrast to h-convergence of the FE mesh, does not target convergence of results.Details of the approach are presented in the following section.

The ferroelectric rod element
The FE formalism and the MVE grid are initially implemented to one-dimensional rod elements.For the sake of demonstration, a three-nodes FE with a length l and the cross section A cross (ξ ) is presented, exposed to a uniform line load n mech (ξ ), see Fig. 5.It is characterized by the nodes A, B, C and its nondimensional local coordinate ξ .Furthermore, a global coordinate system x i is given, also providing a basis for the coordinate transformation of material tensors.Firstly, the primary field variables u 3 and φ are isoparametrically approximated, based on Lagrange polynomials: A one-dimensional FE requires a suitable formulation of the macroscopic constitutive law, only containing axial components σ 3 , D 3 = f (ε 3 , E 3 ).Therefore, a uniaxial state of stress and electric field (USE) is assumed in the following, leading to a reduced representation of Eqs. ( 25) and ( 26 In contrast to the electric field and macroscopic stress, each component of ε p and D i generally is different from zero.The coefficients C 33,red , e 33,red and κ 33,red depend on the ten elastic, piezoelectric and dielectric constants and the orientation of each grain in the MVE.These relations are not given explicitely, constituting extensive equations due to statistical domain orientations and related fully populated matrices in the global coordinate system.The reduced algebraic SLE of the rod element is obtained including the equivalent axial nodal loads: The 3x3 stiffness sub-matrices follow from Eq. ( 41), yet incorporating line integrals and the reduced effective material coefficients of Eq. ( 44), while { F S 3 } and { QS } represent surface loads at boundary nodes.
Table 1 Integration points ξ g and weights w g of Gaussian quadrature  ν 0 = 0.001

Verification of the approach
In this section various aspects of physical and numerical plausibility of the FEM-CM approach are demonstrated in simple examples.A one-dimensional problem turns out to be appropriate for this purpose.Therefore, each BVP in the following incorporates only one unilaterally fixed FE with the length l = 1 m and a constant cross section A cross = 1 m 2 , disregarding any volume forces.Either the rod element of Sect.4.3 is considered or an analogous two-nodes element with linear interpolation functions and boundary nodes A and B. Numerical integration is realized by the Gaussian quadrature.Integration points ξ g and weights w g up to G = 5 IPs are given in Table 1.Furthermore, material data of barium titanate [23] are adopted, see Table 2, and each MVE exhibits M = 75 statistically orientated initially unpoled grains, i.e., . The rod is assumed to be embedded in air with a dielectric constant being three orders of magnitude smaller, thus the electric energy stored in the environment is negligible.

The FEM-CM interfaces
The correct implementation of the FEM-CM interfaces has to be confirmed at first.Generally, for a rod with homogeneous state of stress and thus without volume forces and a constant cross section, one two-nodes rod element with linear interpolation functions is appropriate to provide the exact analytical solution.This issue analogically holds for a USE in a piezoelectric rod element.Accordingly, the interfaces are correctly implemented, if solutions of macroscopic field quantities of the FE approach equal the semi-analytical ones of a pure CM implementation, see Fig. 3. Therefore, a macroscopically stress-free rod is considered, exposed to a cyclic axial electric field with maximum Homogeneous material behaviour is assumed by considering equivalent MVEs at each position.Numerically, the microstructured BVP is solved by means of a two-nodes rod element with the clamped end A and the free end B. The electrical load is provided by appropriate nodal potentials according to Eq. (38).Due to the material homogeneity, the.(45) are constant, whereupon full integration is achieved by a single IP incorporating the corresponding MVE.In Fig. 6 butterfly hystereses are compared based on the FEM-CM calculations of the single IP and a single MVE directly derived from the CM.Apparently, the curves are congruent, implying equivalent residual stresses as well, which gives evidence of a correct implementation of the FEM-CM interfaces.

Satisfaction of electromechanical equilibrium
Material heterogeneity is implemented now by different grain orientations of each IP, in Fig. 7 indicated by different colors.In order to investigate the equilibrium of charges and forces, a polarization process of the three-nodes rod element is considered, followed by a mechanical tensile loading.The former is realized by an increasing nodal potential at the free end up to φ C = 5 • 10 5 V, while φ A remains zero.The Fig. 6 Axial strain of an electrically loaded two-nodes rod element calculated with the FEM-CM approach compared to the result of a single MVE based on the CM mechanical boundary conditions are again fixed clamping at node A and traction-free boundary at node C for the polarization process.Concerning the second scenario, the electric boundary condition at the free end is replaced by the nodal charge Q C = 0, while an increasing tensile force is imposed up to F C 3 = 10 3 MN.Since volume loads are neglected, each nodal force at the polarization process and each nodal charge at tensile loading have to be zero, regarding the balance of forces and charges, respectively.In fact, the macroscopic stress and electric displacement have to vanish globally in their corresponding loading scenarios.Figure 8a shows σ 3 versus the nodal potential at each IP ξ g , while the macroscopic electric displacements are plotted versus the nodal force F C 3 in Fig. 8b.It is obvious that the balance of momentum is locally infringed above the coercive field, inter alia due to residual grain stresses during the polarization process.The same holds for the electrostatic equilibrium during mechanical loading from the beginning.
In contrast to the microstructured homogeneous element, the individual grain orientations at each MVE lead to a statis-tical variance of material parameters and irreversible strain and polarization along the rod axis.The strain ε 3 and electric field E 3 , on the other hand, are spatially linear due to the constraints of Eq. ( 43).This local mismatch causes the violation of the equilibrium conditions.Furthermore, the graphs at the outer IPs g = 1, 3 are basically similar, noticeably differing from the one at g = 2.This issue might be attributed to the fact that the former exhibit the same central distance as well as identical Gaussian weights.
At this point, it has to be verified, if balance equations are globally satisfied.In this context the arithmetic means of IP values are considered as well as their weighted means.The latter is obtained by integration in terms of Gaussian quadrature according to taking the fluctuation of axial macroscopic quantities into account.In Fig. 9 arithmetic and weighted mean values of the results of Fig. 8 are compared to each other.While the former do not vanish, the weighted mean values satisfy the balance equations as expected for a weak formulation of the BVP.

Influence of numerical integration on the evolution of the mesoscopic structure
Due to the fluctuation of material quantities in a microstructured heterogeneous body, it is obvious that the number of IPs required for full integration of the generalized stiffness matrix and irreversible quantities comes along with the polynomial degree of (ξ ) and a irr (ξ ).Since the latter cannot be prescribed a priori, the influence of numerical integration on results is investigated by comparing analyses of two twonodes FEs, differing in the number of IPs, i.e., G = 2 and G = 3.In order to create similar material properties, the Fig. 7 Rod element clamped on one side with inhomogeneous material behavior, represented by different MVEs at the integration points ξ g 123 Fig. 8 Macroscopic stresses of the polarization process (a) and macroscopic electric displacements of the mechanical loading process (b) at the Gaussian integration points ξ g Fig. 9 Comparison of arithmetic and weighted means of integration point values shown in Fig. 8 outer IPs exhibit identical grain orientations in both FEs, see the colors in Fig. 10.
Both curves coincide qualitatively with a slight deviation in the maximum strains.Changing the number of IPs and thus of MVEs effectuates different distributions of polycrystals along the element, consequently having an impact on 33,red (ξ ) and the stiffnes matrices according to Eq. ( 41), finally leading to the deviation observed in Fig. 11.This discrepancy might be non-negligible considering BVPs with large gradients, e.g., at notches or electrode tips [15], motivating the introduction of the MVE grid.

Verification of the MVE grid approach
In order to verify the concept of the MVE grid of Sect.4.2, decoupling the mesoscopic structure from numerical integration and circumventing a correlation with the FE discretization, two FEs with different numbers of IPs are compared to each other.
The three-nodes rod element with an equidistantly distributed MVE grid of R = 5 SPs is considered, starting from a boundary node with ξ (r ) = 1/2.Heterogeneity is realized by different polycrystalline MVEs at each SP r , in Fig. 12 indicated by different colors and orientations of polarization arrows, coming along with fourth-degree Lagrange polynomials L (r ) (ξ ) in Eq. (46).Taking the constant cross section into account, coefficients of the generalized stiffness matrix in Eq. ( 41), exhibiting the largest power in ξ , generally have the polynomial degree of six, whereupon full integration requires at least G = 4, see Fig. 12.While mechanical loads are disregarded, a cyclic potential difference with φ A = 0 V Fig. 10 Two inhomogeneous two-nodes rod elements with different number of integration points and MVEs, respectively Fig. 11 Comparison of butterfly hystereses of two heterogeneous twonodes rod elements, differing in their number of integration points G, see Fig. 10 and |φ C,max | = 10 6 V is imposed.In Fig. 13 the axial strains of both elements with G = 4 and G = 5 are compared to each other, fading out the initial polarization process.In the plots (a) and (b) different positions ξ a/b = ±2/5 have been chosen for the evaluation, not coinciding with neither SPs nor IPs.
In contrast to Fig. 11, the butterly hystereses are congruent at both positions, proving the suitability of the MVE grid approach to decouple the mesoscopic structure from IPs.Furthermore, the axial strains are locally identical, revealing negligible gradients related to an equidistant distribution of SPs.Compared to, e.g., Fig. 11, slight fluctuations are observed in the loading branches of the butterfly hystereses based on the MVE grid approach, whereas remanent and maximum strains remain almost unaffected.
In Fig. 14 the axial stress σ 3 (φ C = 10 6 V) along the rod element with G = 4 is depicted, while Fig. 15 shows the spatial distribution of the electric displacement D 3 (φ C = 10 6 V).Values at SPs and IPs are highlighted along the graphs and the interpolation of SP values is based on Eq. (46).The colored sketches are intended to illustrate the fields along the rod axis.As expected in Sect.5.2, the heterogeneous mesostructure exhibits an obvious fluctuation of the associated variables.This effect is specifically caused by the individual evolutions of e 33,red (ξ ) and ε irr i j (ξ ).Figures 14 and 15 might give the rough impression of an axial symmetry of the graphs with respect to ξ = 0 which, however, may evolve fundamentally different based on different grain orientations in the MVEs.The weighted mean satisfying balance equations is again confirmed in Fig. 14.

Conclusion
A hybrid multiscale approach, denoted as FEM-CM, which is highly efficient regarding the computational investigation of ferroelectric smart structures, is presented.The FEM as a numerical discretization tool is employed in order to solve arbitrary boundary value problems on the macroscale.The constitutive behavior, however, is described on the mesoscale by the CM, constituting a semianalytical homogenization technique, providing effective quantities of a polycrystalline volume element based on its multiscale evolution processes.The switching of unit cells on the microscale is taken into account as well as domain wall motion and grain interaction on mesoscopic scales.
While the methodology is elaborated for general 3D cases, 1D rods with only one FE have been investigated for the sake of demonstration and verification.After proving the appropriateness of the FEM-CM interfaces, as well as mechanical and electrostatical consistencies of the approach, an unphysical coupling between the grain structure and numerical aspects has been revealed if mesoscopic volume elements are connected to integration points of FEs.Therefore, a more sophisticated approach is given by a so-called MVE grid, at which mesoscopic volume elements are located at fixed sampling points in the whole body, being independent from integration points and thus from numerical discretization.
Disregarding spatially resolved simulations on the microand mesoscales, the FEM-CM is extremely efficient, exhibiting advantages in pre-processing and computational cost compared to established multiscale approaches like the FE 2 or FE-FFT methods.Extensions towards ferromagnetic or Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Fig. 1
Fig.1Classifications of the multiscale modeling approach and definition of length scales (only plane states of polarization depicted for the sake of illustration)

Fig. 3 Fig. 4
Fig. 3 Flowchart of the multiscale FE approach based on the CM

Fig. 5
Fig. 5 Ferroelectric three-nodes rod element with volume force and boundary loads