Numerical algorithms based on the Finite Element Method will be specialized for Analysis, Design, and Optimization of Sensors and Actuators (S-A) and their Application to Smart Structures. The S-A based on tangible assets can couple several fields, such as mechanical, electrical, magnetic, and thermal. They are used in many applications, particularly in smart structures, damage monitoring, or aerodynamics. Despite the considerable experience in these studies, the steps addressed are first to develop a thermodynamically consistent formulation for macro-scale to introduce plasticity models; second, to provide the tools to take into account the heterogeneities of multi-scale models for smart materials. The main objective is the development of a research computer code to simulate and study the performance, not only of the S-A themselves but also of the smart structures in which these S-A will be mounted.

Résumé

Les algorithmes numériques basés sur la méthode des éléments finis seront spécialisés dans l'analyse, la conception et l'optimisation de capteurs et d'actionneurs (S-A), ainsi que dans leur application aux structures intelligentes. Les S-A basés sur des actifs tangibles peuvent coupler plusieurs domaines, tels que les domaines mécanique, électrique, magnétique et thermique. Ils sont utilisés dans de nombreuses applications, notamment dans les structures intelligentes, la surveillance des dommages ou l'aérodynamique. Malgré l'expérience considérable de ces études, les étapes abordées consistent d'abord à développer une formulation thermodynamiquement cohérente à l'échelle macro pour introduire des modèles de plasticité; deuxièmement, fournir les outils permettant de prendre en compte les hétérogénéités des modèles multi-échelles pour les matériaux intelligents. L'objectif principal est la mise au point d'un code informatique de recherche permettant de simuler et d'étudier les performances, non seulement des S-A eux-mêmes, mais également des structures intelligentes dans lesquelles ces S-A seront montés.

List of Figures

-Entropy density (J/K m 3 ) r -Heat source (W/m 3 ) s n , s n -Polarization and fracture multipliers (-) q -Heat flux (W/m 2 ) j -Jacobian matrix (m) q, q, q -Stress-like hardening variable (Pa) c p -Specific heat (J/kg K) e

Introduction

In this first chapter, a brief description of smart materials and its main applications of interest is exposed. The wide range of applications determines the motivation and objectives of this thesis. A brief review of the existing literature on the numerical aspects of smart materials is followed right after. The outline of the thesis is presented at the end of the chapter.

1 SECTION 1.1. Motivation

Motivation 1.Smart materials

The concept of smart materials (SMs) comprehends a wide range of materials that can react to external stimuli and modify a particular property of the material or produce a response in another field. In this sense, the term coupled material is typically used to define a SM that has two or more physical fields, such as mechanics or electrics, inherently interacting with each other.

The main difference with traditional materials is that the response or the property change is spontaneous, useful to the purpose of the system, reliable, and predictable (usually proportional) to the intensity of the stimulus, such as external conditions like temperature, light or pressure. Therefore, this sort of materials is suitable for its use as sensors and actuators on a more complex system or structure, not requiring any control system, making a smart structure.

The field of SMs is interdisciplinary since it affects many disciplines of engineering (sensors, actuators, control, information processing, material science) and different branches of physics (electromechanics, thermodynamics, mechanics, fluids, chemical reactions). This interdisciplinarity makes the state of research and development of SMs mostly disparate, depending on the maturity of the applications or the branches involved for a particular material. Namely, piezoelectric, electrostrictive ceramics, piezoelectric polymers, and fiber-optic sensors have reached a high grade of maturity, and they are well-established commercial technologies. On the other hand, micromachined electromechanical systems (MEMS), magnetostrictive materials, shape memory alloys (SMA) and polymers, and conductive polymers are in the early stages of commercialization. The next generation of SMs will introduce chromogenic materials and systems, electro-and magneto-rheological fluids, and biometric polymers and gels [START_REF] Schwartz | Encyclopedia of smart materials[END_REF].

These materials are unparalleled in what they can offer. On top of that, the general advantages of the integration of SMs into more complex devices or structures are their reduced dimensions, lightness, and, usually, their autonomy from a power source to function. Besides, they can be controlled and combined with computer systems. Therefore, SMs are linked to the main objectives of our generation: mobility, energy, environment, and health [START_REF] Bengisu | Materials that Move: Smart Materials, Intelligent Design[END_REF].

Research objectives

This thesis is focused on the numerical aspects of the implementation of smart materials that couple thermal, electric, magnetic, and mechanic linear interactions into the finite element method at different scales. In particular, the research aims to provide:

• A macro-scale, hyperelastic, thermodynamically-consistent finite element formulation for both primary electric constitutive behaviors, conductors and insulators.

• The introduction of a macro-scale non-linear model for the mechanical field on the previous models where the interactions with other fields affect the performance of the solid.

• Explore the possibilities to improve the discretization of electromagnetic fields through differential forms and exterior calculus theory.

• Develop a microscale lattice model for piezoelectric polarizable ceramics with the inclusion of domain switching and crack discontinuities.

1.3 Literature review

Coupling formulations

In a coupled formulation, the presence of electromagnetic fields introduces body forces that must be accounted for in the conservation principles for linear and angular momenta. From a thermodynamics point of view, the addition of electromagnetics produces, in general, the appearance of irreversible phenomena in the form of dissipative terms. Therefore, the energy conservation equation and thermal field are needed to take into consideration the loss of efficiency due to these dissipative terms [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF].

A significant difference must be established between the two main electric constitutive models in terms of conductors versus insulators. The electrons for the conductors are free to move, whereas for the insulators the electrons remain bound to their nucleus. Insulators react to electric fields by orienting dipoles, deformed atoms or molecules with two opposite charges held in place by molecular forces. Thus, these charges are not free to move as it happens with conductors [START_REF] Balanis | Advanced engineering electromagnetics[END_REF], [START_REF] Jackson | Classical electrodynamics[END_REF]. The different behaviors of conductors and insulators influence the possible couplings in the constitutive equations [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF].

In the literature there are numerous scientific articles on the first order coupling interactions and their formulation for insulators of just two or three of the mentioned fields, such as piezoelectric [START_REF] Allik | Finite element method for piezoelectric vibration[END_REF], [START_REF] Damjanovic | Electrostrictive and piezoelectric materials for actuator applications[END_REF], [START_REF] Lezgy-Nazargah | An efficient finite element model for static and dynamic analyses of functionally graded piezoelectric beams[END_REF], [START_REF] Safari | Piezoelectric and acoustic materials for transducer applications[END_REF] and a special mention to [START_REF] Duczek | Anisotropic hierarchic finite elements for the simulation of piezoelectric smart structures[END_REF] for the development of a piezoelectric element based on p-version finite element formulation first introduced in [START_REF] Babuska | The p-version of the finite element method[END_REF]; thermo-electro-elastic as in [START_REF] Ryu | Magnetoelectric properties in piezoelectric and magnetostrictive laminate composites[END_REF], [START_REF] Wang | A finitely long circular cylindrical shell of piezoelectric/piezomagnetic composite under pressuring and temperature change[END_REF], [START_REF] Ferrari | Thermodynamic formulation of the constitutive equations for solids and fluids[END_REF]; or electro-magneto-elastic as in [START_REF] Ramirez | Free vibration response of two-dimensional magnetoelectro-elastic laminated plates[END_REF], [START_REF] Görnandt | Finite element analysis of thermopiezoelectric smart structures[END_REF], [START_REF] Fung | Dynamic analysis of a piezothermoelastic resonator with various shapes[END_REF], [START_REF] Rao | Analysis of distributed thermopiezoelectric sensors and actuators inadvanced intelligent structures[END_REF], [START_REF] Jiang | A new finite element model for piezothermoelastic composite beam[END_REF], [START_REF] Hou | The transient responses of a special non-homogeneous magneto-electro-elastic hollow cylinder for axisymmetric plane strain problem[END_REF] and [START_REF] Pan | Exact solution for simply supported and multilayered magneto-electro-elastic plates[END_REF]; other papers couple all fields but without a very consistent thermodynamic approach as in [START_REF] Li | Magnetoelectroelastic multi-inclusion and inhomogeneity problems and their applications in composite materials[END_REF] or [START_REF] Aboudi | Micromechanical analysis of fully coupled electro-magneto-thermo-elastic multiphase composites[END_REF]. The recent ones that come closest to the objectives of this thesis, not only because of their full-interaction approach but also for their thermodynamic framework are [START_REF] Chen | General solution for transversely isotropic magneto-electrothermo-elasticity and the potential theory method[END_REF], (Pérez-Aparicio et al., 2016).

Regarding the conductors, it is more difficult to find a comprehensive approach to the first order full-interaction between all four fields. Relevant articles on the basics of thermoelectricity (that covers the Peltier, Seebeck, and Thomson effects) include [START_REF] Pérez-Aparicio | Finite element analysis of nonlinear fully coupled thermoelectric materials[END_REF], [START_REF] Rowe | CRC handbook of thermoelectrics[END_REF][START_REF] Disalvo | Thermoelectric cooling and power generation[END_REF], [START_REF] Antonova | Finite elements for thermoelectric device analysis in ansys[END_REF], [START_REF] Zhu | Finite element analysis of miniature thermoelectric coolers with high cooling performance and short response time[END_REF], [START_REF] Snyder | Thermoelectric efficiency and compatibility[END_REF], among others. It is interesting to highlight (Pérez-Aparicio et al., 2016) again as one of the few that coupled the mechanic field in addition to the electric and thermal. The latter also covers some non-linear coupled magnetic effects, such as Ettinghausen or Nernst, that can be reflected in the theoretical formulation of (de [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF].

Discretization of electromagnetic fields

The multi-physics problems for SMs are each obtained by a particular solution procedure and a dedicated computer code. The most efficient approach to such multi-physics problems is to create a master code that allows exchanging information in real time and for every iteration among dedicated mono-physics simulation codes. The latter was accomplished successfully for fluid-structure interactions (e.g., Kassiotis et al., 2011a,b), by using the partitioned solution procedure [START_REF] Felippa | Partitioned analysis of coupled mechanical systems[END_REF] that allows reusing existing codes. This currently seems not to be possible when it comes to the coupling of mechanics and electromagnetism, due to the quite different philosophy in representing different fields and governing global equations on either side, with equilibrium versus global constitutive law.

Thus, the vast majority of works on coupled mechanics and electromagnetism turned toward extending the philosophy specialized for their particular field to all the other fields computations. The current approach by mechanics experts relies upon the standard isoparametric finite element approximations and the simplest case of using a scalar potential for electric and magnetic fields, which can accommodate, for instance, piezoelectric or piezomagnetic devices. This type of approach is further extended to higher order approximations, large displacement, and homogenization techniques always using the finite element approximation for all fields; see [START_REF] Gil | A new framework for large strain electromechanics based on convex multi-variable strain energies: Variational formulation and material characterisation[END_REF], [START_REF] Schröder | Two-scale homogenization of electromechanically coupled boundary value problems[END_REF] and [START_REF] Keip | Two-scale computational homogenization of electro-elasticity at finite strains[END_REF], among others.

A completely different approach is followed by experts in physics and electrical engineering, with a typical attempt by [START_REF] Tonti | The mathematical structure of classical and relativistic physics[END_REF] where the formalism most suitable for electromagnetism is generalized also to describe mechanics. The theoretical formulation of this kind does not consider the local fields, but rather the global quantities set within the corresponding mathematical structure often called De Rham's complex [START_REF] Bochev | Matching algorithms with physics: exact sequences of finite element spaces[END_REF]. The discrete approximation based on such formulation is different from standard finite element method [START_REF] Zienkiewicz | The Finite Element Method[END_REF], and is based on differential forms and exterior calculus.

The mathematical foundation of that kind of discretization is found in [START_REF] Nédélec | Mixed finite elements in R 3[END_REF] or [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] and one of the first applications appears in the literature under the name of edge elements, allowing to calculate variables on edges for a particular kind of problem in electromagnetism, the calculation of eddy-currents (Bossavit, 1988b). The step forward in the right direction is offered by Whitney's elements with a complete choice of unknowns where not only edges but also facets and volumes could allocate variables (Bossavit, 1988a). Since then, several works explored linear tetrahedron edge elements, [START_REF] Albanese | Finite element methods for the solution of 3d eddy current problems[END_REF] or [START_REF] Desbrun | Discrete differential forms for computational modeling[END_REF], prisms, pyramids and hexahedron linear elements [START_REF] Dular | Mixed finite elements associated with a collection of tetrahedra, hexahedra and prisms[END_REF], [START_REF] Gradinaru | Whitney elements on pyramids[END_REF] or [START_REF] Nédélec | Mixed finite elements in R 3[END_REF], and more recently some higher order elements [START_REF] Manges | Generation of tangential vector finite elements[END_REF], [START_REF] Webb | Hierarchal scalar and vector tetrahedra[END_REF], [START_REF] Wang | Curvilinear and higher order 'edge' finite elements in electromagnetic field computation[END_REF] or [START_REF] Yioultsis | The mystery and magic of whitney elements-an insight in their properties and construction[END_REF]. However, many of these elements cannot be fully invariant (as verified for hexahedral elements implementing Whitney's discretization) for distorted meshes as revealed by the basic 3D patch test [START_REF] Zienkiewicz | The Finite Element Method[END_REF], [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF]; the proof is given in [START_REF] Nédélec | Mixed finite elements in R 3[END_REF], where only hexahedral element produced by affine isoparametric transformation can be used; in other words, the element cannot be used for a degenerated parallelepiped. Such condition is always fulfilled for linear tetrahedron since the transformation is always affine unless the Jacobian is negative or zero.

The cell method proposed in [START_REF] Repetto | Global formulation of 3d magnetostatics using flux and gauged potentials[END_REF], [START_REF] Alotto | The Cell Method for Electrical Engineering and Multiphysics Problems: An Introduction[END_REF] or [START_REF] Tonti | The mathematical structure of classical and relativistic physics[END_REF] makes use of different geometrical entities in every element, such as vertices, edges, faces, and volumes in order to provide the best locations for unknown parameters. The advantage of such an approach is its ability to preserve tangential and normal continuity for electric and magnetic variables. This is achieved by cell-based discretization of the global governing equations transforming continuum operators into discrete operators over each cell. The choice of discrete approximation and element connectivity is based entirely on the complete topology of the problem and the connections of the different elements between them. The main disadvantage from the standpoint of coupling is the end product of discretization, which is equivalent to the global constitutive equations, which would be very difficult (if not impossible) to generalize to any other but linear case. Besides, this method complicates meshing since it requires a perfect match between the primal and the dual mesh when allocating the dual variables [START_REF] Ren | On the complementarity of dual formulations on dual meshes[END_REF].

For eventual coupling with mechanics, it is essential to note that the cell method approach is in contrast to the standard finite element philosophy. The latter relies upon the weak formulation of (equilibrium) equations to provide the corresponding discretization [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. Solving such weak form does not guarantee that the equation is ensured at every point of the domain, but rather in an average sense in a partitioned volume. All unknown variables, called degrees of freedom, are placed at nodes; inside the element, the discrete approximation is continuous and offers an important advantage of partition-of-unity to make sure the basic case (patch test) is verified. The finite element approximation typically does not offer the continuity of field derivatives across element boundaries, which was the main obstacle of extending the low-order finite element approximations to electric or magnetic fields.

Similarly, the cell method approach is not easily adapted to mechanics with several attempts to recast all fields within such unified framework typical for physics, such as [START_REF] Bossavit | Discrete magneto-elasticity: a geometrical approach[END_REF] or electromagnetics code like GetDP [START_REF] Dular | Getdp: A general environment for the treatment of discrete problems[END_REF]. Namely, despite similarities in theoretical formalisms, the mechanics primal and dual variables are second order tensors as opposed to electric, magnetic, and thermal fields where they remain vectors. Attempts to generalize the cell method approach to solve the problems with more complicated non-linear or rheological mechanics models have not been recorded, and furthermore, they seem very difficult to imagine. However, a few authors have used the exterior calculus formalism to adapt it to finite element method philosophy such as [START_REF] Angoshtari | Compatible-strain mixed finite element methods for 2d compressible nonlinear elasticity[END_REF] As for theoretical aspects of the problem, the electromagnetics community also evaluated a potential energy-type of variational formulation [START_REF] Penman | Complementary and dual energy finite element principles in magnetostatics[END_REF], [START_REF] Hammond | Calculation of inductance and capacitance by means of dual energy principles[END_REF], [START_REF] Hammond | Dual finite-element calculations for static electric and magnetic fields[END_REF], which is no different from the mechanical approach in, for instance, [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]. Alternatively, other variational formulations explored concern the complementary potential energy approach, which dual to the former variational form, as suggested in [START_REF] Razek | Computation of 3d electrostatic local fields and forces using complementarity of dual formulations. application for capacitance and torque calculations in micromotors[END_REF], [START_REF] Ren | A 3d vector potential formulation using edge element for electrostatic field computation[END_REF] or [START_REF] Golias | Constitutive inconsistency: rigorous solution of maxwell equations based on a dual approach[END_REF]. These latter works are again cast in the standard format of electro-engineering with the formulation trying to bracket the error estimates in term of the constitutive equations. Interestingly enough, one type of error estimates for mechanics proposed in [START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF] also considers the same approach where the potential energy based formulation gives the upper limit while the complementary gives the lower limit.

Micro-model for ferroelectric behavior

There is a noticeable trend in state-of-the-art smart materials applications, which is the miniaturization of the devices. The smaller they are, the more efficient they tend to be due to the control of the microstructure [START_REF] Rowe | Thermoelectrics handbook: macro to nano[END_REF]. Namely, the inclusions of thermoelectric materials in a crystal lattice structure help improving performance by reducing the thermal conductivity [START_REF] White | Energy-harvesting devices: Beyond the battery[END_REF]. Also, this miniaturization has allowed placing an increasing number of sensors and actuators on a structure to monitor and control the current state.

Numerically, the direct study of macro models is no longer valid to describe the behavior of a structure without a homogenization process. The inclusions of devices or materials in the microstructure requires a more detailed study of micro or even nanoscales. Additionally, the more realistic models for electromagnetics and mechanics induce localizations that are easier to capture at a smaller scale.

There is a wide number of references on homogenized solid finite element computation with localized failure and plasticity for mechanics such as [START_REF] Belytschko | A finite element with embedded localization zones[END_REF], [START_REF] Simo | An analysis of strong discontinuities induced by strainsoftening in rate-independent inelastic solids[END_REF], [START_REF] Armero | An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids[END_REF], [START_REF] Ibrahimbegovic | Combined hardening and softening constitutive model of plasticity: precursor to shear slip line failure[END_REF], [START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures: Part i: theoretical formulation and numerical implementation[END_REF], [START_REF] Saksala | Combined continuum damageembedded discontinuity model for explicit dynamic fracture analyses of quasi-brittle materials[END_REF], [START_REF] Do | Dynamics framework for 2d anisotropic continuum-discrete damage model for progressive localized failure of massive structures[END_REF], and [START_REF] Karavelić | Multi-surface plasticity model for concrete with 3d hardening/softening failure modes for tension, compression and shear[END_REF].

Different numerical representations of micro-models are available. In terms of numerical efficiency, a structure in terms of Voronoi cell representation with an irregular cohesive discrete lattice model is quite advantageous. This representation is best suited when trying to replicate the multi-scale behavior of a material composed of different micro-behaviors [START_REF] Ostoja-Starzewski | Lattice models in micromechanics[END_REF]. The cohesive links are one-dimensional finite elements, such as trusses or beams. This representation can replicate the behavior of an equivalent continuum structure [START_REF] Nikolić | Lattice element models and their peculiarities[END_REF].

In the literature, several papers can be found on the inclusion of embedded discontinuity models in lattice-type structures such as [START_REF] Nikolic | Rock mechanics model capable of representing initial heterogeneities and full set of 3d failure mechanisms[END_REF], [START_REF] Schlangen | Fracture simulations of concrete using lattice models: computational aspects[END_REF], [START_REF] Ibrahimbegovic | Microscale and mesoscale discrete models for dynamic fracture of structures built of brittle material[END_REF], [START_REF] Ibrahimbegovic | Embedded discontinuity finite element method for modeling of localized failure in heterogeneous materials with structured mesh: an alternative to extended finite element method[END_REF], and [START_REF] Bui | Enriched timoshenko beam finite element for modeling bending and shear failure of reinforced concrete frames[END_REF]. In terms of coupling mechanics with other physics, there are a few references availables, for instance, poroplastic media with fluids [START_REF] Hadzalic | Theoretical formulation and seamless discrete approximation for localized failure of saturated poro-plastic structure interacting with reservoir[END_REF] and [START_REF] Nikolic | Discrete element model for the analysis of fluid-saturated fractured poro-plastic medium based on sharp crack representation with embedded strong discontinuities[END_REF], or coupling with thermal field [START_REF] Ngo | Model for localized failure with thermoplastic coupling: theoretical formulation and ed-fem implementation[END_REF].

Regarding more elaborated constitutive models in terms of ferroelectrics, there are a hand-ful of references that describe this behavior for macro-scale in a similar fashion as plasticity, like [START_REF] Miehe | Variational principles in dissipative electro-magnetomechanics: A framework for the macro-modeling of functional materials[END_REF], [START_REF] Schröder | A thermodynamically consistent mesoscopic model for transversely isotropic ferroelectric ceramics in a coordinate-invariant setting[END_REF], [START_REF] Mcmeeking | A phenomenological multi-axial constitutive law for switching in polycrystalline ferroelectric ceramics[END_REF], [START_REF] Huber | A crystal plasticity model for ferroelectrics[END_REF][START_REF] Mcmeeking | On the potential energy of a piezoelectric inclusion and the criterion for ferroelectric switching[END_REF]. Other references focus on the behavior of the crystals at a micro-scale level such as [START_REF] Hwang | Ferroelectric/ferroelastic interactions and a polarization switching model[END_REF] or [START_REF] Chen | A micro-electro-mechanical model for polarization switching of ferroelectric materials[END_REF]. A complete description on the typical ferroelectric hysteresis loops can be found in [START_REF] Damjanovic | Hysteresis in piezoelectric and ferroelectric materials[END_REF], including the description of Debye models, such as [START_REF] Palma | Dissipative finite-element formulation applied to piezoelectric materials with the debye memory[END_REF]. References on multi-scale modeling of ferroelectrics and ferromagnetics can be found as well in the literature like [START_REF] Labusch | On the influence of ferroelectric polarization states on the magneto-electric coupling in two-phase composites[END_REF], [START_REF] Daniel | Homogenisation of magneto-elastic behaviour: from the grain to the macro scale[END_REF], [START_REF] Daniel | Reversible magneto-elastic behavior: A multiscale approach[END_REF], and [START_REF] Daniel | A multiscale model for magneto-elastic behaviour including hysteresis effects[END_REF], among others.

Methodology

The main focus of the opening part of the thesis, chapters 2 and 3, is to obtain a thermodynamically consistent formulation suitable for fully-coupled materials and plastic regime. In terms of completeness of the presented developments for thermomechanics and electromagnetic coupling, the recent ones that come closest to this part are the hypo-elastic approach in [START_REF] Pérez-Aparicio | Finite element analysis and material sensitivity of peltier thermoelectric cells coolers[END_REF], (Pérez-Aparicio et al., 2016), (Pérez-Aparicio et al., 2016), but they remain limited to elastic behavior. In contrast, a hyper-elastoplastic response, suitable to obtain the evolution equations for plastic variables, is provided in this work. First, a set of state variables and their corresponding dual ones are defined for each field (thermal, electric, magnetic, and mechanic) in establishing the governing equations based on conservation principles. Second, by settling a full quadratic form of the free-energy potential, the linear constitutive equations can be derived from it consistently, providing a hyper-elastic response instead of a hypo-elastic. Finally, the theoretical formulation is accompanied by the discrete approximation, based upon 3D Finite Elements with thermo-electro-magneto-mechanic degrees of freedom, along with simultaneous solution procedures of the weak form for all governing equations. The primary motivation for the second part of the thesis, chapter 4, is seeking the optimal formulation for coupled electromagnetic and thermomechanical fields, which will be amenable to the material and mechanics point-of-view and the finite element discrete approximations that can be combined with existing approaches and be integrated within the standard computer code architecture. This is especially the case when considering a currently critical issue of seeking to open up the path toward exploration of the engineering materials which can be heterogeneous or experience hardening and permanent damage with microscale point-of-view that proved very successful in mechanics (e.g., [START_REF] Balanis | Advanced engineering electromagnetics[END_REF] or [START_REF] Jackson | Classical electrodynamics[END_REF]). This approach is in contrast with the cell method, where variables are discretized employing global operators as opposed to the localization that requires phenomena like plasticity or damage. In comparison to the current tendency to introduce higher-order approximations, the discrete approximations proposed in this work are low-order (linear), which makes them more suitable for heterogeneous materials and inelasticity regime where the solution, in general, is not smooth. Besides, we recast the novel approach to electrostatics and magnetostatics, where the constitutive equations are an explicit part of the formulation, easily generalizable to more complex behavior like ferroelectricity. We also recast the discrete approximation in the manner that can be combined with finite element codes, by exploiting the partition-of-unity property SECTION 1.5. Outline of the thesis of the finite element discretization.

In the final part, the objective is to obtain a numerical model for beams to replicate the micro behavior of a ferroelectric grain. This model also implements a model for nonlinear viscoplasticity and localized failure for the exponential softening regime. The remanent polarization is described, following [START_REF] Hwang | Ferroelectric/ferroelastic interactions and a polarization switching model[END_REF], as a constant value added on top of the linear electric displacement. A polarization multiplier determines the evolution of the remanent polarization and strain. Each beam pictures a grain of the material. Therefore, the lattice model produces the macro response of a heterogeneous ferroelectric model. This response can be obtained by averaging the internal values of the system, or through a process of homogenization where the stiffness terms are averaged [START_REF] Rukavina | Ed-fem multi-scale computation procedure for localized failure[END_REF].

All the numerical implementations and computations are carried out with the computer code FEAP, a Finite Element Analysis Program [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. This code has been developed at the University of Berkeley, California, by Prof. R.L. Taylor in Fortran. Among the different parts of this thesis, other computer codes have been used, such as Mathematica (Wolfram), Matlab (Mathworks), or Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and postprocessing facilities[END_REF].

Outline of the thesis

The thesis is divided into six chapters, including this one. In Chapter 2, the macroscale formulation is derived for simulating fully-coupled conductors and insulators, as well as the corresponding finite element implementation for the linear-elastic case. A basic model for plasticity is introduced into the previous thermodynamically consistent formulations in Chapter 3, where all the fields play a role in the model criteria, and they all get affected by the plastic regime. Chapter 4 introduces a discretization based on differential forms and exterior calculus that perfectly match Maxwell's equations. A micro-scale lattice-solid model is presented in Chapter 5 for ferroelectric models of ferroelectrics, in conjunction with domain switching and crack propagation. Finally, some concluding remarks, along with possible future works or extensions, are summarized in Chapter 6.

2

Multiphysics coupling formulation

Two consistent thermodynamic-based theoretical frameworks for insulators and conductors, capable of coupling linear interactions of thermal, electric, magnetic, and mechanic fields are presented. The complete set of governing equations for all fields is obtained from the conservation principles. The second principle of thermodynamics is taken into account to introduce the irreversible Joule heating for conductors. A full-interaction for each corresponding dual variable in terms of the defined set of state variables is defined through Helmholtz free-energy potential, which provides that the constitutive equations for corresponding dual variables can be derived consistently. The implementations are carried out in a research version of the research computer code FEAP by using isoparametric 3D hexahedra solid elements. Several numerical simulations are shown in order to illustrate the proposed model and formulation capabilities.

Formulation

Kinematic equations

Four fields are considered in formulating the thermo-electro-magneto-mechanic coupling: displacement u, temperature T , electric potential V and magnetic scalar potential ϕ. The state variables are obtained as the corresponding gradients of these fields. The set of resulting "kinematic" equations can be written as:

ε = 1 2 ∇ ⊗ u + (∇ ⊗ u) T E = -∇V H = -∇ϕ (2.1)
where ε is the strain tensor, E is the electric field, H is the magnetic field, while ∇ = [∂/∂x 1 ∂/∂x 2 ∂/∂x 3 ] T is a convenient notation for the nabla operator of partial derivatives.

With the hypothesis of small displacement gradient theory, we are limited here to strains defined in terms of the symmetric part of the displacement gradient. Hence the operator ∇ sy can be used to define the strain tensor, which can also be written in Voigt notation as:

ε = ∇ sy u ; ∇ sy =                   ∂ ∂x 1 0 0 0 ∂ ∂x 2 0 0 0 ∂ ∂x 3 ∂ ∂x 2 ∂ ∂x 1 0 0 ∂ ∂x 3 ∂ ∂x 2 ∂ ∂x 3 0 ∂ ∂x 1                   (2.2)

Conservation principles

A foreword on electric and magnetic constitutive behavior

One has to deal with two diametrically different constitutive behaviors for electric field in terms of conductors versus insulators. The key difference between those two types of behavior is that the electrons for the conductors are free to move, whereas for the insulators the electrons remain bounded to their nucleus [START_REF] Balanis | Advanced engineering electromagnetics[END_REF], [START_REF] Jackson | Classical electrodynamics[END_REF].

In conductors, the field of interest is the electric flux J, which measures the flow of electrons through a surface. For the simplest linear relation between the electric field and the electric flux defined for conductors through the conductivity coefficient γ, we can write:

J = γE (2.3)
which indicates that the electric flux will clearly go in the same direction as the electric field.

For perfect conductors, the conductivity coefficient γ → ∞ whereas for perfect insulators γ → 0. Most of all real materials fall in-between these two extremes. Namely, an insulator like glass has free electric flux, although it remains many orders of magnitude lower than the one in a conductor such as copper; for instance, the two materials we just compared have a relation of conductivities of the order γ con ≈ 10 20 γ ins [START_REF] Mitchell | An introduction to materials engineering and science for chemical and materials engineers[END_REF]. Thus, the flux in an insulator material is simply considered as zero.

In fact, insulators react to electric fields by orienting dipoles, deformed atoms or molecules with two opposite charges held in place by molecular forces. The electric charge is the physical property of matter that causes it to experience a force when placed within an electromagnetic field. The decomposition of the density of electric charge ρ q = ρ f q + ρ b q into free and bound terms is introduced. Thus, the latter are not free to move as it happens with the former in conductors. The orientation of dipoles is measured by the polarization vector P . Dielectric materials are a particular case of insulators where the dual variable electric displacement D is proportional to E. Therefore, the electric displacement definition can be written as:

D = 0 E + P = 0 (I + χ e ) E = E (2.4)
with 0 as the permittivity of the vacuum, as the permittivity of the material, and I as the second-order identity tensor.

In Figure 2.1, we illustrate the behavior of conductors and dielectric materials with respect to an applied electric field. For conductors, the electrons move through the domain in the direction of E. For the insulator, the dipoles inside the domain are also oriented in the same direction, but net charges appear only at the boundaries, due to cancellation of opposite charges inside the domain [START_REF] Hammond | Electromagnetism for engineers: an introductory course[END_REF]. In contrast to the complexity of electric behavior placed between two extremes of conductors and insulators, the magnetic constitutive behavior is simpler and it does not have a conductor analogy. For this thesis, only diamagnetic materials are taken into account, where the magnetization M is proportional to H. Thus, we define the dual magnetic induction B through linear constitutive relation:

B = µ 0 (H + M) = µ 0 (I + χ h ) H = µH (2.5)
with µ 0 as the permeability of vacuum and µ as the permeability of the material.

Conservation of electric and magnetic flux, and free electric charge

The traditional description of the electromagnetic phenomena is defined by four laws first proposed by Faraday, Ampère and Gauss. These equations feature E and H, and their dual fields, electric displacement D and magnetic induction B. They can jointly be cast in terms of four differential equations, known as Maxwell's equations for electromagnetism [START_REF] Balanis | Advanced engineering electromagnetics[END_REF]. The latter provides, from the macroscopic point-of-view, an adequate description for many electromagnetic phenomena observations on how these fields evolve in the presence of the free electric flux J and free electric charge density ρ f q :

∇ × E = -Ḃ ;

∇ × H = J + Ḋ ; ∇ • D = ρ f q ; ∇ • B = 0 (2.6)
These equations are enough to describe a pure electromagnetic problem. However, if mechanics are involved, the definition of Lorentz's force must be added to the conservation of linear momentum. One can also recover from (2.6) 2 and (2.6) 3 the condition on the conservation of electric charge that constraints ρ f q and J for coupled electromagnetic case according to:

ρf q + ∇ • J = 0 (2.7)
For insulators, we consider quasi-static case for the sake of simplicity, with negligible rate of change of electric and magnetic fields and time derivatives Ḃ, Ḋ and ρf q equal to zero. Therefore, we uncouple electrics and magnetics obtaining two separate problems, one for i) electrostatics:

∇ × E = 0 ;

∇ • D = ρ f q (2.8)
and another for ii) magnetostatics:

∇ × H = J ; ∇ • B = 0 ; ∇ • J = 0 (2.9)
Thanks to this kind of weak coupling, each of these problems can be handled independently from one another. The equations (2.8) 1 and, if J = 0, (2.9) 1 are automatically verified due to the convenience of the kinematics assumptions in (2.1).

For conductors, the only focus in this thesis concerns the electric field and its coupling with thermal and mechanic field. Thus, magnetic field is not going to be calculated. After combining (2.6) 2 and (2.6) 3 as described before, the set of equations to solve is:

∇ × E = 0 ; ∇ • J = 0 (2.10)
In other words, the introduction of the conservation of electric charge automatically corroborates (2.6) 2 and (2.6) 3 under the hypotheses of quasi-static case. Again, (2.10) 1 is verified with the current kinematics assumptions.

Conservation of momentum

The classical approach to conservation of momentum (e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]) is now generalized to account for the electromagnetic field. First, by enforcing the angular momentum conservation we obtain the symmetry of stress tensor (e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]. Second, by postulating the linear momentum conservation principle for a domain, and in the limit case of the domain shrinking to a point, we obtain the local form of momentum conservation equation (de [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF]:

ρ m ü ṗm + Ṡp c 2 ṗeh = (σ C + σ M ) ∇ + b m (2.11)
where ρ m is the mass density, p m the mechanical linear momentum, S p = E × H the Poynting vector, p eh the electromagnetic linear momentum, σ C the Cauchy stress tensor, σ M the Maxwell electromagnetic stress tensor, b m the volume forces and c the speed of light, related to the permeability and permittivity of the vacuum through the following expression:

c 2 = 1 0 µ 0 (2.12)
This expression can be changed to transform the Maxwell stress tensor and the electromagnetic linear momentum into a single corresponding force. The assumed behavior for the materials in this chapter is dielectric and diamagnetic, although the following development is made in the most general manner. The time derivative of the electromagnetic momentum (2.11) can be expanded into:

0 µ 0 ∂ (E × H) ∂t = ∂ (D × B) ∂t -0 µ 0 ∂ (E × M) ∂t - ∂ (P × B) ∂t (2.13)
By using the Maxwell's equations (2.6) and the following identity:

(∇ × u) × v = [v ⊗ u -(u • v ) I] ∇ -u (∇ • v ) + (∇ ⊗ v ) u (2.14)
it is possible to obtain an alternative expression for the first term of the right hand side:

∂ (D × B) ∂t = [D ⊗ E + B ⊗ H -(D • E + B • H) I] ∇ + (∇ ⊗ D) E + (∇ ⊗ B) H -ρ f q E -J × B (2.15)
The last equation can also be transformed by using the following:

(∇ ⊗ D) E = P • E + 1 2 0 E • E I ∇ -(∇ ⊗ E) P (∇ ⊗ B) H = 1 2 µ 0 (B • B) I ∇ -(∇ ⊗ B) M (2.16)
By introducing all these identities into (2.15), we can write:

∂ (D × B) ∂t = D ⊗ E + B ⊗ H - I 2 ( 0 E • E + µ 0 B • B -2B • M) ∇ -(∇ ⊗ E) P -(∇ ⊗ B) M -ρ f q E -J × B
(2.17)

With these results at hand, we return to (2.13), with three terms set apart:

ṗeh = σ M ∇ + b eh (2.18)
where by identification:

b eh • • = (∇ ⊗ E) P + (∇ ⊗ B) M + 0 µ 0 ∂ (E × M) ∂t + ∂ (P × B) ∂t + ρ f q E + J × B ; σ M • • = D ⊗ E + B ⊗ H - I 2 (D • E + H • B -2B • M) (2.
19) The electromagnetic force b eh , which accounts for the distributed electric and magnetic force that can arise due to the possible presence of a non-uniform electrical or magnetic field in the material [START_REF] Ferrari | Thermodynamic formulation of the constitutive equations for solids and fluids[END_REF]. This variable accounts for the Lorentz's force in the last two terms and some others coming from the polarization and magnetization of the material. Finally, the conservation of momentum (2.11) can be expressed in terms of b eh simply as:

ρ m ü = σ C ∇ + b eh + b m (2.20)
Both interpretations in (2.11) and (2.20) can be used indistinctly. To illustrate this feature, (2.11) will be used for conductors and (2.20) for insulators. An additional hypothesis is assumed concerning all time derivatives of electric and magnetic variables. These rate terms are fast enough that we can consider that electromagnetism happens instantaneously, compared to thermal or mechanic field. Thus, the derivatives will be considered null.

To unify formulations for both kind of materials, the stress tensor is σ = σ C + σ M for conductors and σ = σ C for insulators; the body force is b = b m for conductors and b = b eh + b m .

Conservation of energy

The global form of the first principle of thermodynamics can be stated in a closed domain Ω:

∂ ∂t E = W m + W eh + Q (2.21)
where E is the total energy, W m is the mechanical power, W eh is the electromagnetic power and Q is the total heat supplied to that domain. The total energy can be split into potential P and kinetic K energy terms:

E = P + K = ˆΩ e (ε, s, D, B) dΩ + 1 2 ˆΩ ρ m u • u dΩ (2.22)
where e is the scalar potential of the internal energy density that depends on the state variables.

The list of state variables is defined in Table 2.1, along with their corresponding dual variables (with s as the entropy per unit volume). The source of mechanical power inserted within the particular domain Ω can be written as:

Fields

W m • • = ˆΩ b • u dΩ + ˆΓ t n • u dΓ = ˆΩ b • u dΩ + ˆΩ ∇ • (σ u) dΩ (2.23)
where n is unit normal vector. The last term is defined from the boundary traction vector t n by the Cauchy principle t n = σ n; this vector has been transformed into the corresponding volume integral with the divergence theorem. A similar integral transformation can be done for the boundary term of the electromagnetic power source (see [START_REF] Balanis | Advanced engineering electromagnetics[END_REF], and it is defined as:

W eh • • = -ˆΓ (E × H) • n dΓ = -ˆΩ ∇ • (E × H) dΩ (2.24)
We note in passing that the negative sign is in agreement with the above expression representation of total power exiting the volume Ω bounded by the surface Γ. Similar transformation can finally be made for the heat power source that stems from the outgoing heat flux q, along with the heat source r , which results in:

Q • • = ˆΩ r dΩ -ˆΓ q • n dΓ = ˆΩ r dΩ -ˆΩ ∇ • q dΩ
(2.25)

In the limit case of Ω shrinking to a point, we obtain from (2.21) the local form of the first principle, which can be written as:

ė (ε, s, D, B) + ρ m ü • u = b • u + ∇ • (σ u) + r -∇ • q -∇ • (E × H) (2.26)
Furthermore, by substituting the kinematic equations in (2.1) and the equation of motion in (2.20), along with the identity ∇ • (σ u) = (σ ∇) • u + σ • (∇ ⊗ u), we can obtain the reduced form of the first principle or energy conservation:

ė (ε, s, D, B) = σ • ε + r -∇ • q -∇ • (E × H) (2.27)
The final ingredients pertain to provide definition of electric and heat fluxes through the generalized form of Ohm's and Fourier's laws, as suggested in [START_REF] Anatychuk | Modern Thermodynamic Theory of Thermoelectricity[END_REF] and (Pérez-Aparicio et al., 2016):

J = γ E -γ α ∇T ; q = -κ ∇T + Π J (2.28)
where α is the Seebeck coefficient, κ is the thermal conductivity and Π = T α is the Peltier coefficient. For insulators, α and γ are assumed to be zero.

By scalar multiplication of the first two Maxwell's equations of (2.6) with H and E respectively, and by exploiting the following identity:

∇ • (E × H) = H • (∇ × E) -E • (∇ × H) (2.29)
the reduced form of energy conservation (2.27) can be recast in an equivalent format:

ė (ε, s, D, B) = σ • ε + r -∇ • q + E • Ḋ + H • Ḃ (2.30)
The second principle of thermodynamics imposes that the rate of increase of entropy Ṡ should never be smaller than the amount of heat divided by the absolute temperature. For a particular domain Ω we can write:

Ṡ ≥ Q T ; S = ˆΩ s dΩ (2.31)
In the limit case of shrinking the domain to a point, we can obtain the local form of the second principle. In the simplest case of a rigid conductor (corresponds to neglecting all fields but the thermal), the second principle provides the proper definition of dissipation by conduction D c [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], which always remains positive as long as Fourier's law applies:

- 1 T q • ∇T D c • • = T ṡ -(r -∇ • q) ≥ 0 (2.32)
The second principle, combined with the result of the first principle in (2.30) can be used to define the local dissipation, which always remains non-negative:

D • • = T ṡ -ė (ε, s, D, B) + σ • ε + J • E + E • Ḋ + H • Ḃ ≥ 0 (2.33)
In the result above we dropped the dissipation by conduction. We can introduce the free energy potential by appealing to the Legendre transformation (e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], which allows to exchange the roles between the state variables and their dual, s, D, B versus T, E, H:

ψ (ε, T, E, H) = e (ε, s, D, B) -T s -E • D -H • B (2.34)
Deriving with respect to time the last expression we can obtain:

ė = ∂ψ ∂ε • ε + ∂ψ ∂T Ṫ + ∂ψ ∂E • Ė + ∂ψ ∂H • Ḣ + Ṫ s + T ṡ + E • Ḋ + Ė • D + H • Ḃ + Ḣ • B (2.35)
With this result on hand, the dissipation (2.33) can be expressed as:

0 ≤ D = σ - ∂ψ ∂ε • ε -s + ∂ψ ∂T Ṫ -D + ∂ψ ∂E • Ė -B + ∂ψ ∂H • Ḣ + J • E (2.36)
The last term is often referred to as Joule's dissipation or heating, which can be written as:

D J = J • E (2.37)
The dissipation inequality will become an equality providing the set of constitutive equations to be defined in agreement with the chosen free energy potential. Here, we choose a quadratic form of the free energy potential, which can be written as follows for conductors:

ψ (ε, T, E, H) = ψ m + ψ t + ψ e + ψ h + ψ mt + ψ me + ψ mh ; ψ m = 1 2 ε • Cε ; ψ mt = -β • (T -T 0 ) ε ; ψ t = ρ m c p (T -T 0 ) -T ln T T 0 ; ψ me = (E ⊗ E) • ε - 1 2 (E • E) tr (ε)
;

ψ e = - 1 2 E • E ; ψ mh = µ (H ⊗ H) • ε - 1 2 (H • H) tr (ε) ψ h = - 1 2 H • µH ; (2.38)
where C is the elasticity tensor, c p the specific heat, and β = Cα T I the thermal isotropic stress tensor (with α T as the expansion coefficient). Analogously, a quadratic free energy potential for insulators:

ψ (ε, T, E, H) = ψ m + ψ t + ψ e + ψ h + ψ mt + ψ me + ψ mh + ψ te + ψ th + ψ eh ; ψ m = 1 2 ε • Cε ; ψ me = -ε • e e E ; ψ t = ρ m c p (T -T 0 ) -T ln T T 0 ; ψ mh = -ε • e h H ; ψ e = - 1 2 E • E ; ψ te = -(T -T 0 ) π e • E ; ψ h = - 1 2 H • µH ; ψ th = -(T -T 0 ) π h • H ; ψ mt = -(T -T 0 ) β • ε ; ψ eh = -E • υH (2.39)
with e e as the piezoelectric tensor, e h as the piezomagnetic tensor, π e as the pyroelectric vector, π h as the pyromagnetic vector, and υ as the magnetoelectric tensor. The corresponding hyperelastic constitutive equations can easily be obtained by derivatives of such a potential with respect to the state variables. For conductors:

σ • • = ∂ψ ∂ε = Cε -β (T -T 0 ) + E ⊗ E + µH ⊗ H - I 2 ( E • E + µ H • H) ; D • • = - ∂ψ ∂E = {[1 + tr (ε)] E -εE} ; B • • = - ∂ψ ∂H = µ {[1 + tr (ε)] H -εH} ; s • • = - ∂ψ ∂T = ρ m c p ln T T 0 + β • ε (2.40)
For insulators:

σ • • = ∂ψ ∂ε = Cε -β (T -T 0 ) -e e E -e h H s • • = - ∂ψ ∂T = β • ε + ρ m c p ln T T 0 + π e • E + π h • H D • • = - ∂ψ ∂E = e e ε + π e (T -T 0 ) + E + υH B • • = - ∂ψ ∂H = e h ε + π h (T -T 0 ) + υE + µH (2.41)
Again, with these results on hand, we can write from (2.32) the heat equation for an elastic case as:

T ṡ = r -∇ • q (2.42)
where all terms in the expression for dissipation (2.36) are zero due to the constitutive equation definitions except for the Joule heating. Furthermore, by using the entropy s from (2.41) or (2.40), we can rewrite the heat equation for conductors:

ρ m c p Ṫ = r -∇ • q -T β • ε + D J (2.43)
And for insulators, where J is neglected:

ρ m c p Ṫ = r -∇ • q -T β • ε + π e • Ė + π h • Ḣ (2.44)
In summary, the strong form of the equations to be solved fopr conductors regroups the results written in (2.11), (2.10) 2 , (2.6) 4 and (2.43), here restated in tensor notation:

ρ m ü = σ∇ + b ∇ • J = 0 ∇ • B = 0 ρ m c p Ṫ = r -∇ • q -T β • ε + D J (2.45)
And for insulators, (2.20), (2.8) 2 , (2.6) 4 and (2.44):

ρ m ü = σ∇ + b ∇ • D = ρ f q ∇ • B = 0 ρ m c p Ṫ = r -∇ • q -T β • ε + π e • Ė + π h • Ḣ (2.46)

Finite element implementation

In this section we present the details of the discrete approximation constructed by the the finite element method, as a particular case of the Galerkin method. The starting point is provided by the weak form of the conservation equations in (2.45) which can be stated in tensor notation as follows:

-ˆΩ ∇ sy δu • σ + δu • (ρ m ü -b) dΩ + ˆΓ δu • tc dΓ = 0 ; ˆΩ ∇δV • J dΩ -ˆΓ δV j dΓ = 0 ; ˆΩ ∇δT • q -δT ρ m c p Ṫ -r + T β • ε -J • E dΩ -ˆΓ δT q dΓ = 0 (2.47)
where δ• denotes a virtual field or variation. It can be noticed that the conservation equation for the magnetic field has been dropped as already stated above. We can readily obtain the discrete approximations for all the fields, along with their space and time derivative by appealing to separation of variables:

u ≈ N b a U b ; V ≈ N b a V b ; ϕ ≈ N b a ϕ b ; T ≈ N b a T b ; δu ≈ N a w U a ; δV ≈ N a w V a ; δϕ ≈ N a w ϕ a ; δT ≈ N a w T a ; ∇ sy u ≈ B sy b a U b ; ∇V ≈ B b a V b ; ∇ϕ ≈ B b a ϕ b ; ∇T ≈ B b a T b ; ∇ sy δu ≈ B sy a w U a ; ∇δV ≈ B a w V a ; ∇δϕ ≈ B a w ϕ a ; ∇δT ≈ B a w T a ; ∇ sy u ≈ B sy b ȧU b ; ∇ V ≈ B b ȧV b ; ∇ φ ≈ B b ȧϕ b ; Ṫ ≈ N b ȧT b ; ü ≈ N b äU b ; a = 1, 2, . . . , 8 b = 1, 2, . . . , 8
(2.48) where a represent the nodal values of different fields (yet called degrees of freedom), whereas w represent the nodal values of their variations. In last expression, N a denotes the standard isoparametric shape function for node a (e.g., [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], with their gradients gathered in matrix form as:

B a = ∇N a ; B sy a = ∇ sy N a (2.49)
By introducing the corresponding finite element approximations into the weak form of the conservation equations, and switching from tensors to matrices by means of the Voigt notation, we can finally obtain:

- ˆΩ B sy a w U a T σ + N a w U a T ρ m N b äU b -b dΩ + ˆΓ N a w U a T tc dΓ = 0 ; ˆΩ B a w V a T J dΩ -ˆΓ N a w V a j dΓ = 0 ; ˆΩ B a w T a T q -N a w T a ρ m c p N b ȧT b -r + N b a T b β T B sy b ȧU b + J T B b a V b dΩ -ˆΓ N a w T a q dΓ = 0 (2.50)
where we introduce the discrete approximations of dual variables and fluxes:

                 σ = C (6×6) B sy b (6×3) a U b (3×1) -β N b a T b -T 0 + 0 B ⊗ b a V b 2 - 1 2 0 B b a V b T B b a V b ; J = -γ B b a V b -γ α B b a T b ; q = -κ B b a T b + N b a T b α J (2.51)
featuring 1 = (1, 1, 1, 0, 0, 0) T and the subscript ⊗ as the tensor multiplication into Voigt notation; in this case

B ⊗ b = B b ⊗ B b : B ⊗ b =         N b,1 N b,1 N b,2 N b,2 N b,3 N b,3 N b,1 N b,2 N b,2 N b,3 N b,1 N b,3         (2.52)
with N b,i is the derivative of the shape function in the direction i .

By considering that the nodal values of virtual field (∀w) can be picked arbitrarily, it is possible to obtain from (2.50) the final set of non-linear residual equations that need to be solved:

R U a = -ˆΩ B sy a T σ + N a ρ m N b äU b -b dΩ + ˆΓ N a tc dΓ ; R V a = ˆΩ B T a J dΩ -ˆΓ N a j dΓ ; R T a = ˆΩ B T a q -N a ρ m c p N b ȧT b -r + N b a T b β T B sy b ȧU b + J T B b a V b dΩ -ˆΓ N a q dΓ (2.53)
A similar process can be followed for insulators obtaining the following residuals:

R U a = -ˆΩ B sy a T σ + N a ρ m N b äU b -b dΩ + ˆΓ N a tc dΓ ; R V a = ˆΩ B T a D + N a ρ f q dΩ -ˆΓ N a D dΓ ; R ϕ a = ˆΩ B T a B dΩ -ˆΓ N a B dΓ ; R T a = ˆΩ B T a q -N a ρ m c p N b ȧT b -r + N b a T b β T B sy b ȧU b -π T e B b ȧV b -π T h B b ȧϕ b dΩ -ˆΓ N a q dΓ (2.54)
where we introduced the discrete approximations into the dual variables:

σ = CB sy b a U b -β N b a T b -T 0 + e T e B b a V b + e T h B b a ϕ b ; D = e e B sy b a U b + π e N b a T b -T 0 -B b a V b -υB b a ϕ b ; B = e h B sy b a U b + π h N b a T b -T 0 -υB b a V b -µB b a ϕ b (2.55)

Time Discretization

Global solution step by Newmark scheme

The Newmark scheme is used for the time discretization of the global solution step. This scheme requires two parameters γ and β that will determine the numerical damping and order of the scheme (e.g., [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]:

a b,n+1 = a b,n + ∆t ȧb,n + ∆t 2 1 2 -β äb,n + β äb,n+1 ; ȧb,n+1 = ȧb,n + ∆t [(1 -γ) äb,n + γ äb,n+1 ] (2.56)
where we denoted the time step as ∆t = t n+1t n . These equations, often referred to as Newmark equations, are accompanied by the residual equations enforcing the zero value at time step n + 1, which can be explicitly written for conductors as:

R U,(i) a,n+1 = -ˆΩ B sy a T σ (i) n+1 + N a ρ m N b äU,(i) b,n+1 -b n+1 dΩ + ˆΓ N a tc n+1 dΓ ; R V,(i) a,n+1 = ˆΩ B T a J (i) n+1 dΩ -ˆΓ N a jn+1 dΓ ; R T,(i) a,n+1 = ˆΩ B T a q (i) n+1 -N a ρ m c p N b ȧT,(i) b,n+1 -r n+1 + N b a T,(i) b,n+1 β T B sy b ȧU,(i) b,n+1 +J (i),T n+1 B b a V,(i) b,n+1 dΩ -ˆΓ N a qn+1 dΓ (2.57)
and for insulators as:

R U,(i) a,n+1 = -ˆΩ B sy a T σ (i) n+1 + N a ρ m N b äU,(i) b,n+1 -b (i) n+1 dΩ + ˆΓ N a tc n+1 dΓ ; R V,(i) a,n+1 = ˆΩ B T a D (i) n+1 + N a ρ f q,n+1 dΩ -ˆΓ N a Dn+1 dΓ ; R ϕ,(i) a,n+1 = ˆΩ B T a B (i) n+1 dΩ -ˆΓ N a Bn+1 dΓ ; R T,(i) a,n+1 = ˆΩ B T a q (i) n+1 -N a ρ m c p N b ȧT,(i) b,n+1 -r n+1 + N b a T,(i) b,n+1 β T B sy b ȧU,(i) b,n+1 -π T e B b ȧV,(i) b,n+1 -π T h B b ȧϕ,(i) b,n+1 dΩ -ˆΓ N a qn+1 dΓ (2.58)
where the corresponding values for velocities ȧ(i) n+1 and accelerations ä(i) n+1 are obtained by recasting the result (2.56) resulting with:

ȧ(i) b,n+1 = ȧb,n + ∆t (1 -γ) äb,n + γ ä(i) b,n+1 ; ä(i) b,n+1 = a (i) b,n+1 -a b,n -∆t ȧb,n β ∆t 2 + 1 - 1 2β äb,n (2.59)
Thus, the time stepping scheme of this kind will finally render the set of nonlinear algebraic equations. To solve such a nonlinear problem, we use Newton's iterative method where at each iteration (i + 1) we perform the consistent linearization of residual leading to:

R (i+1) a,n+1 = 0 ⇒ R (i) a,n+1 + ∂R a,n+1 ∂a b,n+1 (i) ∆a (i) b,n+1 = 0 (2.60)
where ∆a b are iterative contributions to nodal values of temperature, electric field, along with displacements, velocities and accelerations. At each iterative sweep, we can then perform the corresponding state variable updates according to:

a (i+1) b,n+1 = a (i) b,n+1 + ∆a (i) b,n+1
(2.61)

In the first iteration within each time step we will take the starting guess equal to the converged value at the previous step:

a (0) b,n+1 = a b,n (2.62)
The mechanics part of the residual vector at particular iteration can be further compressed, reducing it to the form presented explicitly in (2.57) or (2.58). Namely, the second term in (2.60) can be reduced to so-called effective tangent stiffness for mechanical part, which is directly used to compute the iterative contributions to displacement increments:

S U,(i) ab = - ∂R U a ∂a b (i) - ∂R U a ∂ ȧb ∂ ȧb ∂a b (i) - ∂R U a ∂ä b ∂ä b ∂a b (i) (2.63)
where time-step subscript n + 1 was dropped to simplify notation. By exploiting the relations between the nodal displacements and its first and second derivatives provided by the Newmark scheme, we can provide the closed form final linearized problem to be solved. More precisely, in view of the Newmark result for constructing discrete approximations for nodal velocities and accelerations in (2.56), we can write:

∂ä b ∂a b = 1 β ∆t 2 ; ∂ ȧb ∂a b = ∂ ȧb ∂ä b ∂ä b ∂a b = γ β ∆t (2.64)
Thus, the final form of the tangent operator for mechanics part can now be written as:

S U,(i) ab = K (i) ab + γ β ∆t D (i) ab + 1 β ∆t 2 M (i) ab (2.65) where K ab = - ∂R a ∂a b is the elastic tangent stiffness matrix, D ab = - ∂R a
∂ ȧb is the damping matrix and M ab = -∂R a ∂ä b is the mass matrix.

Given the solution for displacement increment at iteration (i ), we proceed to compute the displacement updates a (i+1) n+1 by using the result in (2.61).

Numerical examples

Insulators

In this section, we present the results of numerical simulations for coupled materials with conductor or insulator electric behavior. Several applications are studied as examples of the constructed element capabilities. All computations are performed by a research version of the well-known computer code FEAP [START_REF] Zienkiewicz | The Finite Element Method[END_REF].

The insulator materials are tested first. They present a preferential direction due to their internal structure [START_REF] Smith | Smart material systems: model development[END_REF]. As a result, a transverse isotropic model is taken into consideration, model that leads into vectors for π e and π h , diagonal tensors for µ, , υ and β, the last four with a component different than the others (the polarized direction). For the remaining coefficients, a simplification from a full form tensor into a symmetric one with many null entries is necessary; the non-zero coefficients are related among them. All strain related coefficients are expressed in what follows in Voigt notation: β, e e , e h , C.

Introducing the notation of a second order diagonal tensor:

diag (a 1 , a 2 , a 3 ) =   a 1 0 0 0 a 2 0 0 0 a 3   (2.66)
unless otherwise said, the properties extracted from [START_REF] Ferrari | Thermodynamic formulation of the constitutive equations for solids and fluids[END_REF], [START_REF] Ramirez | Free vibration response of two-dimensional magnetoelectro-elastic laminated plates[END_REF] and (Pérez-Aparicio et al., 2016) for the materials studied in this section are presented here: (11.2, 11.2, 12.6) × 10 -9 F/m µ = diag (5, 5, 10) × 10 -6 H/m υ = diag (5.37, 5.37, 2737.5) × 10 -12 s/m β = [1.67, 1.67, 1.96, 0, 0, 0] T MPa/K π e = [58.3, 58.3, 58.3

ρ m = 5700 kg/m 3 c p = 434 J/kg K κ = 2.61 W/m K T 0 = 293 K = diag
] T × 10 -5 C/m 2 K π h = [5, 5, 5] T × 10 -2 kg/s 2 A K e e =   0 0 0 0 0 11.6 0 0 0 0 11.6 0 -4.4 -4.4 18.6 0 0 0   C m 2 e h =
  0 0 0 0 0 5.5 0 0 0 0 5.5 0 5.8 5.8 7 0 0 0

  T C =        
116 77 78 0 0 0 77 116 78 0 0 0 78 78 162 0 0 0 0 0 0 89 0 0 0 0 0 0 86 0 0 0 0 0 0 86

       

GPa

(2.67) These properties correspond to different materials: PZT, Terfenol-D, BaTiO 3 or CoFe 2 O 4 ; however, in this chapter, they artificially belong to the same material just for calculation purposes.

Piezoelectricity

A piezoelectric material is considered having null coefficients and tensors in (2.41) except for C, e e and . Therefore, only electric and mechanic field variables matter in this section. This material works two ways: as actuator, which can induce movements when an electric field is applied, or as generator, for which when a displacement or a force is applied, an electric potential distribution is generated.

V = 10 V 1 x 3 x 1 x 2 3 3 x 2 x 1 Figure 2
.2: Sketch of the boundary conditions considered for the actuator validation example. Three planes of symmetry considered in x 1 = x 2 = x 3 = 0; electric potential prescriptions at the top and bottom planes. Measures in mm.

In Figure 2.2, the boundary conditions for the first validation example can be seen with the piezoelectric working as actuator. Regarding the mechanical field, symmetry conditions have been taken into account so that an eighth of a simple box geometry with dimensions 6 × 6 × 2 mm is represented. As for the electric field, V = 10 V at the top and ground at the x 3 symmetry plane have been set.

1.31E-09 2.63E-09 3.94E-09 5.25E-09 6.57E-09 0.00E+00 7.88E-09 Displacement in x direction -4.66E-09 -3.73E-09 -2.79E-09 -1.86E-09 -9.31E-10 -5.59E-09 0.00E+00 Displacement in x direction 1.67E+00 3.33E+00 5.00E+00 6.67E+00 8.33E+00 0.00E+00 1.00E+01 Electric Potential (V) 1 3
Figure 2.3: Displacement in directions x 1 and x 3 and electric potential distributions for the actuator validation piezoelectric case.

The voltage difference between the top and the bottom faces is generating a component of E in direction x 3 while the others remain zero. The piezoelectric coefficient couples the electric with the mechanic field generating displacements. Since the x 3 direction is the polarized, the induced displacements in all directions are proportional to coefficients e e,31 , e e,32 and e e,33 respectively. As the boundary conditions do not prevent these movements, no stresses appear and the equality Cε = e e E holds.

The appearance of strains modifies the electric displacement D field. However, no change is expected in E since the electric potential is not affected by the coupling in this particular example.

In Figure 2.3, a summary of all relevant distributions calculated by FEAP is shown; these distributions are linear as expected. All displacements are proportional to the piezoelectric coefficients mentioned before (x 2 direction distribution not displayed but equal to the one of x 1 ). Note that u 3 is negative while u 1 and u 2 are positive, due to the numerical values assigned to e e,31 and e e,32 . The also linear V cause a constant distribution of E 3 , not shown in the figure for simplicity.

The complementary case, for which the piezoelectric material is used as generator is also run for another validation. The problem is set with the same boundary conditions as in Figure 2.2, but instead of V = 10 V, a vertical uniform displacement of a hundreth of the x 3 length is prescribed at the top.

-9.84E-06 -7.87E-06 -5.90E-06 -3.94E-06 -1.97E-06 -1.18E-05 0.00E+00 Displacement in x direction (m) 1.67E-06 3.33E-06 5.00E-06 6.67E-06 8.33E-06 0.00E+00 1.00E-05 Displacement in x direction (m) 2.92E+03 5.84E+03 8.76E+03 1.17E+04 1.46E+04 0.00E+00 1.75E+04 Electric potential (V) 1 3
Figure 2.4: Displacement in directions x 1 and x 3 and electric potential distributions for the generator validation piezoelectric case.

The displacement u 3 is linear due to the applied prescription (see Figure 2.4), while u 1 appears because of the crossed coefficients C 13 and C 23 being negative due to the Poisson ratio. The piezoelectric coupling generates a linear V distribution and therefore a constant E 3 . Analogously to the previous example, the boundary conditions enforce a zero D, therefore e e ε = -E.

Another example, not so straightforward, is run under the boundary conditions prescribed of Figure 2.5. This time a quarter of the box is being modeled, due to the non-symmetric mechanic boundary conditions: at the bottom all the degrees of freedom, including the electric potential, are set to zero. The other boundary conditions remain equal, that is, with symmetry in planes x 1 = x 2 = 0 and voltage of 10 V at the top.

In this case, the plane x 3 = 0 is clamped, concentrating all stresses there as can be seen in Figure 2.6. Both u 1 and u 2 displacements are skewed due to the boundary conditions. To fully appreciate the movement, the deformed configuration is also given.

The electric potential may seem linear, but it can be appreciated in the figure that near the bottom the isolines are closer among themselves than at the top; near the right free edge the isolines are also closer and not straight. Hence, the E 3 component is not constant and for E 1 some concentrations appear near the lower frontal edge (the distribution for E 2 is the same). The value of the latter is an order of magnitude lower than the former.

V = 10 V x x x 2 3 1 2 3 3 x 2 x 1 Figure 2
.5: Sketch of the boundary conditions for the third piezoelectric example. Two planes of symmetry considered in x 1 = x 2 = 0; electric potential prescriptions at the top and bottom planes; the latter face clamped. Measures in mm.

Piezo-electro-magnetism

First, a simple piezomagnetic case with the same configurations as before is studied; the only non-zero tensors are C, e h and µ. Instead of electric potential for the boundary conditions, magnetic potential is applied. The obtained results give similar values to the previous ones, but with different numbers due to the difference in e h and µ with respect to e e and . This coincidence is due to the similarities between the electric and magnetic field.

A more complex piezo-electro-magnetic material is simulated ; therefore all the previous piezoelectricity and piezomagnetism tensors and also υ are active now. With this hypothetical material, the same example as in Figure 2.5 but with ϕ = 10 C/m at the top and ϕ = 0 C/m at the bottom is run.

In Figure 2.7, a mosaic with the most important magnitudes calculated by FEAP can be seen. The displacements come from the superposition of electric and magnetic fields; for instance, the electric field generates positive u 1 while the magnetic field induces negative ones. Since e h,31 is higher than e e,31 , the sum of both displacements is a negative one. Otherwise, both fields generate negative displacements for u 3 , so that all contribute are in the same direction. The deformed configuration amplified by 2 × 10 4 times, and it is given in the central figure for better understanding of the displacements.

As for the electric and magnetic fields, their distributions are basically the same: for E 1 and H 1 a concentration appears near the bottom front edge due to the boundary conditions, and for E 3 and H 3 the concentration is near the bottom plane; in the free edge a homogeneous distribution can be seen. The higher value of µ with respect to is what causes the almost insignificant variation of H 3 and a negligible H 1 . For the electric field, E 1 is only one order of magnitude lower than E 3 . Since H is the predominant field, the u 3 distribution is more regular than the one displayed in Figure 2.6. 1.67E+00 3.33E+00 5.00E+00 6.67E+00 8.33E+00 0.00E+00

1.00E+01 Electric potential (V) -3.58E+01 4.38E+01 1.23E+02 2.03E+02 2.83E+02 -1.16E+02 3.62E+02 Electric field in x direction (V/m) -6.50E+03 -5.94E+03 -5.39E+03 -4.83E+03 -4.27E+03 -7.06E+03 -3.72E+03 Electric field in x direction (V/m) -3.69E+05 -2.85E+05 -2.01E+05 -1.18E+05 -3.43E+04 -4.52E+05 4.92E+04 Stress in x direction (Pa) -3.42E+05 -2.70E+05 -1.97E+05 -1.25E+05 -5.33E+04 -4.14E+05 1.88E+04
Stress in x direction (Pa)

1 3 1 3
Figure 2.6: Displacement in directions x 1 and x 3 and deformed configuration amplified 10 5 .

Thermoelasticity

For the thermoelastic behavior the only tensors or coefficients not null from the constitutive (2.41) are β, C and ρ m c p . A first benchmark example is calculated under the same geometry and mechanical boundary conditions as in Figure 2.2. Temperature boundary conditions are required instead of the electric ones: the value for T at the top and the bottom is T 0 + 20 K and the other faces have adiabatic boundary conditions. This example is run in FEAP for static case obtaining the results in Figure 2.8. Since T is both a degree of freedom and a dual variable, the regular benchmark result could only be obtained with constant temperature. In this case, a linear distribution for all displacements is obtained and the distributions are the same as expected from analytical results. All displacements are positive since T is greater than the reference temperature. Despite the fact that that the expansion coefficient is isotropic α T = 6.16 × 10 -6 1/ • C, β is not due to the elasticity tensor C.

An opposite benchmark result could be obtained imposing a constant and homogeneous value for ε = 0.01 and keeping only the boundary for T = T 0 at the bottom, imposing the adiabatic condition also at the top side so that the temperature can change. The imposed displacement is introduced at the top first with an increasing ramp until t = 1 s and then returning to zero with the same pace until t = 2 s; the movement at the bottom is restrained. In Figure 2.9 left the transient distributions of T on the top surface are plot for two different heights measured along the coordinate x 3 . When ε is positive T significantly decreases up to 0.25 s and after remains constant until the ramp sign is changed, stabilizing the temperature in a steady-state that can be calculated analytically. When ε is negative a smooth transition to a T symmetric increase is observed. The effect is much less pronounced when the height is small (l = 1 mm), since it is analytically demonstrated that it depends on l 2 . This analytic solution is obtained solving a more easy version of the heat equation with constant ε, assuming that Poisson modulus is null and unidimensionality. The temperature expression derived is:

-5.12E-08 -4.09E-08 -3.07E-08 -2.05E-08 -1.02E-08 -6.14E-08 0.00E+00 Displacement in x direction (m) 1.67E+00 3.33E+00 5.00E+00 6.67E+00 8.33E+00 0.00E+00 1.00E+01 Electric potential (V) 1.67E+00 3.33E+00 5.00E+00 6.67E+00 8.33E+00 0.00E+00 1.00E+01 Magnetic potential (C/m) -4.57E-09 -3.66E-09 -2.74E-09 -1.83E-09 -9.14E-10 -5.48E-09 0.00E+00 Displacement in x direction (m) -5.14E+02 -3.72E+02 -2.29E+02 -8.63E+01 5.64E+01 -6.57E+02 1.99E+02 Electric field in x direction (V/m) -6.29E+03 -5.31E+03 -4.33E+03 -3.36E+03 -2.38E+03 -7.27E+03 -1.40E+03 Electric field in x direction (V/m) 7.20E-02 2.11E+00 4.14E+00 6.18E+00 8.22E+00 -1.96E+00 1.03E+01 7.20E-02 2.11E+00 4.14E+00 6.18E+00 8.22E+00 -1.96E+00 1.03E+01 -5.03E+03 -5.02E+03 -5.00E+03 -4.99E+03 -4.98E+03 -5.04E+03 -4.97E+03 Magnetic field in x direction (A/m)
T (x 3 ) = T 0 β 3 ε3 2κ x 3 (x 3 -2l) + T 0 (2.68)
In the right figure, results for this solution (solid line) and for the numerical FEAP (circles) along the height are shown; the bottom two lines are at an instant right before 1 s, the top ones right before 2 s, all of them when T is stabilized. It is clear that this temperature varies quadratically with the height and that again the effect of the thermoelasticity coupling is much more evident for the large height. In this example we consider a transient problem for which different coupling effects can be presented, thus providing the sort of qualitative benchmark result. Namely, we compute a thermo-electro-magneto-mechanical coupling in 3D domain in the form of a bar, with x 3 dimension significantly larger than the other two (1 × 1 × 100 mm); this dimension is chosen for preferential polarization and magnetization. The boundary conditions are chosen taking into account the propagation to be studied along the bar: a displacement varying in time is imposed at the bar left end, while the right end is free. We also impose adiabatic boundaries for the thermal and isolated contours for the electromagnetic fields (see Figure 2.10). The choice of Newmark scheme parameters γ = 0.5 and β = 0.25 is in agreement with the trapezoidal second order scheme.

The bar is set in motion by a sinusoidal displacement pulse imposed on the left end of maximum amplitude u 3,max = 10 mm; this pulse propagates through the bar, with a period much lower than that of the bar. This displacement produces the corresponding local change in strain and in its time derivative, which further generates the couplings in (2.41). In Figure 2.11 we present the computed response for the displacement in the longitudinal direction of the bar, for the electric potential, for the magnetic potential and for the temperature in the middle x 3 = l/2, through the transient sequence. The time scale used for all of them is the step size ∆t = 5 × 10 -8 s.

Regarding the evolution of the mechanical field (first figure), this problem can be examined as a wave propagation: when the wave arrives to the right end, the reflected wave keeps the same sign due to the free end of the bar. However, when the reflected wave arrives to the left, the sign changes due to the built-in boundary. Since there is no dissipation implemented in the formulation, the wave conserves its maximum value and its period along the simulation. This displacement wave creates perturbations in the other dual variables through the tensors e e , e h and β, creating variations in the other degrees of freedom.

The temperature is analyzed first; this degree of freedom is proportional to the strain rate and also depends on the direction of the wave through (2.44). Thus, when a wave of positive displacement comes from the left to the right (inducing a positive heat gradient), T experiences first a decrease and second an increase, driven by first a negative strain rate and second a positive one in the same wave. After the wave bounces off the free end, the heat gradient is now negative, then T first increases and second decreases. The contrary happens when the displacement wave is negative.

-1.68E+04 7.40E+03 3.16E+04 5.58E+04 8.00E+04 -4.10E+04 1.04E+05 Magnetic Potential (C/m) -1.22E+06 -1.14E+06 -1.06E+06 -9.80E+05 -9.01E+05 -1.30E+06 -8.21E+05 Electric potential (V)
Figure 2.12: Contour plot for electric and magnetic potential at t = 1.7 × 10 -5 s.

As for the electric and magnetic fields, their response cannot be interpreted directly since the absence of essential boundary conditions prevents the uniqueness of the solution. In Figure 2.12, the waves in a contour plot from FEAP for both V and ϕ can be appreciated. The images are taken at the instant when the maximum value of u 3 is at the bar middle. The shape of both waves is equal to that of the elastic field, but the zones outside this wave influence have values different than zero.

Elastic pulse under complete coupling for conductor formulation

A similar example is tested for the conductor formulation. The chosen material properties are shown in Table 2.2. This time, the choice of Newmark scheme parameters γ = 1.5 and β = 1 results with the highest numerical dissipation (see [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF]. This displacement produces the corresponding local change in the time derivative of the strain, which further generates the structural heating defined in (2.44) and leads to a temperature change. These temperature changes also lead to change in the electric field, which is affected through the constitutive behavior of the material (2.28) assuming J = 0 due to the isolated boundary conditions.

Property Units Value Young Modulus E GPa 200 Mass density ρ m kg/m 3 7.8 × 10 3 Specific heat c p J/kg K 1.2 × 10 6 Seebeck coefficient α V/K 2 × 10 -4 Thermal expansion coef. α T 1/K 1.5 × 10 -5 Thermal conductivity κ W/K m 0.15 Electric conductivity γ A/V m 10 6
Table 2.2: Properties of the material used in validation example

Regarding the evolution of the mechanical field, the maximum displacement is gradually reduced due to the numerical dissipation by the Newmark parameters, which produces a similar effect than the dissipation by conduction for non-homogenize temperature distribution along the bar. In addition, the wave length is increased in the course of this wave propagation.

The evolution of other two fields are directly affected by the evolution of displacement field. The temperature in particular is proportional to the time derivative of the strain or strain rate, which changes the sign within the displacement pulse. Moreover, when the wave passes by the middle point, the strain will be positive or negative, which depends upon the direction of the wave and the sign of the displacement. As it can be seen in Figure 2.13, a wave passing from left to right in the ascending part generates at first a positive strain rate, but once passed the maximum displacement, the same wave produces negative strain rate. The voltage evolution is directly affected by the evolution of the temperature, but with the sign changed as stated in (2.28).

Macroscopic inelastic behavior

In this chapter, a macroscopic isotropic hardening model for plasticity is introduced for the formulations developed in the previous chapter. The plastic evolution equations are derived consistently, as well as the tangent thermoelectro-magneto-elastoplastic modulus through dual variable computations. The latter plays a crucial role in ensuring fast convergence properties of the finite element computations with the proposed coupled plasticity model. The thermal, electric, and magnetic fields play an essential role in defining the criteria for the activation of the plastic regime. Several numerical simulations are presented in order to illustrate the proposed model and formulation capabilities for providing an enhanced formulation of a relevant practical application.

Plasticity under thermo-electro-magneto-mechanics coupling

We seek to develop here the general form of associative plasticity for a coupled problem as the ones presented in the previous chapter. Dissipation has already been defined in (2.33). We now turn to the case when plasticity is activated. In the framework of linear kinematics defined in (2.1), it is generally accepted to use the additive decomposition of the state variables into elastic and plastic parts:

ε = ε e + ε p ; s = s e + s p ; D = D e + D p ; B = B e + B p (3.1)
The internal energy potential is now defined in terms of the elastic part of the state variables ε e , s e , D e , B e . For describing eventual hardening phenomena, which makes the model more predictive than perfect plasticity [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], we also need to add a new state variable that can monitor the progress of plasticity threshold, here chosen in terms of isotropic hardening variable ζ. As in the previous chapter, we introduce the free energy potential ψ by means of the Legendre transformation (e.g., [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], exchanging the roles between the state variables and their duals, s e , D e , B e versus T, E, H:

ψ (ε e , ζ, T, E, H) = e (ε e , ζ, s e , D e , B e ) -T s e -E • D e -H • B e (3.2)
The time derivative of the last expression leads to:

ė = ∂ψ ∂ε e • εe + ∂ψ ∂ζ -q ζ + ∂ψ ∂T Ṫ + ∂ψ ∂E • Ė + ∂ψ ∂H • Ḣ + Ṫ s e + T ṡe + E • Ḋe + Ė • D e + H • Ḃe + Ḣ • B e (3.3)
where q is stress-like variable thermodynamically conjugate to ζ, which controls the evolution of plasticity threshold. With this result on hand, the dissipation in (2.33) can be expressed as:

0 ≤ D = σ - ∂ψ ∂ε e • εe -s e + ∂ψ ∂T Ṫ -D e + ∂ψ ∂E • Ė -B e + ∂ψ ∂H • Ḣ + q ζ + σ • εp + T ṡp + E • Ḋp + H • Ḃp + J • E (3.4)
By considering the elastic process, where (possibly non-zero values of) internal variables remain frozen enforcing that εp = 0, ζ = 0, ṡp = 0, Ḋp = 0, Ḃp = 0, and putting aside the Joule's dissipation, we can conclude that the local plastic dissipation remains equal to zero. The dissipation inequality will become an equality providing the set of constitutive equations to be defined in agreement with the chosen free energy potential. Here, we choose a quadratic form of that potential, which can be written for conductors as follows:

ψ (ε e , ζ, T, E, H) = ψ m + ψ p + ψ t + ψ e + ψ h + ψ mt + ψ me + ψ mh ; ψ m = 1 2 ε e • Cε e ; ψ h = - 1 2 H • µH ; ψ p = 1 2 ζKζ ; ψ mt = -β • (T -T 0 ) ε e ; ψ t = ρ m c p (T -T 0 ) -T ln T T 0 ; ψ me = (E ⊗ E) • ε e - 1 2 (E • E) tr (ε e )
;

ψ e = - 1 2 E • E ; ψ mh = µ (H ⊗ H) • ε e - 1 2 (H • H) tr (ε e )
(3.5) For insulators:

ψ (ε e , ζ, T, E, H) = ψ m + ψ p + ψ t + ψ e + ψ h + ψ mt + ψ me + ψ mh + ψ te + ψ th + ψ eh ; ψ m = 1 2 ε e • Cε e ; ψ me = -ε e • e e E ; ψ p = 1 2 ζKζ ; ψ mh = -ε e • e h H ; ψ t = ρ m c p (T -T 0 ) -T ln T T 0 ; ψ te = -(T -T 0 ) π e • E ; ψ e = - 1 2 E • E ; ψ th = -(T -T 0 ) π h • H ; ψ h = - 1 2 H • µH ; ψ eh = -E • υH ψ mt = -(T -T 0 ) β • ε e (3.6)
The constitutive equations can then be obtained from derivatives of such a potential, and written in tensor notation:

σ • • = ∂ψ ∂ε e = C (ε -ε p ) -β (T -T 0 ) + E ⊗ E + µH ⊗ H - I 2 ( E • E + µ H • H) ; q • • = - ∂ψ ∂ζ = -Kζ ; D e • • = - ∂ψ ∂E = {[1 + tr (ε e )] E -ε e E} ; B e • • = - ∂ψ ∂H = µ {[1 + tr (ε e )] H -ε e H} ; s e • • = - ∂ψ ∂T = ρ m c p ln T T 0 + β • (ε -ε p ) (3.7)
For insulators:

σ • • = ∂ψ ∂ε e = C (ε -ε p ) -β (T -T 0 ) -e e E -e h H q • • = - ∂ψ ∂ζ = -Kζ ; s e • • = - ∂ψ ∂T = β • (ε -ε p ) + ρ m c v ln T T 0 + π e • E + π h • H D e • • = - ∂ψ ∂E = e e (ε -ε p ) + π e (T -T 0 ) + E + υH B e • • = - ∂ψ ∂H = e h (ε -ε p ) + π h (T -T 0 ) + υE + µH (3.8)
By assuming that such constitutive equations also remain valid in an inelastic process, the dissipation can now be defined as the sum of plastic and Joule terms. Furthermore, the plastic dissipation can be split into mechanic, thermal and electromagnetic terms defined explicitly as:

D = D p + D J ; D p = q ζ + σ • εp D p m + T ṡp D p t + E • Ḋp + H • Ḃp D p eh (3.9)
With these results on hand, we can rewrite (2.33) in terms of generalizing the heat equation for inelastic case which can be written as:

T ṡ = r -∇ • q + D ⇒ T ṡe = r -∇ • q + D p m + D p eh + D J (3.10)
where the results (3.1) and (3.9) are used. Furthermore, by using the constitutive equation for s e in (3.7) 5 or (3.8) 3 , depending on the formulation, we can write explicitly:

T ṡe • • = ρ m c p Ṫ + T β • (ε -εp ) ; For conductors T ṡe • • = ρ m c p Ṫ + T β • (ε -εp ) + π e • Ė + π h • Ḣ ; For insulators (3.11)
By introducing (3.11) into (3.10) and using (3.7), the generalized heat equation for conductors can finally be restated as:

ρ m c p Ṫ = r -∇ • q -T β • (ε -εp ) + D J + D p m + D p eh (3.12)
Analogously, taking into account that J = 0, for insulators:

ρ m c p Ṫ = r -∇ • q -T β • (ε -εp ) + π e • Ė + π h • Ḣ + D p m + D p eh (3.13)
In summary, the strong form of the equations to be solved for conductors regroups the results written in (2.10) 2 , (2.6) 4 , (2.11) and (3.12), here restated in tensor notation:

ρ m ü = ∇ • σ + b ∇ • J = 0 ∇ • B = 0 ρ m c p Ṫ = r -∇ • q -T β • (ε -εp ) + D J + D p m + D p eh (3.14)
Similarly, regrouping the results for inductors in (2.8) 2 , (2.6) 4 , (2.20) and (3.13), here restated in tensor notation:

ρ m ü = ∇ • σ + b ∇ • D = ρ f q ∇ • B = 0 ρ m c p Ṫ = r -∇ • q -T β • (ε -εp ) + π e • Ė + π h • Ḣ + D p m + D p eh (3.15)

Evolution equations for internal variables

We here use the generalized form of von Mises criterion that builds upon the proposal in [START_REF] Landis | Non-linear constitutive modeling of ferroelectrics[END_REF] for perfect plasticity, which allows to determine if the material is in elastic or plastic regime based upon its stress state and electric field values. The first generalization of such criterion concerns the possibility to include the thermal field, leading to:

φ(σ, q, T, E) • • = 3 2 dev(σ) 2 σ 2 y (T ) + E 2 E 2 0 + H 2 H 2 0 -1 = 0 (3.16)
where

dev(A) = A -1 3 [tr(A)]
is the deviatoric part of the tensor A, A = tr A T A is the Euclidean norm of A, σ y (T ) is the temperature dependent yield stress, and E 0 and H 0 reference values for electric and magnetic fields.

The same criterion can be recast in a more standard format (e.g., [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]) by making the yield stress dependent upon thermal and electric field variations, which also allows to account for potential hardening effect in terms of isotropic hardening. Such a generalized von Mises criterion can be written as:

φ (σ, q, T, E, H) = dev(σ) - 2 3 [σ y (T, E, H) -q (T, E, H)] = 0 ; q = -K (T, E, H) ζ (3.17)
where the yield stress and the isotropic hardening modulus are now assumed dependent upon thermal and electric fields, according to the following dependency:

σ y (T, E, H) = σ y 0 [1 -ω t (T -T 0 )] [1 -ω e ( E -E 0 )] [1 -ω h ( H -H 0 )] K (T, E, H) = K 0 [1 -ω t (T -T 0 )] [1 -ω e ( E -E 0 )] [1 -ω h ( H -H 0 )] (3.18)
where ω t , ω e , and ω h are some coefficients to define the influence of the temperature and the electric fields respectively in σ y and K.

By using the principle of maximum dissipation [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], we will pick among all plastically admissible values of dual variables (which satisfy φ(σ, q, T, E, H) = 0) the one which will maximize the plastic dissipation. We can further recast such a constrained optimization problem in terms of min-max problem of unconstrained minimization where the plastic admissibility constraint is enforced by using the method of Lagrange multipliers leading to plastic Lagrangian:

max φ(σ,q,T,E,H)=0 D p (σ, q, T, E, H) ⇔ min ∀(σ,q,T,E,H) max γ>0 L p σ, q, T, E, H, γ (3.19)
where the Lagrange multiplier is equal to the plastic multiplier γ. The plastic Lagrangian will incorporate the yield criterion constraint φ = 0, which can be written as:

L p σ, q, T, E, H, γ • • = -D p (σ, q, T, E, H) + γp φ (σ, q, T, E, H) (3.20)
The solution to such constrained minimization problem can be obtained by using the Kuhn-Tucker optimality conditions, which allows to write:

0 • • = ∂L p ∂σ = γ ∂φ ∂σ -εp ; 0 • • = ∂L p ∂q = γ ∂φ ∂q -ζ ; 0 • • = ∂L p ∂T = γ ∂φ ∂T -ṡp ; 0 • • = ∂L p ∂E = γ ∂φ ∂E - Ḋp ; 0 • • = ∂L p ∂H = γ ∂φ ∂H - Ḃp ; γ ≥ 0 ; φ ≤ 0 ; γ φ = 0 (3.21)
Using the first two previous equations, D p m can further be simplified:

D p m = γ φ=0   σ • ν dev(σ) - 2 3 (σ y -q)   + γ 2 3 σ y = γ 2 3 σ y (3.22)
where ν • • = ∂φ/∂σ = dev(σ)/ dev(σ) .

Finite element implementation

When plasticity is introduced, the computation is divided into two phases: local and global. In the local phase, the current values of the plastic variables are computed inside every element, whereas the values of the degrees of freedom remain frozen. No assembly operation is performed in this phase because the plastic variables are not continuous among elements. In the global phase, the plastic variables calculated in the local phase are used to compute the stiffness matrix and the force vector for each element. Then, they are assembled into the finite element mesh via the connections of nodes. Finally, the new values of the degrees of freedom are solved for the whole geometry. This process is repeated until the variables are converged.

Local iterative solution for plasticity with isotropic hardening

This local computation has to be performed for every Gauss point in agreement with the corresponding result a (i) n+1 obtained in the global phase. For proposed generalization of von Mises plasticity, we first define the yield function for coupled problem of this kind:

φ n+1 a (i) n+1 = dev(σ n+1 ) - 2 3 (σ y (T n+1 , E n+1 , H n+1 ) -q n+1 ) (3.23)
The elastic trial step is computed first, with the zero value of plastic multiplier γ p,tr n+1 = 0. This result with:

ε p,tr n+1 = ε p n ; ζ tr n+1 = ζ n ; s p,tr n+1 = s p n ; D p,tr n+1 = D p n ; B p,tr n+1 = B p n (3.24)
These results are exploited in computation of the corresponding trial values of stress and hardening stress-like variable for conductors:

σ tr n+1 = C ε n+1 -ε p n -β (T n+1 -T 0 ) + E ⊗ n+1 - 1 2 E n+1 • E n+1 ; q tr n+1 = q n = -K ζ n (3.25)
and for insulators:

σ tr n+1 = C ε n+1 -ε p n -β (T n+1 -T 0 ) -e e E n+1 -e h H n+1 ; q tr n+1 = q n = -K ζ n (3.26)
To ensure that the elastic step is acceptable, the trial value of yield function must remain negative or zero

0 ≥ φ tr n+1 = dev σ tr n+1 - 2 3 σ y (T n+1 , E n+1 ) -q tr n+1 (3.27)
If such a condition is not verified, the step is plastic. The latter implies that all the internal variables must be recomputed by integrating the evolution equations in (3.21):

ε p n+1 = ε p n + γ p n+1 ν n+1 ; ζ n+1 = ζ n + 2 3 γ p n+1 ; s p n+1 = s p n + γ p n+1 2 3 ω t (σ y 0 + K 0 ζ n+1 ) ; D p n+1 = D p n + γ p n+1 E E 2 3 ω e (σ y 0 + K 0 ζ n+1 ) ; B p n+1 = B p n + γ p n+1 H H 2 3 ω h (σ y 0 + K 0 ζ n+1 ) (3.28)
The plastically admissible value of stresses can be computed as:

σ n+1 = σ tr n+1 -γ p n+1 C ν n+1 q n+1 = -K ζ n+1 (3.29)
We note that the final value of plastic multiplier γ p n+1 can be obtained from φ n+1 = 0, by taking into account that:

ν n+1 = ν tr n+1 ; dev (σ n+1 ) = dev σ tr n+1 -2µ γ p n+1 (3.30)
where µ is the shear modulus. We finally obtain:

γ p n+1 = φ tr n+1 2µ + 2K 3 (3.31)

Consistent tangent thermo-electro-elastoplastic modulus

The computation of the consistent tangent thermo-electro-elastoplastic modulus is carried by generalizing the scheme first proposed for classical plasticity in [START_REF] Ibrahimbegovic | Classical plasticity and viscoplasticity models reformulated: theoretical basis and numerical implementation[END_REF]. At converged values of local computation, we take a new iterative sweep brought by an increment of total strain. The corresponding linearized form of local problem is obtained in terms of dual variables:

L (k) n+1         ∆σ n+1 ∆T n+1 ∆E n+1 ∆H n+1 ∆γ p n+1 ∆q n+1         =                    -ε e n+1 + ε e n + ∆ε (i) n+1 -γ p n+1 ∂φ ∂σ n+1
-s e n+1 + s e n + ∆s

(i) n+1 -γ p n+1 ∂φ ∂T n+1 -D e n+1 + D e n + ∆D (i) n+1 -γ p n+1 ∂φ ∂E n+1 -B e n+1 + B e n + ∆B (i) n+1 -γ p n+1 ∂φ ∂H n+1 φ n+1 ζ n+1 -ζ n -γ p n+1 ∂φ ∂q n+1                    (3.32) The tangent matrix L (k)
n+1 is defined for the most general case as:

L • • = L aa (13×13) L ab (13×2) L T ab (2×13) L bb (2×2) (3.33)
where the submatrices are:

L aa • • =              ∂ 2 ψ ∂ε e ∂ε e -1 + γ ∂ 2 φ ∂σ ∂σ γ ∂ 2 φ ∂σ ∂T γ ∂ 2 φ ∂σ ∂E γ ∂ 2 φ ∂T ∂σ γ ∂ 2 φ ∂T 2 - ∂ 2 ψ ∂T 2 γ ∂ 2 φ ∂T ∂E γ ∂ 2 φ ∂E ∂σ γ ∂ 2 φ ∂E ∂T γ ∂ 2 φ ∂E ∂E - ∂ 2 ψ ∂E ∂E γ ∂ 2 φ ∂H ∂σ γ ∂ 2 φ ∂H ∂T γ ∂ 2 φ ∂H ∂E γ ∂ 2 φ ∂H ∂H - ∂ 2 ψ ∂H ∂H              ; L ab • • =             ∂φ ∂σ 0 ∂φ ∂T 0 ∂φ ∂E 0 ∂φ ∂H 0             ; L bb • • =     0 ∂φ ∂q ∂φ ∂q γ ∂ 2 φ ∂q 2 - ∂ 2 ψ ∂q 2 -1    
(3.34) The subsequent iterative values can then be obtained with simple additive updates:

σ (k+1) n+1 = σ (k) n + ∆σ n+1 T (k+1) n+1 = T (k) n + ∆T n+1 E (k+1) n+1 = E (k) n + ∆E n+1 H (k+1) n+1 = H (k) n + ∆H n+1 γ p,(k+1) n+1 = γ p,(k) n + ∆γ p n+1 q (k+1) n+1 = q (k) n + ∆q n+1 (3.35)
With the converged values of internal variables in the local computation phase, we can write the result in (3.32) as:

L n+1         ∆σ n+1 ∆T n+1 ∆E n+1 ∆H n+1 γ p n+1 q n+1 -q n         =          ∆ε (i+1) n+1 ∆s (i+1) n+1 ∆D (i) n+1 ∆B (i) n+1 0 0          (3.36)
Further, by using the static condensation method (e.g., [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF] we can reduce the size of this system:

L aa -L ab L -1 bb L ba     ∆σ ∆T ∆E ∆H     n+1 =     ∆ε ∆s ∆D ∆B     (i+1) n+1 (3.37)
Due to the zero column and row of L ab and L ba respectively, the only relevant term of L -1 bb for the multiplication in (3.37) is the first one:

L • • = -L -1 bb (1, 1) = b a 2 (3.38)
where:

a • • = L bb (1, 2) = L bb (2, 1) = ∂φ ∂q ; b • • = L bb (2, 2) = γ ∂ 2 φ ∂q 2 - ∂ 2 ψ ∂q 2 -1 (3.39)
so that this multiplication is simplified to:

-L ab L -1 bb L ba = L ξ ξ T (3.40)
where ξ is a vector containing the first column of L ab . By substituting (3.40) into (3.37), after some further simplifications, we obtain:

L aa + L ξ ξ T Ĉ     ∆σ ∆T ∆E ∆H     n+1 =     ∆ε ∆s ∆D ∆B     (i+1) n+1 (3.41)
featuring Ĉ as a consistent tangent thermo-electro-elastoplastic modulus for the given variables. By exchanging the roles between variables σ and ε, we can obtain the final form of tangent modulus for the global phase C; namely, with the split of Ĉ:

Ĉ = Ĉaa (6×6) Ĉab (6×7) Ĉba (7×6) Ĉbb (7×7) (3.42)
we can easily construct C which can be written as:

C = Ĉ-1 aa (6×6) -Ĉ-1 aa Ĉab (6×7) Ĉba Ĉ-1 aa (7×6) Ĉbb -Ĉba Ĉ-1 aa Ĉab (7×7) ⇒ C     ∆ε ∆T ∆E ∆H     n+1 =     ∆σ ∆s ∆D ∆B     n+1 (3.43)

Global step

This step is similar to the one in the previous chapter. The weak form of the equations is set out first. Then, the approximations for the degrees of freedom and the derivatives are introduced into the formulation. Finally, the virtual degrees of freedom are isolated, obtaining accordingly the residuals. The modifications with respect to the previous section include a new residual expression for conductors:

R U,(i) a,n+1 = -ˆΩ B sy a T σ (i) n+1 + N a ρ m N b äU,(i) b,n+1 -b n+1 dΩ + ˆΓ N a t C n+1 dΓ ; R V,(i) a,n+1 = ˆΩ B T a J (i) n+1 dΩ -ˆΓ N a j n+1 dΓ ; R T,(i) a,n+1 = ˆΩ B T a q (i) n+1 -N a ρ m c p N b ȧT,(i) b,n+1 -r n+1 +N b a T,(i) b,n+1 β T B sy b ȧU,(i) b,n+1 - εp,(i) n+1 + J (i) n+1 + Ḋp n+1 T B b a V,(i) b,n+1 -γp n+1 2 3 σ y dΩ -ˆΓ N a q n+1 dΓ (3.44)
and for insulators:

R U,(i) a,n+1 = -ˆΩ B sy a T σ (i) n+1 + N a ρ m N b äU,(i) b,n+1 -b (i) n+1 dΩ + ˆΓ N a t c n+1 dΓ ; R V,(i) a,n+1 = ˆΩ B T a D (i) n+1 + N a ρ f q,n+1 dΩ -ˆΓ N a D n+1 dΓ ; R ϕ,(i) a,n+1 = ˆΩ B T a B (i) n+1 dΩ -ˆΓ N a B n+1 dΓ ; R T,(i) a,n+1 = ˆΩ B T a q (i) n+1 -N a ρ m c p N b ȧT,(i) b,n+1 -r n+1 +N b a T,(i) b,n+1 β T B sy b ȧU,(i) b,n+1 - εp,(i) n+1 -π T e B b ȧV,(i) b,n+1 -π T h B b ȧϕ,(i) b,n+1 -γp n+1 2 3 σ y - Ḋp,T n+1 B b a V,(i) b,n+1 - Ḃp,T n+1 B b a ϕ,(i) b,n+1 dΩ -ˆΓ N a q n+1 dΓ (3.45)
The stiffness terms are obtained by deriving the residuals with respect to the nodal variables as showed in (2.63). If plasticity is activated for any time step, the consistent elastoplastic modulus should be used instead of the regular ones.

Numerical example: pulsed Peltier cells under plasticity constraint

This formulation can be used in a vast number of practical applications, such as Peltier cells, where we need to account for plasticity phenomena in metallic materials. In this numerical example, we simulate the behavior of half a thermocouple of a pulsed Peltier cell, representing a cooler based on the Peltier effect under the action of electric flux. This cell consists of four different materials: thermoelectric material (here chosen Bi 2 Te 3 ), copper, alumina, and tin-lead solder.

A complete simulation for elastic regime can be found in (Pérez-Aparicio et al., 2016). The material properties for this problem are presented in Table 3.1. Moreover, some of these Finally, the alumina is considered to remain elastic material (setting very high value for yield stress) because the plasticity begins well after the ultimate stress is reached. The same values for the Newmark parameters are chosen as in the last example. The finite element mesh constructed with the hexaedral 8-node finite elements used for computing the solution to this problem is shown in Figure 3.1. Different materials are presented in different colors. The boundary conditions are shown in Figure 3.1 right. We fix the displacements in direction x 1 and x 3 at the respective end faces, along with the symmetry condition in the plane x 2 = 0. The temperature is prescribed at the top face T t = 50 • C (called hot side) along with a zero heat flux imposed at the bottom (cold side) where the temperature T b is left free. The electric flux is introduced at the bottom left end of the copper while ground voltage is assumed at the top right end copper. The applied loading is presented in Figure 3.2 top. First, an electric flux J ss = I op /A is introduced, where I op is an optimal electric current that maximizes the temperature difference T t -T b in steady-state (more details in Pérez-Aparicio et al., 2016)) with A as the transversal area of the copper. This electric flux value is kept fixed until evolution for all the degrees of freedom have been stabilized, and a steady-state has been reached. After this moment, the flux is incremented P times, where P = J/J ss is the pulse gain, which implies automatically that T b decreases until a minimum is reached, but with an overheating penalty reached right afterwards. The steady-state temperature is reached eventually after reimposing P = 1 if no other pulse is introduced. Figure 3.2 provides the representation of the maximum von Mises stress values calculated in Gauss points for each material in both elastic and plastic regime. In the steady-state, the maximum stresses are similar since the yield stress has not been reached yet. However, with the increase of electric flux, the stress increments seen in (Pérez-Aparicio et al., 2016) due to accumulation of Joule heating, are no longer the same.

Property

When the pulse is introduced reaching the yield stress σ y of Bi 2 Te 3 , the subsequent stress-rate increase is reduced, due to the corresponding increase in plastic strains. These changes result in a noted decrease of the maximum stresses computed with activated plasticity constraint compared to the elastic case, specially in the copper and in the thermoelectric material. Once the steady-state gain is restated, the asymptotic stresses are reduced as a result of the plastic deformation in the copper and the solder. Moreover, due to the stresses redistribution, enforced by plasticity constraint, the stresses are increased in the Bi 2 Te 3 . Thus, a new study of the introduction of another pulse after the application of the first is needed. In Figure 3.3 the nodal von Mises projection distribution is shown at the end of the pulse t = 130 s, when the maximum stresses take place, along with a detailed zoom of the cold side, the most stressed part of the structure. The alumina has been omitted in order to show more clearly the differences with respect to the plastic materials.

In both cases the maximum is located at the bottom since this side has a bigger difference with respect to the reference temperature than the hot side. The maximum values are rather similar in both regimes. This is due to the influence of alumina that remains elastic, which stabilizes the corresponding stress distribution inside with no marked changes.

The first noticeable difference is that the stress distribution of the lower values is quite different. This can be appreciated specially for the copper, where the part on the bottom left has quite reduced values of stress, and in the Bi 2 Te 3 where they have increased in average in a homogenization process. This further explains what was already observed in Figure 3.2, as in general in the copper the stresses in steady-state will be reduced and the opposite for the thermoelectric material.

It is interesting to note (as shown in Figure 3.4 bottom) that the temperature evolution is not very affected by the presence of plasticity unlike the stress distribution. The latter is due to rather small contribution of the plastic dissipation towards total power sources, especially 4

Discrete approximation for electromagnetics

In this chapter, different multi-field variational formulations for electrostatics and magnetostatics are proposed, which can provide an optimal discrete approximation of any particular vector field. The proposed formulations are constructed by appealing to mechanics point of view amenable to using general constitutive equations, which is quite different from electrostatics and magnetostatics formulations typical of physics and electrical engineering focusing on the corresponding global form suitable only for the linear case. In particular, the formulations we propose can be combined with mixed discrete approximations that can ensure the continuity of tangential component of the electric or magnetic field and normal component of electric displacement and magnetic flux even for low order interpolations. The choice of this kind is quite different from currently the favorite choice of high order finite element interpolations used for coupling electromagnetism with mechanics. The discrete approximation is based upon Whitney's interpolations representing the vector fields in terms of corresponding differential forms, with electric and magnetic fields as one-form and electric displacement and magnetic flux as two-form. The implementation of interpolations of this kind is made for 3D tetrahedron elements with non-standard approximation parameters defined not only at vertices (for zero-form) but at edges (for one-form) and facets (for two-form). The results of several numerical simulations are presented to illustrate the performance of different formulations proposed herein. In this section we present our main objectives of generalizing the electromagnetics formulation for the case that can include heterogeneous materials, and further account for combined constitutive behavior with mechanics. Hence, we seek to establish a clear analogy with mechanics point-of-view in formulating the boundary value problems in electrostatics in the vein of the mechanical problem formulation. Consider a domain Ω, with two different parts distinguished: the interior of domain Ω and its boundary Γ. In addition, the boundary can be separated into two different parts, depending on the primary or dual variable imposed. For electrostatics, on Γ V we impose the essential or Dirichlet boundary condition through scalar electric potential V , and on Γ D the natural or Neumann boundary condition through the normal component of D [START_REF] Balanis | Advanced engineering electromagnetics[END_REF], [START_REF] Jackson | Classical electrodynamics[END_REF]. The complete set of equations that describe the electrostatic phenomena for such domain can then be written as follows:

Field equations:

Kinematics E = -∇V Constitutive D = E Equilibrium ∇ • D = ρ f q    in Ω
Boundary conditions:

Essential BC V = V in Γ V Natural BC D • n = D in Γ D (4.1)
We indicated in the first equation in (4.1), that a "kinematics" equation links V and E in the equivalent manner as displacements are related to strains in mechanics. This equation automatically implies the electrostatic version of Faraday's law (2.8) 1 . The second equation in (4.1) shows the constitutive relation between E and D through the permittivity tensor, in a similar way that stresses are related to strains through the elasticity tensor. The final equation balancing the divergence of D with ρ f q corresponds to mechanics equilibrium equation relating the divergence of the stress tensor with the volume force.

A complementary problem formulation in terms of dual variables can also be posed. It will be illustrated for magnetostatics by using magnetic vector potential A, equivalent to Airy's function for the stress tensor σ. The corresponding complementary single-field formulation can be written as:

Field equations:

Compatibility B = ∇ × A Constitutive H = νB Equilibrium ∇ × H = J    in Ω
Boundary conditions:

Essential BC A = A in Γ A Natural BC n × H = J in Γ H (4.2)
where ν = µ -1 is the reluctancy of the material. The key difference from related problems in mechanics is the choice of vector potential and the curl operator instead of the gradient to define B . Such a choice allows to automatically satisfy (2.9) 2 .

Variational formulations

In this section, three variational formulations of the electrostatic problem are presented providing the basis for the corresponding choice of the finite element discrete approximation. The first two are single-field formulations: one is based on total potential energy with scalar potential V as the main variable, and the other is based on the total complementary potential energy with D. A mixed formulation is proposed as well based on the Hellinger-Reissner principle for mechanics adapted to electrostatics with two independent variables V and D. A magnetic formulation is provided in the last place to highlight the analogy between the complementary energy formulation of electrostatic and the corresponding formulation for magnetostatics.

Single field electrostatic formulation with scalar potential

The simplest approach for the electrostatic problems is by using a single independent field as electric scalar potential V [START_REF] Hammond | Dual finite-element calculations for static electric and magnetic fields[END_REF], [START_REF] Vogel | On some mixed variational principles in electroelastostatics[END_REF]. The remaining fields are obtained by using the strong form equations (4.1) 1 and (4.1) 2 . Namely, the weak form of balance equation in (4.1) 3 leads to variational formulation of the electrostatics:

Π V (V ) = ˆΩ 1 2 ∇V • ∇V -V ρ f q dΩ + ˆΓD V D dΓ (4.3)
where D is the electric displacement imposed at the Neumann boundary. Any trial functions must satisfy the Dirichlet boundary condition with V | Γ V = V , which should be applied at least at one point in order to guarantee a unique solution. Minimizing this potential leads to the corresponding variational equation:

G V V (V ; V * ) • • = ˆΩ ∇V * • ∇V dΩ -ˆΩ V * ρ f q dΩ + ˆΓD V * D dΓ = 0 (4.4)
where V * is the variation of V , which should take zero value at the Dirichlet boundary

V * | Γ V = 0.
The standard choice of finite element discrete approximation is suitable for this variational formulation (Moreno-Navarro et al., 2017).

Hellinger-Reissner hybrid formulation for electrostatics

It is possible to further weaken one of the equations in (4.1) in order to establish the variational formulation of the problem by using two independent fields, V and D. This is the Hellinger-Reissner type of variational formulation [START_REF] Vogel | On some mixed variational principles in electroelastostatics[END_REF], for which the mechanics equivalent is obtained featuring the complementary energy, which can be written as:

Π HR (V, D) = ˆΩ - 1 2 -1 D • D -∇V • D -V ρ f q dΩ + ˆΓD V D dΓ (4.5) SECTION 4.2

. Variational formulations

Keeping these two fields V and D independent allows to improve representation of the electric displacement field with respect to the one obtained from the single field formulation in (4.3) above. The variational equations corresponding to this potential can be obtained as:

G HR V (V, D; V * ) • • = -ˆΩ ∇V * • D dΩ -ˆΩ V * ρ f q dΩ + ˆΓD V * D dΓ = 0 ; G HR D (V, D; D * ) • • = -ˆΩ D * • -1 D dΩ -ˆΩ D * • ∇V dΩ = 0 (4.6)
The first variational equation recovers the strong form of the Gauss' law in (4.1) 3 while the second recovers the constitutive relation in (4.1) 2 .

Complementary energy dual formulation with vector potential for electrostatics

With this type of formulation, we use two different vector potentials, U and S, in order to fully define the electric displacement field [START_REF] Hammond | Dual finite-element calculations for static electric and magnetic fields[END_REF], [START_REF] Ren | A 3d vector potential formulation using edge element for electrostatic field computation[END_REF], [START_REF] Razek | Computation of 3d electrostatic local fields and forces using complementarity of dual formulations. application for capacitance and torque calculations in micromotors[END_REF]. Namely, we compute D as the following decomposition:

D = S + ∇ × U (4.7)
This decomposition allows to replace both potentials in (2.8) 2 and then treat (2.8) 1 in weak form. It is possible to establish the variational formulation for this total complementary potential energy:

Π D (U, S) = - ˆΩ 1 2 -1 ∇ × U • ∇ × U + 1 2 -1 S • S + -1 ∇ × U • S dΩ + ˆΓV V S • n dΓ - ˆΓV n × E • U dΓ (4.8)
Computing the corresponding variations with respect to the potentials, the weak form can be obtained as a result:

G D U (U, S; U * ) • • = -ˆΩ ∇ × U * • -1 ∇ × U dΩ -ˆΩ ∇ × U * • -1 S dΩ - ˆΓV n × E • U * dΓ = 0 ; G D S (U, S; S * ) • • = -ˆΩ S * • -1 S dΩ -ˆΩ S * • -1 ∇ × U dΩ + ˆΓV V S * • n dΓ = 0 (4.9)
with E and V as the imposed electric field and electric scalar potential.

Complementary energy dual formulation with single vector potential for magnetostatics

The dual formulation of this kind is yet easier to construct for the magnetostatic problem defined in (4.2), since no source term will appear [START_REF] Penman | Complementary and dual energy finite element principles in magnetostatics[END_REF]. The complementary energy variational formulation can then be written as:

Π A (A) = ˆΩ 1 2 µ -1 ∇ × A • ∇ × A dΩ -ˆΩ A • J dΩ + ˆΓH A • J s dΓ (4.10)
The weak form is obtained by minimizing this variational formulation with respect to the single field to obtain:

G A A (A; A * ) • • = ˆΩ ∇ × A * • µ -1 ∇ × A dΩ -ˆΩ A * • J dΩ + ˆΓH A * • J s dΓ (4.11)
where A * is the corresponding variation.

Finite element discrete approximation with Whitney's element basis

The main goal of Whitney's elements is to construct discrete approximation of vector fields by using differential forms [START_REF] Bamberg | A course in mathematics for students of physics[END_REF]. These differential forms, yet called i -forms (where i = 0, 1, 2 or 3), are linear functionals of vector fields which can be used to fully define the particular vector field. The unknown variables will the be the values of these differential forms, to be used as the finite element method degrees of freedom. The main advantage of using this kind of approximation is the ability of constructing them in somewhat intrinsic manner, regardless of particular choice of reference frame. Another big advantage of using differential forms is that every degree of freedom remains a scalar, even though the discrete approximation represents a vector field. Thus, the interpolation functions ought to be of vectorial character, which allows to preserve boundary conditions in a natural way. More importantly, such discrete approximation can be used to enforce the corresponding continuity across element boundary and thus improve the result accuracy.

The preferred formulation in mechanics in terms of vector fields [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], can be recast in terms of differential forms. The latter are linear functionals of the vector field, which are associated with both geometric and physical entities. Thus, for the discrete approximation constructed by the finite elements, we can choose:

• 0-form associated with vertices, which preserves point continuity.

• 1-form associated with edges, which preserves continuity of the electric field tangential component:

e i = ˆC E • dl (4.12)
where dl is the vector associated to element edge C.

• 2-form associated with faces, which preserves continuity of the electric displacement normal component over the facet:

d i = ˆS D • n dS (4.13)
where n is the unit exterior normal and dS is an infinitesimal element of the surface.

• 3-form associated with volumes.

In order to recover the continuum vector fields from differential forms, special interpolation functions are used. These functions are scalars for 0-and 3-forms, and vectors for 1-and 2-forms. In order to reconstruct a particular field, a linear combination of all the corresponding geometric entities is established for every element. For instance, for the voltage, the electric field and the electric displacement:

V = nv a=1 0 N a v a ; E = ne a=1 1 N a e a ; D = n f a=1 2 N a d a (4.14)
where 0 N a , 1 N a and 2 N a are respectively, 0-, 1-and 2-form interpolation functions for node, edge and facet a, whereas n v is the number of vertices, n e is the number of edges and n f is the number of facets per element.

The choice of differential forms is especially suitable for describing a unified mathematical structure of Maxwell's equations in terms of the exterior derivative (e.g. [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF]). Here, we use an operator defined for each i -form to obtain (i + 1)-form, which is expressed as follows: 0-form

∇ -→ 1-form ∇× --→ 2-form ∇• -→ 3-form (4.15)
We show in Fig. 4.1 operators acting on the different i -forms and the relations among them. Such a scheme is referred to as Tonti's diagram, e.g. [START_REF] Bossavit | Discrete magneto-elasticity: a geometrical approach[END_REF]. The exterior derivative for electrostatics with scalar potential formulation reduces to: The use of differential forms can be represented by discrete approximation by using Whitney's finite elements (e.g. (Bossavit, 1988a)), which have the advantage to unify the traitment of line, surface and volume integral in the sense of partition-of-unity. In Fig. 4.2, we show the isoparametric reference element used for tetrahedral mesh, where we indicate the different locations of nodes, edges and faces and their corresponding orientations for Whitney's element. There are four vertex, six edge and four facet unknowns for each Whitney's tetrahedral element. Center of the reference frame {ξ, η, ζ} is located at local node number 1 and each side starting at this node is of unit length. The interpolation functions for 0-form are:

V D 0 3 1 2 E D = ǫE ∇ ∇• ρ f q
0 N 1 = 1 -ξ -η -ζ ; 0 N 2 = ξ ; 0 N 3 = η ; 0 N 4 = ζ (4.16)
where the subscript a = 1, . . . , 4 is the local number of node and the superscript "0" denotes the corresponding index for a particular i -form. Each of these interpolation functions take unit value on their corresponding node and 0 at all other nodes, which make the nodal value at every node independent of the others. This will guarantee the field continuity between adjacent elements; Finally, these shape functions are differentiable inside the element domain, and guarantee the partition-of-unity property:

nv i=1 0 N i (ξ, η, ζ) = 1 ∀(ξ, η, ζ) ∈ Ω e (4.17)
This isoparametric element is used to construct the different tetrahedron elements in the mesh by distorting the reference element. Interpolation functions for 0-form are used to create the mapping from reference coordinates to physical coordinates through:

x(ξ) = nv i=1 0 N i (ξ)x i ⇒ x(ξ, η, ζ) = nv i=1 0 N i (ξ, η, ζ)x i y (ξ, η, ζ) = nv i=1 0 N i (ξ, η, ζ)y i z(ξ, η, ζ) = nv i=1 0 N i (ξ, η, ζ)z i (4.18)
where x i , y i and z i are the corresponding physical coordinates for every element node. The isoparametric mapping of element geometry can be characterized by the corresponding Jaco-bian matrix j , which can be written as:

j =        
        =    (x 2 -x 1 ) (y 2 -y 1 ) (z 2 -z 1 ) (x 3 -x 1 ) (y 3 -y 1 ) (z 3 -z 1 ) (x 4 -x 1 ) (y 4 -y 1 ) (z 4 -z 1 )    (4.19)
We note that the Jacobian matrix for a tetrahedron element has constant entries that depend only on the element nodal coordinates. The Jacobian matrix provides direct connection between the gradients with respect to natural and the one with respect to physical coordinates; denoting the latter as ∇ x and the former as ∇ ξ , we can easily show (e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]) that:

∇ x = j -1 ∇ ξ (4.20)
In order to create the gradients of the 0-form interpolation functions defining the electric field from scalar potential, we make use of this new definition of the nabla operator in order to write the corresponding discrete approximation:

0 B i = j -1 ∇ ξ 0 N i (4.21)
where 0 B i is the gradient of 0-form interpolation functions for node i .

The edge interpolation functions are constructed as the linear combination of 0-forms and its gradients, resulting with vector values. For example, for an edge a between nodes i and j, such a vector base function is defined as:

1 N a = 1 N i→j = ϕ a 0 N i 0 B j -0 N j 0 B i (4.22)
where two different notations are adopted, referring to either edge a or to the edge from node i to node j, which illustrates the edge orientation more clearly. The incidence coefficient ϕ a takes values either +1 or -1, depending on the agreement between the edge orientation on the reference element (see Fig. 4.2) with the corresponding edge in the global mesh. The convention sign for global edge orientation is taken arbitrarily. In this work, the unique global orientation is obtained by following the global node numbers: for an edge between nodes with global numbers i and j, the positive edge orientation goes from i to j when i < j and negative in the opposite case. Thus, it follows that 1 N i→j = -1 N j→i .

By following this convention, we can write explicitly 1-form interpolation functions for the six edges of a tetrahedron element:

1 N 1 = 1 N 1→2 = ϕ 1 j -1 (1 -η -ζ, ξ, ξ) T ; 1 N 2 = 1 N 1→3 = ϕ 3 j -1 (η, 1 -ξ -ζ, η) T ; 1 N 3 = 1 N 1→4 = ϕ 5 j -1 (ζ, ζ, 1 -ξ -η) T ; 1 N 4 = 1 N 2→3 = ϕ 2 j -1 (-η, ξ, 0) T ; 1 N 5 = 1 N 2→4 = ϕ 4 j -1 (ζ, 0, ξ) T ; 1 N 6 = 1 N 3→4 = ϕ 6 j -1 (0, -ζ, η) T (4.23)
As depicted in (4.15), the exterior derivative of 1-form is the curl operator and consequently, the derivation of the curl of these 1-forms 1 C a is needed. This definition is obtained by applying the curl operator to 1 N a :

1 C a = 1 C i→j = ∇ x × 1 N i→j = j -1 ∇ ξ × 1 N i→j = 2ϕ a 0 B i × 0 B j (4.24)
where ϕ a as the incidence coefficient and it is defined as for the 1-form. This convention is in agreement with the right-hand screw rule since the resulting flow is solenoidal with respect to edge a.

The facet interpolation functions are a combination of 0-form and the cross product of the gradients. Facet a is defined with positive orientation if the set of nodes {i , j, k} defines the flow normal to the facet by the right-hand rule. The corresponding vector base functions are constructed as:

2 N a = 2 N i→j→k = 2 0 N i 0 B j × 0 B k + 0 N j 0 B k × 0 B i + 0 N k 0 B i × 0 B j (4.25)
These the facet interpolations can be described as the combination of 0-forms and curls of 1-form, following the idea of the mathematical structure to define (i + 1)-forms:

2 N a = 2 N i→j→k = 0 N i 1 C j→k + 0 N j 1 C k→i + 0 N j 1 C k→i (4.26)

Whitney's element implementation

The finite element code chosen for this numerical implementation is FEAP [START_REF] Taylor | Feap-a finite element analysis program[END_REF]. The mesh is prepared in pre-processing to allocate all kinds of variables in different kinds of nodes which leads to the enhancement of the standard computer architecture not with respect to the choice of degrees of freedom but with respect to the corresponding shape functions (using different shape functions for different kind of nodes). Taking advantage of the 14-node isoparametric elements, the nodes are classified into vertex, edge and facet nodes as sketched in Fig. 4.3, where the numbering of the local nodes for the isoparametric tetrahedra is also displayed. Depending on the geometrical entity, the corresponding degree of freedom is put on nodes 1 to 4, 5 to 10 or 11 to 14. This ensures the continuity of these variables among elements. For the chosen Whitney's interpolations, the position of the nodes does not matter as long as the node is placed on the corresponding edge or facet.

Implementation of single-field scalar potential for electrostatics

The definition of the discrete approximations for scalar potential, put on vertices according to [START_REF] Ren | A 3d vector potential formulation using edge element for electrostatic field computation[END_REF] and [START_REF] Razek | Computation of 3d electrostatic local fields and forces using complementarity of dual formulations. application for capacitance and torque calculations in micromotors[END_REF]:

V ≈ 4 i=1 0 N i (x)v i (4.27)
where v i are the degrees of freedom or unknowns that the code is going to calculate. Introducing this approximation into the weak form (4.4), we obtain the residual corresponding to the scalar potential, which is written in matrix form:

R V a • • = ˆΩ 0 B T a 4 i=1 0 B i v i dΩ -ˆΩ 0 N a ρ f q dΩ + ˆΓD 0 N a D dΓ ; a = 1, 2, 3, 4 (4.28)
By consistent linearization of this residual, we obtain the corresponding tangent matrix, which can also be written in matrix notation:

K V ab = -ˆΩ 0 B T a 0 B b dΩ ; a, b = 1, 2, 3, 4 (4.29)
Note that this matrix and the residual pertain only to vertex nodes a and b, whereas they are equal to zero for edge and facet nodes. The central problem to solve for the unknown nodal values of scalar potential can be written as:

Kv = f (4.30)
where f is the force vector. For this formulation, each entry of this vector is defined as:

f a = -ˆΩ 0 N a ρ f q dΩ + ˆΓD 0 N a D dΓ (4.31)

Implementation of Hellinger-Reissner formulation

The same methodology is taken for Hellinger-Reissner formulation, this time introducing a new approximation for the electric displacement:

D ≈ 14 i=11 2 N i (x)d i (4.32)
where d i are the degrees of freedom associated with D. Introducing this approximation in (4.6) the corresponding residuals can be obtained: 11,12,13,14 ;

R D a = -ˆΩ 2 N T a -1 14 i=11 2 N i d i dΩ -ˆΩ 2 N T a 4 j=1 0 B j v j dΩ ; a =
R V a = -ˆΩ 0 B T a 14 i=11 2 N i d i dΩ -ˆΩ 0 N a ρ f q dΩ + ˆΓD 0 N a D dΓ ; a = 1, 2, 3, 4 (4.33)
Then, linearizing the previous equations, the problem to solve is:

K d v = 0 f (4.34)
where the stiffness matrix K can be split into four different sub-matrices:

K ab = K DD ab K DV ab K V D ab 0 (4.35)
The entries of that stiffness matrix are:

K DD ab = ˆΩ 2 N T a -1 2 N b dΩ ; a, b = 11, 12
, 13, 14 11, 12, 13, 14 ; b = 1, 2, 3, 4 ; 2, 3, 4 ; b = 11, 12, 13, 14 (4.36) And the force vector f for node a is:

K DV ab = ˆΩ 2 N T a 0 B b dΩ ; a =
K V D ab = ˆΩ 0 B T a 2 N b dΩ ; a = 1,
f a = -ˆΩ 0 N a ρ f q dΩ + ˆΓD 0 N a D dΓ ; a = 1, 2, 3, 4 (4.37)

Implementation of complementary energy dual formulation for electrostatics with two vector potentials

Introducing the split and separate discrete approximations for U and S, the definition of D can be obtained as:

D ≈ 10 i=5 1 C i (x)u i + 14 j=11 2 N j (x)s j (4.38)
where u i and s j are the degrees of freedom or unknowns that we are going to compute. Note that even though the potentials U and S are vectors, their corresponding degrees of freedom (as differential forms) are scalars, since we will use the vector-type interpolations (namely, one-forms). Moreover, the orientations of chosen unknown differential forms are defined by edge-curl of 1-forms and face shape functions 2-forms.

Introducing the last definition into the weak forms above, the final residuals can be written in matrix notation; namely, for D-form: 11, 12, 13, 14 (4.39) Furthermore, the stiffness matrix can be split in different sub-matrices for D-form:

R U a = ˆΩ 1 C T a -1 10 i=5 1 C i u i dΩ + ˆΩ 1 C T a -1 14 j=11 2 N j s j + ˆΓ1 1 N T a (n × E 0 ) dΓ ; a = 5, 6, . . . , 10 R S a = ˆΩ 2 N T a -1 10 i=5 1 C i u i dΩ + ˆΩ 2 N T a -1 14 j=11 2 N j s j dΩ - ˆΓ2 2 N T a n V dΓ ; a =
K ab = K UU ab K US ab K SU ab K SS ab (4.40)
The entries of this stiffness matrix are: 5, 6, . . . , 10 ; b = 11, 12, 13, 14 11, 12, 13, 14 ; b = 5, 6, . . . , 10 11, 12, 13, 14 (4.41) And the components of the force vector f are defined as:

K UU ab = -ˆΩ 1 C T a -1 1 C b dΩ ; a, b = 5, 6, . . . , 10
K US ab = -ˆΩ 1 C T a -1 2 N b dΩ ; a =
K SU ab = -ˆΩ 2 N T a -1 1 C b dΩ ; a =
K UU ab = -ˆΩ 2 N T a -1 2 N b dΩ ; a, b =
f 1,a = ˆΓ1 1 N T a n × E dΓ f 2,a = - ˆΓ2 2 N T a n V dΓ (4.42)

Gauging a vector potential and boundary conditions for complementary energy formulation

The previous formulation makes use of a vector potential U. The rotational part of the electric displacement is obtained from such potential. To provide a specific value of electric displacement D, there is not a unique possibility of U, since any irrotational potential to be added to this potential would not affect the electric displacement value, but only act as an integration constant. This non-uniqueness can be removed if the vector potential U is gauged. A previous work in [START_REF] Stark | On the boundary conditions for the vector potential formulation in electrostatics[END_REF] has a complete relation of different references on the existing methods to gauge such potential. The method chosen in this thesis is the identification of a tree set. This method is based on the topological aspect of the discretization employed to define a set of linearly dependent differential forms, which have to be eliminated in order to return a unique vector potential. In the finite element method, the removal of those variables implies prescribing the degrees of freedom on the tree set corresponding to the edges of the tetrahedral mesh.

The details for identifying one of the tree sets are given in [START_REF] Hale | A logic for identifying the trees of a graph[END_REF], including a method to determine all possible tree combinations through the incidence matrix. A tree set must contain n n -1 edges, and they cannot close a surface. For simplicity, the value imposed in those edges is zero.

In this thesis, an algorithm to select a tree set has been designed based on a matrix that contains as many rows as edges. In each row, there are the initial and final global numbers of the vertex nodes of the corresponding edge. The method begins evaluating the two global node numbers of the first row. If they are different, the edge is accepted in the tree set. Then, the first node number is replaced by the second one every time it appears in another matrix entry. The method is repeated for every row until the tree has n n -1 branches. If the mesh is regular, like the one drawn in Fig. 4.4, a more specific and systematic algorithm to select the tree can be used, as stated in [START_REF] Kameari | Three dimensional eddy current calculation using edge elements for magnetic vector potential[END_REF]. This algorithm starts selecting all the edges in the line defined by the intersection of the planes y = z = 0. Then, in the plane z = 0, all the edges in y direction that start in every vertex node of the previous line are picked. Finally, every vertical line starting from all vertex nodes in plane z = 0 are selected as well.

The complementary energy usually provides superior accuracy of the electric displacement, but presents a difficulty with respect to the standard V-formulation in not readily providing all other fields of interest, such as voltage potential. Moreover, an additional difficulty concerns any problem where both Dirichlet and Neumann boundary conditions have to be imposed. Namely, with the chosen differential forms as unknowns, a particular edge, the topological entity that contains the variable to solve, leads to an ambiguous situation that it can belong to both kinds of boundary surfaces, Dirichlet and Neumann. The difficulty pertains to deciding for every boundary face which edge should take care of the Dirichlet boundary condition. Moreover, for the remaining edge that handles the Neumann boundary, the residual term ought to be distributed over the edges of the triangle that are not considered in the Dirichlet boundary.

Implementation of dual formulation of magnetostatics with single vector potential

A similar approach to the last formulation for magnetostatics can be taken by introducing the chosen approximations for magnetic vector potential A:

A ≈ ne i=1 1 N i (x)a i (4.43)
where a i are the degrees of freedom or unknowns that the code is going to calculate. Introducing the last definition into (4.11), the corresponding residuals can be written as:

R A a = ˆΩ 1 C T a µ -1 ne i=1 1 C i a i dΩ -ˆΩ 1 N T a J dΩ + ˆΓH 1 N T a J s dΓ ; a = 5, 6, . . . , 10
(4.44)
By consistent linearization of this residual, we obtain the system similar to the one in (4.30). The corresponding tangent matrix can be written in matrix notation: 5, 6, . . . , 10 (4.45) and the force vector f component is defined as:

K A ab = -ˆΩ 1 C T a µ -1 1 C b dΩ ; a, b =
f a = -ˆΩ 1 N T a J dΩ + ˆΓH 1 N T a J s dΓ (4.46)

Numerical simulations and validation tests

In this section, we present the results of several numerical simulations that can illustrate performance of different formulations. First, we choose several validation test examples, which are compared against known analytic solution. A more complex test case where analytic solution is not available is presented afterwards, including the results of some practical applications solved with different formulations.

Validation tests

The following numerical examples are called validation tests due to the fact that an analytic solution exists for each one of them. This is a necessary condition to satisfy in order to prove that the different formulations can all solve the problems where both analytical and numerical results fully agree. However, this is not always fully sufficient, if the problem does not allow to test each and every aspect of the proposed formulation.

In order to find an analytic solution, we assume that E is the negative gradient of V and the constitutive relation for D in (2.4) are written for isotropy and that the permittivity tensor can be simplified to I. With these hypotheses, the electrostatic problem can be fully defined in terms of the scalar potential:

∇ 2 V (x, y , z) = -ρ f q (4.47)
This equation is the well known Poisson equation, which reduces to the Laplace equation when ρ f q = 0. The latter is much easier to solve and represents a sort of patch test.

Laplace equation, imposed scalar potential

The geometry of this first example is a cube that has the following boundary conditions turning into a one-dimensional problem:

∂V ∂x x=0 = ∂V ∂x x=lx = ∂V ∂y y =0 = ∂V ∂y y =ly = 0 ; V (x, y , 0) = 0 ; V (x, y , l z ) = V z (4.48)
The lateral faces do not allow electric displacement flowing out of the cube, whereas top and bottom faces have imposed and constant scalar potentials. The solution for the scalar potential:

V (z) = V z l z z (4.49)
V is a linear distribution in z. Hence, from this expression it is possible to calculate the electric displacement as:

D = -∇V =   0 0 -V z /l z   (4.50)
The only non-zero component is the third one taking a constant value. The selected dimensions for this example are l x = l y = l z = 2 × 10 -3 m, the imposed V z = 20 V and = 15 × 10 -12 F/m.

As can be seen in Fig. 4.5, the computed distribution of voltage potential is indeed linear with bottom and top values 0 and 20 V, respectively. The computed value for D z = -1.5 × 10 -7 C/m 2 is homogeneous as expected from (4.50).

3.33E+00 6.67E+00 Laplace equation, imposed scalar potential and electric displacement

1.00E+01 1.33E+01 1.67E+01 -1.16E-15 2.00E+01 VOLTAGE (V) -1.60E-07 -1.54E-07 -1.49E-07 -1.44E-07 -1.39E-07 -1.65E-07 -1.33E-07 ELEC. DISP. Dz (C/m²) 3.33E+00 6.67E+00 1.00E+01 1.33E+01 1.67E+01 -1.16E-15 2.00E+01 VOLTAGE (V) -1.60E-07 -1.54E-07 -1.49E-07 -1.44E-07 -1.39E-07 -1.65E-07 -1.33E-07 ELEC. DISP. Dz (C/m²)
In this example, imposing the Neumann boundary condition is the main difference from the previous example. Thus, there is only a slight variation in the chosen boundary conditions with respect to the previous example:

∂V ∂x x=0 = ∂V ∂x x=lx = ∂V ∂y y =0 = ∂V ∂y y =ly = 0 ; V (x, y , 0) = 0 ; ∂V ∂z = -D z (4.51)
This computed solution can still be simplified to unidimensional and linear in V :

V (x, y , z) = -D z z (4.52)
The material properties and measures are the same as in the last example, with D z = 1.5 × 10 -7 C/m 2 . This gives V (x, y , l z ) = -20 V on top and a constant electric displacement of value D z = 1.5 × 10 -7 C/m 2 as shown in Fig. 4

.6. -1.67E+01 -1.33E+01 -1.00E+01 -6.67E+00 -3.33E+00 -2.00E+01 2.95E-15 VOLTAGE (V) 1.40E-07 1.45E-07 1.50E-07 1.55E-07 1.60E-07 1.35E-07 1.65E-07 ELEC. DISP. Dz (C/m²) -1.67E+01 -1.33E+01 -1.00E+01 -6.67E+00 -3.33E+00 -2.00E+01 2.98E-15 VOLTAGE (V) 1.40E-07 1.45E-07 1.50E-07 1.55E-07 1.60E-07 1.35E-07 1.65E-07 ELEC. DISP. Dz (C/m²)

Poisson equation, imposed scalar potential

In this validation case, a source term is activated to induce non-homogeneous field values. The boundary conditions considered in this section are the same as those considered in the first numerical example with:

∂V ∂x x=0 = ∂V ∂x x=lx = ∂V ∂y y =0 = ∂V ∂y y =ly = 0 ; V (x, y , 0) = 0 ; V (x, y , l z ) = V z (4.53)
For a problem of this kind, the solution can be obtained as a superposition of a Poisson equation solution with homogeneous boundary conditions and a Laplace problem taking into account those conditions. The final solution can be expressed as:

V (x, y , z) = V z l z z + ∞ p=1,3,... 4 ρ f q l 2 z p 3 π 3 sin pπz l z (4.54)
where 50 terms in the summation have been taken into account, ρ f q = 0.01 and the other coefficients are the same as in Section 4.4.1 except V z = 80 V, taken higher to see clearly the asymmetry of the problem. In Fig. 4.7 left, it can be appreciated in the voltage distribution that the maximum is concentrated a little bit above the middle vertical section with a theoretical maximum value of approximately 374 V, obtained from the maximization of (4.54). The closer to this theoretical value, for the mesh plotted with 1000 elements, is V-formulation. The distribution of D z can be obtained by deriving (4.54), and the summation gives an almost linear response with values at the bottom of D z = -1.06 × 10 -5 C/m 2 and at the top of D z = 9.37 × 10 -6 C/m 2 .

One can also notice the existence of slight oscillations near the edges of the mesh in Fig. 4.7 right. In the case of V-formulation, the electric displacement can only be constant for every element. Taking into account that the analytical solution is linear, the exact solution can never be reached, although the more refined the mesh is, the closer to the analytical solution will be. The number of elements surrounding a particular point in the mesh is essential to the plot since the stresses plotted in the figure are smoothened. The edges are influenced by half of the elements within the same row. Regarding H-R, even though D field can be indeed linear, it is limited by (4.6) as it has to be compatible with the gradient subspace of voltage. Therefore, the oscillations now appear along all iso-stress lines since tetrahedra within the same row now have different values.

In Fig. 4.8, a convergence study on the energy by mesh refinement has been performed. Even though V-formulation is closer to the maximum expected value of V , H-R is closer than V-formulation to the analytical solution for energy (1/2 E • D = 8.98 × 10 -9 J) at an equal number of elements, as could be expected from a mixed method, improving the accuracy of the electric displacement calculation.

Tubular geometry

This example is also among validation examples, but somewhat more demanding with respect to element distortion, resulting from a change in geometry turning the parallelepiped shape into a hollowed cylinder or a tube. This modification allows to determine if there is any mesh dependence of the solution with respect to element distortion. Also, it allows to fully verify proposed formulation for distorted mesh.

z x y V 1 V 2 1 2 1 or D Figure 4
.9: Schematic geometry representation of a quarter of the hollowed cylinder. Symmetry planes in x = 0 and y = 0. Dimensions in cm.

The geometry of the problem is defined in Fig. 4.9 along with the boundary conditions and dimensions. Only a quarter of the geometry is represented since it is an axisymmetric problem. On the interior lateral face, a voltage V 1 = 0 V is imposed; on the exterior lateral face, either voltage V 2 = 20 V or D • n = D is imposed depending on the kind of problem. An electric charge density ρ f q is imposed in the material as well for the Poisson case. The properties for the material of this example remains the same as the previous ones except for the permittivity, in agreement with the dielectric material of the next case = 4 0 I. The analytical solution for the proposed problem can be obtained by taking into account the axisymmetry and cylindrical coordinates. Thus, both Laplace and Poisson equations will now be affected by the change of coordinates. Since the solution V (r ) only depends on the radius, the equation to solve reduces to an Euler-Cauchy ordinary differential equation with a particular term in the case of the Poisson equation:

r 2 d 2 V dr 2 + r dV dr = r 2 ρ f q (4.55)
The solution to the previous equation is of the form:

V (r ) = C 1 ln (r ) + C 2 + ρ f q r 2 4 (4.56)

Parallel plate capacitor simulation

This numerical example simulates a parallel plate capacitor that consists of two electrodes, represented by conductor plates charged, and a dielectric cylinder with = 4 0 I. Fig. 4.12 shows the dimensions for a quarter of the geometry since symmetry conditions on planes x = 0 and y = 0 allow to reduce the problem.

Despite the fact that in the graphical illustration in this figure only two different solid materials are presented, for the actual computation we have to introduce vacuum or air as a third material connecting the two solid materials. Hence, the permittivity of this added material is that of the vacuum 0 . We note in passing that such a solution is typical of electromagnetic spectrum, which does not need matter to propagate, and yet the finite element method does need the mesh with a connection between nodes to solve such problem.

x y z V = 3000 V 2 0.6 1 Figure 4.12: Schematic geometry representation of a quarter of the parallel plate capacitor; two charged conductor plates and a dielectric cylinder. Symmetry planes in x = 0 and y = 0. Dimensions in cm.

Boundary conditions can be imposed numerically by setting all nodes of the two conductor plates to the corresponding potential: top to V = 3000 V and bottom to ground potential or V = 0 V. One can observe that the presence of the dielectric with a higher permittivity than that of the air and the fact that the conductor plates are finite, breaks the unidimensional nature of this problem. If the plates as well as the cylinder were infinite, the electric field and displacement would be completely regular and vertical. This problem has been extracted from [START_REF] Ren | A 3d vector potential formulation using edge element for electrostatic field computation[END_REF] to demonstrate that potential and complementary energy formulations, after refinement of the mesh, tend to the same solution. The purpose of this problem is to compare our formulations with an electrostatic example already solved, although in the cited reference there are not explicit results for electric field or displacement, but rather a qualitative representation of both fields. This representation is replicated in Fig. 4.13 indicating the flow orientation of electric field.

V-Formulation

Hellinger-Reissner

Electric displacement D Electric f ield E Figure 4.13: Electric field and displacement vector representation calculated for the parallel plate capacitor example. First row, V-formulation; second row, Hellinger-Reissner.

As expected, the main difference between E (first column) and D (second column) in that figure is that the presence of the dielectric cylinder with different permittivity creates a discontinuity in electric field as opposed to the continuity of electric displacement. This fact can be appreciated in the magnitude difference of the arrows inside the cylinder and outside. This presence of the dielectric creates also a distortion of both fields, bending them.

The arrows in the air material surrounding the capacitor bend around the two electrodes as it is expected from an electrostatic problem, from higher to lower potentials and the intensity of the fields lower as they come out from symmetry planes rather than external sides. The total energy of the system is represented in Fig. 4.14 for both V-formulation and Hellinger-Reissner formulations. It is important to note that both formulations approach the same system energy value when sufficient mesh refinement is carried out, although for this particular case, a finer mesh was impossible to generate due to our low computational resources. As the number of elements increase, energy values for V-formulation come from higher to lower as opposed to Hellinger-Reissner. This is a consequence of using the potential or the complementary energy for the calculations.

For the same number of elements, Hellinger-Reissner formulation takes more computation time due to the presence of an additional degree of freedom. However, the accuracy is higher energy-wise since no big improvement is obtained with higher number of elements. Therefore, for a low number of elements, one can conclude that Hellinger-Reissner formulation is preferable.

5

Multi-scale approach for coupled theory

In this chapter, we present the development of a 3D lattice-type model at micro-scale based upon the Voronoi-cell representation of material microstructure. This model can capture the coupling between mechanic and electric fields with non-linear constitutive behavior for both. More precisely, for electric part we consider the ferroelectric constitutive behavior with the possibility of domain switching polarization, which can be handled in the same fashion as deformation theory of plasticity. For mechanics part, we introduce the constitutive model of plasticity with the Armstrong-Frederick kinematic hardening and a localized failure for the first three modes that can represent the cracks. This model is used to simulate a complete coupling of the chosen electric and mechanics behavior with a multiscale approach implemented within the same computational architecture.

Basics of a lattice-type model

To create a model of this type, the domain to study is divided into polyhedral regions. The circumcenter1 of each polyhedron constitutes one of the ends of the lattice beam. Two adjacent Voronoi cells are held together with a lattice beam which is chosen as a beam [START_REF] Nikolić | Lattice element models and their peculiarities[END_REF], [START_REF] Bui | Enriched timoshenko beam finite element for modeling bending and shear failure of reinforced concrete frames[END_REF].

Figure 5.1: Left, two adjacent Voronoi cells held together by the corresponding lattice model. Right, the assumed shape of the beam and the approximation of the circular cross-section.

In Figure 5.1 left, two adjacent cells, and the corresponding lattice beam are pictured. As a property of the circumcenter, the beam is perpendicular to the joint face of these two cells. Also, the face intersects in the middle of the beam, dividing it into two halves. This is convenient to reproduce the fracture of the element, since a crack on the lattice beam means that two adjacent Voronoi cells are separating.

The transversal area associated with this beam is the area of the joint face, a polygon, as depicted in Figure 5.1 right. The transversal area is assumed to be an equivalent circle with the same area to simplify the calculation of the beam properties. The lattice beams define the dual mesh of the Voronoi cell. Delaunay algorithms for the creation of tetrahedron are used to create the mesh of beams that will be used to solve the problem.

Notation for beam model

The selected model for the beam is in agreement with Timoshenko hypotheses, where the transversal sections can be rotated with respect to the normal of the neutral line. The degrees of freedom for the mechanic field are three displacements and three rotations:

u = (u 1 , u 2 , u 3 , ϕ 1 , ϕ 2 , ϕ 3 ) T (5.1)
The generalized strain measures are ordered in a vector:

ε = (ε 1 , γ 12 , γ 13 , κ 1 , κ 2 , κ 3 ) T (5.2)
The corresponding internal force vector:

F = (N 1 , N 2 , N 3 , M 1 , M 2 , M 3 ) T (5.3)
Regarding the electric field, the only degree of freedom is the electric potential V . For the beam model, only the axial component x 1 is relevant; thus, by definition, the electric field E 1 and the electric displacement D 1 are scalars along any lattice. Assuming the homogeneity of any variable within every section, we can define the electric charge as Q 1 = D 1 A, with A as the cross-section area of the beam.

A point on the beam can also be represented in global coordinates as:

x = (x 1 , x2 , x3 ) T (5.4)
where the tilde overscript ∼ is used in this section to describe that the corresponding variable is expressed in the global orthonormal basis. This basis is composed of base unit vectors e G i for every direction i . The point can also be expressed in beam local coordinates as:

x = (x 1 , x 2 , x 3 ) T (5.5)

Notice that the overscript omitted from the notation means that the vector is expressed in local coordinates. This local basis is comprised of base unit vectors e i for every direction i . The axial direction of the beam is always the first one. Since the beam section is considered circular, transverse directions are taken randomly and normal to the axial direction, in agreement with the right-hand rule, as illustrated in Figure 5.2.

The rows of the transformation matrix T are the base unit vectors e i expressed in global coordinates:

T =   ← ẽT 1 → ← ẽT 2 → ← ẽT 3 →   (5.6)
This transformation matrix is used to change the basis in which every variable is denoted. In general, a tensor A of any order is expressed in global coordinates as:

Ãijkl... = T Ii T Jj T Kk T Ll . . . A IJKL...
(5.7)

Kinematic equations

Kinematic equations relate displacement and voltage degrees of freedom with generalized strain defining the state or primal variables. For a solid continuum model with the hypothesis on small displacements, the corresponding equations that define the infinitesimal strain tensor and the electric field vector are:

ε = 1 2 ∇ ⊗ u + (∇ ⊗ u) T = ∇ sy u E = -∇V (5.8)
To obtain the reduced the model for a Timoshenko beam, we need to take into account that every variable only changes in axial direction. Thus, the only relevant strain components are the axial strain ε and transversal shears γ 12 and γ 13 , as well as the axial electric field E 1 , each one constant within the beam cross-section. Besides, as a consequence of the additional degrees of freedom for the mechanic field, curvatures κ 1 , κ 2 , and κ 3 have to be calculated. These assumptions will further specialize when considering the Timoshenko kinematic model [START_REF] Bui | Enriched timoshenko beam finite element for modeling bending and shear failure of reinforced concrete frames[END_REF], the continuum version of (5.8) reduces to:

          ε 1 γ 12 γ 13 κ 1 κ 2 κ 3 E 1           =                      ∂ ∂x 1 0 0 0 0 0 0 0 ∂ ∂x 1 0 0 0 -1 0 0 0 ∂ ∂x 1 0 1 0 0 0 0 0 ∂ ∂x 1 0 0 0 0 0 0 0 ∂ ∂x 1 0 0 0 0 0 0 0 ∂ ∂x 1 0 0 0 0 0 0 0 ∂ ∂x 1                                u 1 u 2 u 3 ϕ 1 ϕ 2 ϕ 3 V          
(5.9)

Conservation equations

Conservation principles are used for defining the equilibrium in a material, which are defined in terms of partial differential equations. The solid version of these equations is expressing equilibrium of mechanical forces and electric field in the absence of external charges. For this beam model, the main hypotheses assumed are that body forces are neglected and a quasi-static approach is taken. Therefore, the solid version of the equations is:

σ ∇ = 0 ; ∇ • D = 0
(5.10)

The adaptation for the beam model includes the hypotheses in the previous section regarding one-dimensionality and homogeneity of variables in a transversal section. Thus, the simplified version of equilibrium equations for beams reduces to:

∂N i ∂x 1 = 0 ; ∂M i ∂x 1 = 0 ; ∂Q 1 ∂x 1 = 0
(5.11)

Constitutive equations for piezoelectricity

Constitutive equations will relate the primal variables with stress or flow variables as thermodynamic duals through the given defined material properties. In the case of piezoelectric materials, in addition to the conventional relationship between the same kind of variables, electrical variables affect the behavior of the mechanical ones and vice-versa. The constitutive equations can be derived from the free-energy potential as defined in [START_REF] Moreno-Navarro | Linear elastic mechanical system interacting with coupled thermo-electro-magnetic fields[END_REF], from which we can retrieve:

σ = Cε -e e E D = E + e e ε
(5.12)

The coefficients used in (5.12) are the fourth-order stiffness tensor C, the third-order piezoelectric tensor e e , and the second-order permittivity tensor , all of them obtained for assumed material properties of piezoelectricity.

To reduce such constitutive model to beams, we need to account for the simplification made in the last sections. First, only a few components of primal variables are relevant, so are the corresponding material tensor entries. Second, the addition of rotational degrees of freedom for the mechanic field introduces new structural variables into the stiffness tensor. Finally, only axial strain affects the axial electric variables and vice-versa. With these assumptions in hand, we can write the beam-model version of (5.12): (5.13) where

          N 1 N 2 N 3 M 1 M 2 M 3 Q 1           =            ȆA 0 0 0 0 0 -e 11 A 0 k c GA 0 0 0 0 0 0 0 k c GA 0 0 0 0 0 0 0 GJ 0 0 0 0 0 0 0 ȆI 0 0 0 0 0 0 0 ȆI 0 e 11 A 0 0 0 0 0 1 A                      ε 1 γ 12 γ 13 κ 1 κ 2 κ 3 E 1          
Ȇ = E (1 -ν) / [(1 + ν) (1 -2ν)]
is the first term of the stiffness tensor (with E as the Young's modulus and ν as the Poisson's coefficient), k c is the shear correction factor, G is the shear modulus, J is the polar moment of inertia, I is the moment of inertia, 1 is the permittivity in axial direction, and e 11 is the piezoelectric coefficient. Notice that inertia moments are the same for x 2 and x 3 directions, and k c ≈ 0.9 since the transversal area is considered a circle.

Ferroelectric model

The beams used in this chapter have a more elaborated constitutive model than the one presented in (5.13). This model introduces non-linear relations in mechanics (plasticity and softening) and electricity (ferroelectricity). A complete description of these constitutive models is presented herein.

Polarization P is a macroscopic magnitude that accumulates the microscopic electric dipole moments in a material. These dipoles are either generated in the presence of an electric field or as a consequence of a particular microstructure. While the former is proportional to the electric field applied, the latter is a permanent value, usually called remanent polarization P r , with the only possibility of a change of orientation or domain switching [START_REF] Balanis | Advanced engineering electromagnetics[END_REF].

In general, there are three types of macroscopic behavior regarding polarization: dielectric, paraelectric and ferroelectric. For the first one, the polarization is linear to the electric field applied; the second one depends on the electric field as well, but the dependency is non-linear; finally, ferroelectrics have a dependency on the electric field and, superposed, the remanent polarization causing hysteresis phenomena to arise. Ferroelectricity implies always a coupling with the mechanical field, since only piezoelectric materials can be ferroelectrics. In fact, they also have to be pyroelectrics, but we are not considering in this chapter that kind of coupling [START_REF] Said | Ferroelectrics and their applications[END_REF].

Below Curie Temperature

Above Curie Temperature The appearance of this remanent polarization can be explained from a microstructure point-of-view. In this chapter, only polarizable tetragonal cell materials are being discussed [START_REF] Jong | A database to enable discovery and design of piezoelectric materials[END_REF]. A representative unit cell of this type is sketched in Figure 5.3 for PbTiO 3 depending on the temperature. Below the Curie temperature T c , the microstructure is tetragonal, and it behaves as ferroelectric; above this temperature, the microstructure is cubic, the remanent polarization is gone, and its behavior is paraelectric.

The atoms in the center (Ti, in gray) and the corners (Pb, in black) are charged positively while the ones in the middle of the faces are charged negatively (O, in white). Above the Curie temperature, the cell is perfectly symmetric and cubic, with Ti in the center. Below the Curie temperature, the cell reaches minimum potential energy if Ti is ousted from the center in a direction towards one of the oxygen atoms. This eccentricity creates a microscopic electric dipole due to the misalignment compared with the cubic structure. Notice as well that the cell has been enlarged in the same direction of the misalignment since Ti repels the top O. If the position of Ti switches to any of the other five equilibrium positions, the remanent polarization and the cell enlargement will change as well.

A strong electric field has to be applied in the desired direction to switch the central atom into any of the other minimum positions. Also, this position is modified if a compressive stress is applied in the direction of the current remanent polarization or traction in any of the corresponding transversal directions. The angle between the old and the new P r vector determines the two types of switching: 180 • and 90 • . An oriented electric field can cause both switches, whereas stress can only generate a 90 • switch.

The selected criteria to determine whether a switch occurs is extracted from [START_REF] Keip | A ferroelectric and ferroelastic microscopic switching criterion for tetragonal ferroelectrics[END_REF] where the combination of both electric field and stress is taken into account. It is an energy criterion at a microscopic level, and several possibilities have to be evaluated at the same time for both kinds of switch. A switch can take place if the following conditions meet:

E • ∆P r 2E • P s ≥ 1 , for 180 • switch E • ∆P r E • ε s + σ • ∆ε r σ • ε s ≥ 1 , for 90 • switch (5.14)
where the increments ∆P r and ∆ε r are the tensors of change for remanent polarization and strain respectively, and P s and ε s are the spontaneous values of remanent polarization or strain induced by internal microstructure. The above conditions need to be evaluated for each of the five alternative directions of polarization.

This change in the microstructure alters the isotropy of the material, making it transversally isotropic with a preferential direction in the remanent polarization orientation. In [START_REF] Labusch | On the influence of ferroelectric polarization states on the magneto-electric coupling in two-phase composites[END_REF], the authors provide a comprehensive description of the enthalpy for magnetoelectro-mechanical polarizable materials at a microscale model. We obtain the corresponding constitutive equations through the derivation of the enthalpy with respect to the dual variables. Since this definition of enthalpy is quadratic, the constitutive coefficients are constant. A switch in polarization makes these coefficients to change their value.

For the beam lattice of our model, the coefficients of (5.13) are also subjected to change with every polarization switch as in the previous solid model, which implies a change in Ȇ, G, e 11 , and 1 . The factor that modifies these properties from its non-polarized to its polarized value is the polarization multiplier m through the following linear distributions: where the superscripts np and p stand for non-polarized and polarized values of each coefficient and |m| is the absolute value of m. The values for remanent polarization and remanent strain are also a function of this multiplier as follows:

Ȇ = Ȇnp + |m| Ȇp -Ȇnp G = G np + |m| (G p -G np
P r 1 = m P s 1 ε r 1 = ε np 1 + |m| ε p 1 -ε np 1 (5.16)
The manifestation of these remanent variables modify the constitutive equation (5.13) in the following way:

          N 1 N 2 N 3 M 1 M 2 M 3 Q 1           =            ȆA 0 0 0 0 0 -e 11 A 0 k c GA 0 0 0 0 0 0 0 k c GA 0 0 0 0 0 0 0 GJ 0 0 0 0 0 0 0 ȆI 0 0 0 0 0 0 0 ȆI 0 e 11 A 0 0 0 0 0 1 A                      ε 1 -ε r 1 γ 12 γ 13 κ 1 κ 2 κ 3 E 1           +           0 0 0 0 0 0 P r 1 A          
(5.17)

Two models to represent the switching of polarization and to define the polarization multiplier are presented here. For the first one, the only switch considered is the 180 • since beams can only take into account axial variations in voltage. Therefore, and since 90 • switches are not possible to simulate, the strain is not a factor for m, and a compressive force cannot depolarize the beam in axial direction.

A more refined version of this model requires a small variation to take into account the effect of compression and traction in the bar. The idea is to include a third state of the switch of no remanent polarization as the replacement for 90 • switching. This way, every beam can take three states: positive, negative, and no axial remanent polarization or zero polarization. This can be enforced as a constraint in the same manner as modeling the inextensible beam deformation [START_REF] Medić | Beam model refinement and reduction[END_REF].

A numerical difficulty arises when the switch in polarization is implemented as a Heaviside step once one of the conditions in (5.14) is met. This abrupt change can lead to oscillating residual norms in the finite element method computations due to the uncertainty of the beams getting simultaneously polarized. Here, a smooth correction to this Heaviside step is proposed for these two models through the hyperbolic tangent function. The idea is to make a bijective function between electric field and electric displacement, removing any possible uncertainty and introducing a smooth change of slope. Once the polarization switches, the constitutive model changes to represent hysteresis behavior.

Switching model 1

The flowchart for the switching model 1 is sketched in Figure 5.4. In this model, every beam starts with the zero switch-state s n = 0, where subscript n denotes that the variable is taken at time t n . It will remain in this state until the beam reaches the coercive electric field E c or -E c . This forces the switch-state to change in the next time step to either s n+1 = 1 or s n+1 = -1, denoting positive or negative polarization respectively. After the beam is polarized, the beam can only switch to the opposite polarized state (180 • ) once it reaches the opposite value of the coercive electric field, i.e., -E c for s n = 1 and E c for s n = -1.

The three states of polarization are implemented in this model, although it can be noticed that once a beam is polarized either positive or negative, it cannot return to the zero-state. This first model is quite simple; yet, it allows us to get a good approximation in macroscopic electric displacement D. Although the switch-state changes suddenly, the polarization multiplier is chosen as a smooth function, continuous and differentiable in terms of the hyperbolic tangent. This multiplier depends on the current axial electric field in the beam E 1,n+1 and the switch from the previous time step s n as follows:

s n = 1 s n = 0 s n = -1 E n+1 < -E c E n+1 > E c s n+1 = -1 s n+1 = 1 yes yes s n+1 = s n no no E n+1 < -E c E n+1 > E c
m n+1 (E 1,n+1 , s n ) = 1 -s n 2 tanh a (s n ) E c (E 1,n+1 -E c ) + 3 + 1 + s n 2 tanh a (s n ) E c (E 1,n+1 + E c ) -3 (5.18)
The hyperbolic tangent is a function with odd symmetry and works as a smoothed version of a Heaviside step, with horizontal asymptotic values of -1 and 1 for the former instead of 0 and 1 for the latter. In our case, the chosen form of the argument is a/ |x 0 | (x ∓ |x 0 |) ± 3. The last number has been selected in agreement with tanh (3) =tanh (-3) ≈ 0.995 to translate the function so that when |x| > |x 0 | we can consider that the function has already arrived at the other asymptotic value.

In (5.18), the normalized parameter a (s n ) controls how smooth or how sharp the curve is as can be appreciated in Figure 5.5 left, where the value of s n = -1 has been adopted. The higher the value of a is, the more similar to a Heaviside step the polarization multiplier is. The decision for taking a value of a representative of a crystal is based on when m starts to change from one asymptotic value to the next one. For this model, we adopted a value of a = 200 for s n = ±1 that corresponds to m ≈ ±0.995 when E 1,n+1 = ∓0.97 E c . If s n = 0 then a = 244 to adjust the final part of the curve to the one with polarization as can be seen in Figure 5.5 right.

The polarization multiplier curves for the three switch-states are drawn in Figure 5.6 left, where the typical hysteresis loop for ferroelectric materials can be observed [START_REF] Hwang | Ferroelectric/ferroelastic interactions and a polarization switching model[END_REF]. The choice of parameters for this model makes the behavior of variables affected by polarization very similar to the Heaviside step. Hence, the derivative with respect to the electric field is comparable to the Dirac delta function, as can be seen in the detail of Figure 5.6 right.

Switching model 2

Switching model 2 is capable of representing the depolarization of the beams since the switching criteria include the dependency on mechanical stress. The flowchart for this model is sketched in Figure 5.7.

Every beam starts with the zero switch-state s n = 0. A 90 • switching condition is tested to check if the beam is getting polarized in the main axial directions:

s n+1 = -1 ; if -E 1,n+1 E c + N 1,n N c > 1 s n+1 = 1 ; if E 1,n+1 E c + N 1,n N c > 1 (5.19)
This condition is based on (5.14) combined with the beam hypotheses. Once the beam gets polarized into one of the axial directions, both 180 • and 90 • switchings can happen. The former can be reached only with a significant change in electric variables. As for the latter, only mechanical variables are relevant since electric variables are non-zero in any transversal direction. The switching criterion is the simplified version of (5.14). For s n = 1:

s n+1 = -1 ; if E 1,n+1 < -E c s n+1 = 0 ; if N 1,n < -N c (5.20)
And for s n = -1:

s n = 1 ; if E 1,n+1 > E c s n = 0 ; if N 1,n < -N c (5.21)
One can notice that the criteria, once the beam is polarized in any of the axial directions, are decoupled. I.e., 90 • switching depends only on the mechanic stress, and 180 • switching depends only on the electric field. Observe as well that the axial force is considered at the previous time step to simplify the recursiveness of the calculation of local variables. An iterative method to obtain the current axial force could be developed, but the assumption is that the time step is small enough to consider major differences between previous and current values.

The polarization multiplier for this model is created similarly to the one in the previous model by using the hyperbolic tangent. This time, for every polarization state, there is a different expression:

m n+1 (E 1,n+1 , N 1,n ) = 1 2 1 + tanh a (0) N c (N 1,n + N c ) -3 tanh a (1) E c (E 1,n+1 + E c ) -3 ; if s n = 1 ; m n+1 (E 1,n+1 , N 1,n ) = 1 2 1 + tanh a (0) N c (N 1,n + N c ) -3 tanh a (1) E c (E 1,n+1 -E c ) + 3 ; if s n = -1 ; m n+1 (E 1,n+1 , N 1,n ) = sign (E 1,n+1 ) 2 1 + tanh a (0) |E 1,n+1 | E c + N 1,n N c -1 + 3 ; if s n = 0
(5.22)

In Figure 5.8, all three polarization multiplier situations are drawn, with E 1,n+1 and N 1,n as the independent variables. The change of value for m takes place near the coercive values E c or N c for the axial polarization state with the same smoothened pattern as in the previous yes yes model. As for the depolarized beam, the combined criterion creates a diagonal line of values that nullify the argument of the hyperbolic tangent. Notice that near E 1,n+1 = 0 there is a jump for the polarization multiplier, but it should not affect the stability of the element since once m n+1 ≈ ±1, the polarization state switches to the other corresponding criteria.

s n+1 = 0 s n+1 = s n s n+1 = -1 s n+1 = 1 E 1,n+1 < -E c E 1,n+1 > E c s n = -1 s n = 1 s n = 0 N 1,n < -N c no no yes no E 1,n+1 E c + N 1,n N c > 1 -E 1n+1 E c + N 1,n N c > 1

Mechanical model

Viscoplastic model

The proposed mechanic model will change once the material has reached the yield limit. Here, we propose the one-dimensional viscoplasticity with Armstrong-Frederick non-linear kinematic hardening (e.g. [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]. This model is implemented independently for the axial strain and both shear strains. The description of this constitutive model starts with the hypothesis of additive decomposition of the strains into elastic, plastic, and, for the case of axial strain, remanent parts induced by polarization:

ε 1 = ε e 1 + ε vp 1 + ε r 1 ; γ 12 = γ e 12 + γ vp 12 ;
γ 13 = γ e 13 + γ vp 13 ;

(5.23)

where the bar indicates the regular part of the strains, superscript e denotes the elastic part, v p the viscoplastic part, and r the remanent part. The curvatures κ i are kept elastic; in other words, they are not affected by the viscoplastic or ferroelectric model. Also, the electric charge is split into elastic and remanent (saturated polarization) parts free energy potential ψ is decomposed into three decoupled terms:

Q 1 = Q e 1 + P r 1 . The -N c + → A xi al fo rc e N 1 -E c + → E le c t r ic F ie ld E 1 -1 0 1 Pol. multiplier m -N c + → A xi al fo rc e N 1 E c + → E le c t r ic F ie ld E 1 -1 0 1 Pol. multiplier m N c + → A xi al fo rc e N 1 0 + → E le c t r ic F ie ld E 1 -1 0
ψ 1 ε e 1 , ζ 1 , E 1 = 1 2 ε e 1 Ȇ ε e 1 + 1 2 ζ 1 H is,1 ζ 1 -e 11 ε e 1 E 1 - 1 2 E 1 1 E 1 ; ψ 2 γ e 12 , ζ 2 = 1 2 γ e 12 k c G γ e 12 + 1 2 ζ 2 H is,2 ζ 2 ; ψ 3 γ e 13 , ζ 3 = 1 2 γ e 13 k c G γ e 13 + 1 2 ζ 3 H is,3 ζ 3 (5.24)
The dissipation is defined through the second principle of thermodynamics, and applying the Legendre transformation to the internal energy to exchange the role of Q 1 and E 1 :

D 1 = -ψ1 - 1 A ∂ Q e 1 E 1 ∂t + N 1 A ε1 + Q 1 A Ė1 = N 1 A εvp 1 + εr 1 + q 1 ζ1 + E 1 Ṗ r 1 A ; D 2 = -ψ2 + N 2 A γ12 = N 2 A γvp 12 + q 2 ζ2 ; D 3 = -ψ3 + N 3 A γ13 = N 3 A γvp 13 + q 3 ζ3
(5.25)

In order to get the previous expressions, the following definitions have been established:

N 1 = A ∂ψ 1 ∂ε e 1 ; N 2 = A ∂ψ 2 ∂γ e 12 ; N 3 = A ∂ψ 3 ∂γ e 13 ; Q e 1 = -A ∂ψ 1 ∂E 1 ; q 1 = - ∂ψ 1 ∂ζ 1 ; q 2 = - ∂ψ 2 ∂ζ 2 ; q 3 = - ∂ψ 3 ∂ζ 3 (5.26)
where q i is the stress-like hardening variable in direction i . The viscoplastic part of the dissipation is extracted from (5.25) by removing the remanent part of the strain and the polarization terms because they have different evolution equations as described in the previous section. The yield functions are established as follows:

φ i = |N i -χ i A| -N i,y -q i A (5.27)
where N i,y is the yield limit in direction i , and χ i is the back-stress in direction i , defined for the Armstrong-Frederick model as the implicit differential equation:

χ1 = H l,1 εvp 1 -H n,1 ζ1 χ 1 ; χ2 = H l,2 γvp 12 -H n,2 ζ2 χ 2 ; χ3 = H l,3 γvp 13 -H n,3 ζ3 χ 3
(5.28)

The principle of maximum viscoplastic dissipation can be applied to obtain the evolution of the viscoplastic variables by introducing the corresponding Lagrangian. The main difference with plasticity is that φ i > 0 are admissible but with a penalization proportional to the inverse of the viscosity coefficient η i as follows:

L 1 = -D vp 1 + φ 1 η 1 ; L 2 = -D vp 2 + φ 2 η 2 ; L 3 = -D vp 3 + φ 3 η 3 (5.29)
where the Macaulay brackets are an operator defined as φ = φ if φ > 0 and φ = 0 otherwise, penalizing positive values of the yield function. The evolution of the viscoplastic variables can then be calculated:

∂L 1 ∂N 1 = -εvp 1 + φ 1 η 1 sign (N 1 -χ 1 A) = 0 ⇒ εvp 1 = γ1 sign (N 1 -χ 1 A) ; ∂L 1 ∂q 1 = -A ζ1 + A φ 1 η 1 = 0 ⇒ ζ1 = γ1 ; ∂L 2 ∂N 2 = -γvp 12 + φ 2 η 2 sign (N 2 -χ 2 A) = 0 ⇒ γvp 12 = γ2 sign (N 2 -χ 2 A) ; ∂L 2 ∂q 2 = -A ζ2 + A φ 2 η 2 = 0 ⇒ ζ2 = γ2 ; ∂L 3 ∂N 3 = -γvp 13 + φ 3 η 3 sign (N 3 -χ 3 A) = 0 ⇒ γvp 13 = γ3 sign (N 3 -χ 3 A) ; ∂L 3 ∂q 3 = -A ζ3 + A φ 3 η 3 = 0 ⇒ ζ3 = γ3
(5.30)

where we denoted the viscoplastic multiplier γi = φ i /η i .

Embedded discontinuity for localized failure

A model of embedded discontinuity is triggered once the element reaches the ultimate stress followed by softening response. A jump in displacement is created in the middle of the beam, which constitutes a weak link that undergoes localized failure. Therefore, localized plastic deformation is limited to the discontinuity, whereas other parts of the beam lattice are unloading. The description of this model starts with the inclusion of a jump at x 1 = x 1 in the displacements and voltage:

u i (x 1 , t) = u i (x 1 , t) + u i (t) H x 1 -x 1 V (x 1 , t) = V (x 1 , t) + V (t) H x 1 -x 1 (5.31)
where the bar and the double bar indicate, respectively, the regular part and the jump in displacement in direction i and voltage. With H(x), we denote the Heaviside step function, taking value one for positive argument or zero otherwise. The discontinuity is placed at x 1 = x 1 . Having this result in mind, the total axial and shear strains, and the electric field are calculated as:

ε 1 = ε 1 + ε 1 = ε 1 + u 1 δ x 1 -x 1 ; γ 12 = γ 12 + γ 12 = γ 12 + u 2 δ x 1 -x 1 ; γ 13 = γ 13 + γ 13 = γ 13 + u 3 δ x 1 -x 1 ; E 1 = E 1 + E 1 = E 1 + V δ x 1 -x 1
(5.32) with δ (x) as the Dirac delta function. When softening plasticity is activated, the forces and the electric charge need to be regular for its derivative to comply with the equilibrium with N a as the shape function of node a and a u i a , a ϕ i a or a V a as the nodal values at node a of the displacements, rotations and voltage. These shape functions are the standard interpolating shape functions for one dimensional beam elements:

N 1 (x 1 ) = 1 - x 1 L ; N 2 (x 1 ) = x 1 L (5.38)
where L is the length of the beam and x 1 is measured from node 1 in the axial direction of the beam as depicted for the local frame in Figure 5.2. Thus, coordinate x 1 for node 1 x 1,1 = 0 and for node 2 x 1,2 = L. The derivatives of these shape functions are constant:

B 1 = -1 L ; B 2 = 1 L (5.39)
The interpolation of the crack is simulated with a model of embedded discontinuity with the incompatible mode method. Thus, the influence of the opening does not propagate to neighboring elements. Also, the calculation of the crack opening can be performed in each element. The new shape function M and its derivative G are defined as:

M • • = -N 2 + H x 1 -x 1 ; G • • = G + δ x 1 -x 1 = -B 2 + δ x 1 -x 1 (5.40)
Notice that G denotes the regular part of the derivative, which has to ensure the orthogonality with constant stress mode as described in [START_REF] Ibrahimbegovic | A modified method of incompatible modes[END_REF]. Choosing the middle of the beam as the point where it breaks already guarantees that the shape function is orthogonal. The strain vector and the electric field can be calculated using (5.9) and adding the corresponding incompatible mode:

ε 1 = 2 a=1 B a a u 1 a + Gα 1 ; γ 12 = 2 a=1 B a a u 2 a -N a a ϕ 3 a + Gα 2 ; γ 13 = 2 a=1 B a a u 3 a + N a a ϕ 2 a + Gα 3 ; κ i = 2 a=1 B a a ϕ i a ; E 1 = - 2 a=1 B a a V a (5.41)
For a particular time step t n+1 , denoted with subscript n + 1, the computation is divided into two phases: local and global. Newton-Raphson iterative method, capable of solving the non-linearities present in the formulation, is implemented to make both phases converge. First, in the local stage, all the plastic, softening, and ferroelectric internal variables are made to converge. In the global stage, these converged values are used to verify the equilibrium equations, and to obtain new values of displacements, rotations, and voltage, with the assembly of the stiffness matrix and the force vector. If the equilibrium equation is not verified, the local variables are recalculated with the new values of the degrees of freedom followed by the global phase.

Local computation

At iterative sweep (j + 1), the values of the previous iteration for nodal displacements, rotations, and voltage are known: a u i ,(j) a,n+1 , a ϕ i ,(j) a,n+1 , a V i ,(j) a,n+1 . First, the electric field is calculated in order to obtain the value of the polarization multiplier. The previous stored value of axial force N 1,n is used as well if the switching model 2 is activated.

A trial step, denoted by superscript tr, is assumed where all plastic and softening variables remain unchanged from the converged values in the previous time step:

ε vp,tr 1,n+1 = ε vp 1,n ; γ vp,tr 12,n+1 = γ vp 12,n ; γ vp,tr 13,n+1 = γ vp 13,n ; ζ tr i,n+1 = ζ i,n q tr i,n+1 = q i,n ; α tr i,n+1 = α i,n ; ζ tr i,n+1 = ζ i,n ; q tr i,n+1 = q i,n
(5.42)

The computation is followed by the calculation of the axial strain, modified by polarization material properties in (5.15), and the constitutive terms in (5.16). The trial forces, moments, and the electric charge are computed with the previous values as follows:

N tr 1 = ȆA 2 a=1 B a a u 1 a,n+1 -ε vp 1,n -ε r 1,n+1 + Gα 1,n + e 11 A 2 a=1 B a,n+1 a V a ; N tr 2 = k c GA 2 a=1 B a a u 2 a,n+1 -N a a ϕ 3 a -γ vp 12,n + Gα 2,n ; N tr 3 = k c GA 2 a=1 B a a u 3 a,n+1 + N a a ϕ 2 a -γ vp 13,n + Gα 3,n ; M tr 1 = GJ 2 a=1 B a a ϕ 1 a,n+1 ; M tr 2 = ȆI 2 a=1 B a a ϕ 2 a,n+1 ; M tr 3 = ȆI 2 a=1 B a a ϕ 3 a,n+1 ; Q tr 1 = e 11 A 2 a=1 B a a u 1 a,n+1 -ε vp 1,n -ε r 1,n+1 + Gα 1,n -1 A 2 a=1 B a a V a,n+1 (5.43)
The path to follow next is sketched in the diagram of Figure 5.9. If the element has not broken yet, a verification on if the yield limit N i,y being lower or not than the fracture limit N i,f is made. Plasticity can only occur if this verification holds; otherwise, the material is considered brittle enough so that elastic regime is followed directly by softening.

For the brittle case, the softening criteria is verified with trial values:

φ tr i,n+1 = Ntr i,n+1 -N i,f -q tr i,n+1 A (5.44)
If φ tr i,n+1 ≤ 0 the trial step is accepted and all local variables are guessed correctly. Otherwise, the crack switch-state is turned on for future time steps and softening regime activates.

For the case where the material can plastify, the yield criteria is verified with trial values: φ tr i,n+1 = N tr i,n+1χ tr i,n+1 A -N i,yq tr i,n+1 A (5.45) If φ tr i,n+1 ≤ 0, the hypothesis of elastic regime is correct and all local variables remain the same as in the previous time step. Otherwise, there is still uncertainty on wether plasticity or softening have to be active. To resolve this indecision, the softening criteria is tested with plastic variables: (5.46) where the plastic force at the discontinuity is defined as in (5.43) but with the updated value of viscoplastic strain. If φ pl i,n+1 ≤ 0, the viscoplastic step is accepted; otherwise, the viscoplastic variables are reset to the values of the previous step and the softening variables are calculated.

φ pl i,n+1 = Npl i,n+1 -N i,f -q tr i,n+1 A
If the previous crack switch-state indicates that the beam has already been broken, the only regime to check is softening, as in the brittle case. Therefore, the same procedure is followed.

Viscoplastic regime

After checking that φ tr i,n+1 > 0 and φ pl i,n+1 ≤ 0, the evolution of the viscoplastic variables described in (5.30) is discretized for implicit time integration scheme. The following development is an extension to the one developed in [START_REF] Ibrahimbegovic | Classical plasticity and viscoplasticity models reformulated: theoretical basis and numerical implementation[END_REF]. This process is analogous for the other directions, taking into account the role equivalence of the variables for every direction. For direction 1:

ε vp 1,n+1 = ε vp 1,n + γ 1,n+1 sign (N 1,n+1 -χ 1,n+1 A) ; ζ 1,n+1 = ζ 1,n + γ 1,n+1 ; q 1,n+1 = -H is ζ 1,n+1 ;
χ 1,n+1 = χ 1,n + H l,1 γ 1,n+1 sign (N 1,n+1χ 1,n+1 A) -H n,1 γ 1,n+1 χ 1,n+1 ; N 1,n+1 = N tr 1,n+1 -ȆA γ 1,n+1 sign (N 1,n+1χ 1,n+1 A)

(5.47)

Substracting chi 1,n+1 multiplied by the area to N 1,n+1 , introducing the null term γ 1 H n,1 -N 1,n+1 + N tr 1,n+1 -ȆA γ 1 and simplyfying:

(N 1,n+1χ 1,n+1 A) = 1 + γ 1,n+1 H n,1 N tr 1,n+1χ 1,n Aγ 1,n+1 ȆA + H l,1 A sign (N 1,n+1χ 1,n+1 A) + H n,1 (N 1,n+1χ 1,n+1 A) + H n,1 ȆAγ 1,n+1 sign (N 1,n+1χ 1,n+1 A)

(5.48) By applying the identity a = |a| sign (a) to the previous equation, we obtain:

|N 1,n+1 -χ 1,n+1 A| + γ 1,n+1 ȆA + H l,1 A + H n,1 |N 1,n+1 -χ 1,n+1 A| + H n,1 ȆAγ 1,n+1 sign (N 1,n+1 -χ 1,n+1 A) = 1 + γ 1,n+1 H n,1 N tr 1,n+1 -χ 1,n A sign 1 + γ 1,n+1 H n,1 N tr 1,n+1 -χ 1,n A
(5.49)

By taking into account that all constitutive coefficients are positive, the following relation with the sign operator holds:

sign (N 1,n+1χ 1,n+1 A) = sign 1 + γ 1,n+1 H n,1 N tr 1,n+1χ 1,n A (5.50) Therefore, the terms multiplying both signs in (5.49) must be equal. Introducing the definition of viscoplastic multiplier, the yield criteria enters the equation:

φ i,n+1
= |N 1,n+1χ 1,n+1 A| -N i,yq 1,n+1 = γ 1,n+1 η i ∆t (5.51)

where ∆t is the current time step. The norm |N 1,n+1χ 1,n+1 A| in (5.49) can be replaced now. Regroupping:

1 + γ 1,n+1 H n,1 N tr 1,n+1 -χ 1,n A -N 1,y + H is A ζ 1,n -γ 1,n+1 ȆA + H l,1 A +H is A + η 1 ∆t + H n,1 N 1,y + H is A ζ 1,n + γ 1,n+1 H n,1 η 1 ∆t + ȆA + H is A = 0 (5.52)
The Newton-Raphson iterative method is used to solve γ 1,n+1 in the previous quadratic equation. This method requires a good guess of the initial value so that the positive value is found. A good guess can be found in [START_REF] Hadzalic | Theoretical formulation and seamless discrete approximation for localized failure of saturated poro-plastic structure interacting with reservoir[END_REF], where the following expression is obtained for explicit back-stress Armstrong-Frederick:

γ (0) 1,n+1 = φ tr i,n+1 ȆA + H is A + H l,1 A + η i ∆t -H n,1 Aχ 1,n
(5.53)

Softening regime

When softening regime is activated, i.e. with φ tr i,n+1 > 0 for the brittle and already broken case or with φ pl i,n+1 > 0 otherwise, the crack switch-state in direction i is activated: s i,n+1 = 1. The development is further described in detail only for direction 1, given that the same approach can be extrapolated to any of the other directions. First, the evolution of softening variables and forces is described as:

α 1,n+1 = α 1,n + γ 1,n+1 sign Ni,n+1 ζ 1,n+1 = ζ 1,n + γ 1,n+1 q 1,n+1 = N 1,f A 1 -exp -ζ 1,n+1 N 1,f A G 1,f N1,n+1 = Ntr 1,n+1 + ȆA G (α 1,n+1 -α 1,n ) = Ntr 1,n+1 -ȆA B 2 γ 1,n+1 sign N1,n+1
(5.54)

Manipulating the last equation to introduce the norm of the force at the discontinuity:

Ni,n+1 + ȆA B 2 γ 1,n+1 sign Ni,n+1 = Ntr i,n+1 sign Ntr i,n+1

(5.55)

Considering that all terms multiplying the sign operators are always positive, the following relations hold: sign Ni,n+1 = sign Ntr i,n+1

Ni,n+1 + ȆA B 2 γ 1,n+1 = Ntr i,n+1

(5.56)

Therefore, trial values can be used to eliminate recursiveness. For plastic softening, φ i,n+1 has to be forced to zero. Thus:

φ i,n+1 = φ tr i,n+1
-ȆA B 2 γ 1,n+1 + q 1,n+1q 1,n A = 0 (5.57)

Introducing explicitly the expression for q 1,n+1 and q 1,n and simplifying:

φ tr i,n+1 -ȆA B 2 γ 1,n+1 + N 1,f exp -ζ 1,n N 1,f A G 1,f -exp -ζ 1,n -γ 1,n+1 N 1,f A G 1,f = 0
(5.58)

To obtain the value of γ 1,n+1 , an iterative method is applied as in the previous example:

φ (k+1) i,n+1 = φ (k) i,n+1 + ∂φ (k) i,n+1 ∂γ (k) 1,n+1 ∆γ (k)
1,n+1 = 0 (5.59)

A good initial guess is γ (0)

1,n+1 = 0. The method stops whenever the desired convergence is reached.

Global computation

The starting point for the global implementation is provided by the weak form of the equilibrium equations (5.11), where the approximations (5.37) and (5.41) have been introduced: ˆL 0 B a w u i a N i,n+1 dx 1w u i a Ni Γ N = 0 ;

ˆL 0 B a w ϕ i a M i,n+1 dx 1w ϕ i a Mi Γ M = 0 ;

ˆL 0 B a w V a Q 1,n+1 dx 1 -w V a Q1 Γ Q = 0
(5.60) with w u i a , w ϕ i a , and w V i a as the virtual nodal values. These nodal values are zero if a Dirichlet boundary condition is applied to node a. Otherwise, the residual, which is obtained by isolating the virtual degrees of freedom, has to be null:

R u i a,n+1 = ˆL 0 B a N i,n+1 dx 1 -Ni Γ N = 0 ; R ϕ i a,n+1 = ˆL 0 B a M i,n+1 dx 1 -Mi Γ M = 0 ; R V a,n+1 = ˆL 0 B a Q 1,n+1 dx 1 -Q1 Γ Q = 0
(5.61)

where the forces, moments and electric charge are defined as:

N 1 = ȆA 2 a=1 B a a u 1 a,n+1 -ε vp 1,n+1 -ε r 1,n+1 + Gα 1,n+1 + e 11 A 2 a=1
B a a V a,n+1 ;

N 2 = k c GA 2 a=1
B a a u 2 a,n+1 -N a a ϕ 3 a,n+1γ vp 12,n+1 + Gα 2,n+1 ;

N 3 = k c GA 2 a=1
B a a u 3 a,n+1 + N a a ϕ 2 a,n+1γ vp 13,n+1 + Gα 3,n+1 ;

M 1 = GJ 2 a=1
B a a ϕ 1 a,n+1 ;

M 2 = ȆI 2 a=1
B a a ϕ 2 a,n+1 ;

M 3 = ȆI 2 a=1
B a a ϕ 3 a,n+1 ;

Q 1 = e 11 A 2 a=1 B a a u 1 a,n+1 -ε vp 1,n+1 -ε r 1,n+1 + Gα 1,n+1 -1 A 2 a=1 B a a V a,n+1
(5.62)

In order to solve the nonlinearities of the previous equations, Newton-Raphson method is used: R (5.63)

At each iterative sweep, we can then perform the corresponding state variable updates according to: a (5.64)

In the first iteration within each time step, the starting guess is assumed to be the converged value at the previous step: a The derivative term of (5.63) is the stiffness matrix, constructed differently, depending on the local regime the bar is. words, when the beam is close to the change of polarization. This model may seem trivial but it is capable of capturing in a lattice model the behavior of a heterogeneous material, as we can see in the next example.

Macro-model subjected to cyclic electric field applied in vertical direction

The geometry of the numerical example in this section is a cube, with sides of 20 cm as can be seen in Figure 5.11. On the left, the Voronoi representation of the cube is pictured, where the different colors show the boundary faces of the Voronoi cells. On the right, the corresponding mesh of lattice beams to solve. The material properties are extracted from [START_REF] Hwang | Ferroelectric/ferroelastic interactions and a polarization switching model[END_REF], and given in Table 5.2. respectively. The voltage at the top is triangular, starting with a value of 0, increasing to V max = 0.2 MV, then decreasing to -V max , and finally returning to V max . In planes x 1 = 0, x 2 = 0, and x 3 = 0, the corresponding perpendicular displacement is prescribed to simulate symmetry boundary conditions for one eighth of the specimen. The Neumann boundary condition is also set introducing a constant compressive force. Rotation degrees of freedom are left free. All switches are set initially to zero.

The electric and mechanic variables calculated with the beam model have to be expressed in the global frame and averaged with the volume to interpret the corresponding macro response of the material by using: Figure 5.12: Averaged vertical electric displacement and strain obtained for the numerical example

In Figure 5.12, the averaged response for the numerical example is plotted, where the typical hysteresis loops for electric displacement and strain can be observed. The curves in black are without the compressive force whereas the blue ones have that boundary condition. The difference lies in the amount of beams that get polarized. The compressive force makes it more difficult to get polarized as depicted in (5.19).

The averaged strain is affected as well by the introduction of the compressive force. The strain is almost an offset of the one without the force. The shape of the butterfly loop is not as perfect as in the single-crystal simulation. Namely, the minimum stress is not close to zero, and the spontaneous polarization is a little bit higher than the expected value, according to [START_REF] Hwang | Ferroelectric/ferroelastic interactions and a polarization switching model[END_REF]. We can see here one of the limitations of this method is that cannot capture transversal effects with accuracy, but the axial variables are close to the ones of the continuum model.

Conclusions

In this chapter, some concluding remarks from the previous chapters are drawn, highlighting the contributions of this thesis along with the possible future extensions.

Conclusions

Several novelties for solving the coupled thermo-electro-magneto-mechanic problems for conductor and insulators materials are presented in this work. In particular, a complete formulation for thermo-electro-magneto-elastoplastic behavior has been developed consistently by using conservation principles along with the definition of a free-energy potential leading to the corresponding constitutive relations, providing a hyper-elastic response (rather than a hypoelastic) with the full coupling between all state variables and their duals. This formulation has been implemented in the most general 3D framework, by using an 8-node hexahedral finite element for constructing semi-discretization, along with the global phase of time discretization by the Newmark scheme and local computation of the internal variables for plasticity.

This development provides a sound basis to analyze practical examples, which either did not include any of the main interactions to the fully-coupled formulation or have not implemented a plastic model in a thermo-electro-magneto-elastic framework.

The chosen numerical simulations have tested the model capabilities to represent firstorder couplings in transient and static cases. The results illustrate the wealth of coupled responses simulations, including the transient fully-coupled pulse propagation in a 3D setting, or the more industry-oriented Peltier cells.

This thesis also provides a first step in establishing the firm link between electromagnetic and mechanical point of view, introducing several novelties on the way for solving numerically electrostatic and magnetostatic problems. In particular, energy-based variational formulations for several single-field and mixed multi-field formulation have been developed instead of the constitutive-based global equations currently used cell method, the preferred tool of physics and the electrical engineering community. This energy formulation allows studying induced heterogeneities typical of inelastic material behavior and not only linear constitutive laws typical of classical works on electromagnetics. Although a full development in this thesis is mostly provided for electrostatics, it is shown as well that an equivalent magnetostatic formulation can be recovered with an adequate change of the corresponding field variables.

The numerical implementation of discrete approximation is based upon Whitney's elements. In this thesis, we have developed in detail the approach granting the partition-of-unity to the chosen discrete approximation in the tetrahedron element, which applies to the corresponding shape functions. Moreover, the proposed discrete approximation improves the accuracy of the electromagnetic fields, since it preserves either tangential or normal continuity of primal and dual variables, respectively. This method can also be extended to higher-order interpolation.

The numerical examples in Chapter 4 provide full validation of the proposed discrete approximations based upon tetrahedra elements. They also illustrate the full scope of the problems and possible applications that can be solved with this kind of element, first for the case of electrostatics and second by analogy for magnetostatics. Of particular interest is the superior result is achieved in the parallel plate capacitor simulations establishing accuracy even with a small number of elements for here proposed Hellinger-Reissner formulation in comparison with any other approach.

Chapter 5 presents a theoretical micro-model formulation that couples ferroelectricity, viscoplasticity, and localized failure within the same element in a consistent manner. Each of the elements simulates a single crystal of the whole macrostructure that can be polarized in the axial or transversal direction. This feature makes the element suitable for multi-scale computations since it can replicate the behavior of the macro response of ferroelectric materials by introducing the properties of a single crystal.

The introduction of an embedded discontinuity model could be crucial in applications under significant or cyclic loadings of any kind, not only mechanical. Potential current and future applications of this type include microsensors that measure outside their usual range, overcharged actuators, charging solid-state batteries that undergo large temperature gradients.

The numerical examples presented in this chapter show a good agreement between experiments and the response of the model. In some cases, the polarization switching model 1 could be enough to capture the axial response accurately, although model 2 is not much more expensive in terms of computation resources.

Future work

The work presented in this thesis establishes the basis for future work on the research of numerical methods on smart materials and further refinement on the models introduced here.

For instance, the macroscopic model in Chapters 2 and 3 can be improved with more accurate 3D plastic models analogous to the ones in Chapter 5. Also, the introduction of an embedded discontinuity model would provide the behavior past the ultimate stress. Regarding Whitney's discretizations, this thesis presents an implementation for finite elements that can be completed by combining mechanical degrees of freedom in the standard finite element fashion. Also, we can still adapt the discretizations of other fields to Whitney's to obtain more accuracy.

Finally, the microscopic model introduced in Chapter 5 can be used to test more interesting and realistic applications that involve different fields and possible crack openings. For example, microsensors in a structure, or solid-state batteries under extreme thermal gradients.

One of the struggles of the formulations in this thesis is the characterization of the parameters that coupling introduces. Experiments with several fields are essential to check the validity of the results in this thesis, the hypotheses assumed, and the numerical importance of plastic parameters.

2. 1

 1 Different behavior of conductors (left) and dielectric (right) electric materials. 14 2.2 Sketch of the boundary conditions considered for the actuator validation example. Three planes of symmetry considered in x 1 = x 2 = x 3 = 0; electric potential prescriptions at the top and bottom planes. Measures in mm. . . . 27 2.3 Displacement in directions x 1 and x 3 and electric potential distributions for the actuator validation piezoelectric case. . . . . . . . . . . . . . . . . . . . . . 27 2.4 Displacement in directions x 1 and x 3 and electric potential distributions for the generator validation piezoelectric case. . . . . . . . . . . . . . . . . . . . . . 28 2.5 Sketch of the boundary conditions for the third piezoelectric example. Two planes of symmetry considered in x 1 = x 2 = 0; electric potential prescriptions at the top and bottom planes; the latter face clamped. Measures in mm. . . 29 2.6 Displacement in directions x 1 and x 3 and deformed configuration amplified 10 5 . 30 2.7 Displacement in directions x 1 and x 3 and deformed configuration amplified 10 5 . 31 2.8 Displacement in directions x 1 and x 3 and temperature . . . . . . . . . . . . . 31 2.9 Left, transient temperature distribution at the top of the geometry of a thermoelastic case; right, temperature distribution along the vertical direction for time slightly lower than t = 1, 2 s. . . . . . . . . . . . . . . . . . . . . . . . 32 2.10 Scheme of the geometry considered and the pulse applied to the left end of the bar. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.11 Computed displacement, electric potential, magnetic potential and temperature in the bar middle point x 3 = l/2. . . . . . . . . . . . . . . . . . . . . . 33 2.12 Contour plot for electric and magnetic potential at t = 1.7 × 10 -5 s. . . . . 34 2.13 Computed displacement, voltage and temperature in the bar middle point x 3 = l/2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 v 3.1 Left, mesh used in this numerical example, each material represented in a different color: Bi 2 Te 3 in red, copper in blue, solder in green and alumina in yellow. Right, boundary conditions imposed in this problem. . . . . . . . . . 49 3.2 Maximum von Mises stress for solder (red), copper (black) and Bi 2 Te3 (blue) when an electric square pulse of gain P = 3.5 and duration ∆t p = 5 s is applied at t = 125 s. Full line, computed result with elastoplastic response; dashed line, computed results for elastic response. . . . . . . . . . . . . . . . . . . . 50 3.3 Von Mises stress [MPa] contour in deformed configuration for elastic (left) and plastic regime (right). All materials but the alumina. Detailed zoom at the bottom side of the Peltier cell. . . . . . . . . . . . . . . . . . . . . . . . 51 3.4 Top: electric flux pulse gain introduced in the problem. Bottom: temperature at the cold side for elastic and plastic regime. . . . . . . . . . . . . . . . . . 52 4.1 Tonti diagram for electrostatic problem with scalar potential. . . . . . . . . . 59 4.2 Locations of unknowns at nodes, edges and facets in a tetrahedral element and the corresponding positive orientations, with the reference frame placed at node 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.3 Schema of the location of the nodes depending on their location on a regular 15-node isoparametric tetrahedron. . . . . . . . . . . . . . . . . . . . . . . . 62 4.4 Tree set (continuous lines) for regular mesh. Every hexahedron represents six tetrahedra in cubic disposition. Every point is a vertex node. . . . . . . . . . 66 4.5 Voltage and electric displacement distributions calculated for the first validation example. First row, V-formulation; second row, Hellinger-Reissner. . . . 69 4.6 Voltage and electric displacement distributions calculated for the second validation example. First row, V-formulation; second row, Hellinger-Reissner. . . 70 4.7 Voltage and electric displacement distributions calculated for the third validation example. First row, V-formulation; second row, Hellinger-Reissner. . . . 71 4.8 Energy convergence for both V-formulation and Hellinger-Reissner, and the exact solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.9 Schematic geometry representation of a quarter of the hollowed cylinder. Symmetry planes in x = 0 and y = 0. Dimensions in cm. . . . . . . . . . . . . . 72 4.10 Voltage and norm of electric displacement distributions calculated for the fourth validation example, with both lateral faces with V boundary condition. First row, V-formulation; second row, Hellinger-Reissner; third row, complementary energy formulation. . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.11 Left, energy convergence for V-formulation, Hellinger-Reissner and Complementary energy formulation compared with the exact solution. Right, computation time for all the formulations tested. . . . . . . . . . . . . . . . . . . . 75 vi 4.12 Schematic geometry representation of a quarter of the parallel plate capacitor; two charged conductor plates and a dielectric cylinder. Symmetry planes in x = 0 and y = 0. Dimensions in cm. . . . . . . . . . . . . . . . . . . . . . . 4.13 Electric field and displacement vector representation calculated for the parallel plate capacitor example. First row, V-formulation; second row, Hellinger-Reissner. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Total electrostatic energy by number of elements of the mesh with two formulations: V-formulation (V-form) and Hellinger-Reissner (H-R). . . . . . . 5.1 Left, two adjacent Voronoi cells held together by the corresponding lattice model. Right, the assumed shape of the beam and the approximation of the circular cross-section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Local frame representation for a beam. The local direction 1 is oriented from local node 1 to 2. Local directions 2 and 3 are randomly chosen, bu remain perpendicular to the beam axis. . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Sketch of a tetragonal unit cell for PbTiO 3 below (tetragonal) and above (cubic) Curie temperature. Pb in black, O in white and Ti in gray. . . . . . . 5.4 Flowchart to select the switch-state for the next time step s n+1 based on its previous value s n and the current value of electric field E 1,n+1 . . . . . . . . 5.5 Left, smoothed polarization multiplier m curves for negative switch-state s n = -1 with different values of parameter a = {10, 25, 60, 160, 400, 1000}. Right, detail of curves with chosen value a = 200 for |s n | = 1 and a = 244 for s n = 0. 5.6 Left, smoothed polarization multiplier m curves for all three previous state of polarization s n with the chosen a = 200. Right, detail of the corresponding derivatives of m with respect to the electric field. . . . . . . . . . . . . . . . 5.7 Flowchart to select the switch for the next time step based on the previous state of switch s n , the current value of electric field E n+1 , and the previous value of axial force N 1,n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Smoothed polarization multiplier m curves for switching model 2 with different values of the previous switching state s n = 1 for top left, s n = -1 for top right, and s n = 0 for bottom figure. . . . . . . . . . . . . . . . . . . . . . . 5.9 Diagram of decision for the current local regime depending on the previous crack switch-state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Electric charge and strain obtained for the beam of the numerical example. . 5.11 Mesh used for the numerical examples. Left, Voronoi cell representation; right, lattice model, dual to Voronoi . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 Averaged vertical electric displacement and strain obtained for the numerical example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables 2.1 List of the state variables and their corresponding dual variables for coupled thermoelasticity and electromagnetism . . . . . . . . . . . . . . . . . . . . . 18 2.2 Properties of the material used in validation example . . . . . . . . . . . . . 35 3.1 Material properties for every material used in the Peltier Cell . . . . . . . . . 49 5.1 Material properties for the material used in the single crystal polarization simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.2 Material properties for the material used in the macro-model. . . . . . . to a boundary surface (-) b -Body force vector (N/m 3 ) t -Traction vector (Pa/m2) e -Internal energy density (J/m 3 ) s

  (C/m 2 K) π h -Pyromagnetic tensor (kg/s 2 A K) υ -Magnetoelectric tensor (s/m)

  Figure 2.1: Different behavior of conductors (left) and dielectric (right) electric materials.

Figure 2

 2 Figure 2.7: Displacement in directions x 1 and x 3 and deformed configuration amplified 10 5 .

  Figure 2.8: Displacement in directions x 1 and x 3 and temperature .

Figure 2

 2 Figure 2.9: Left, transient temperature distribution at the top of the geometry of a thermoelastic case; right, temperature distribution along the vertical direction for time slightly lower than t = 1, 2 s.

Figure 2

 2 Figure 2.10: Scheme of the geometry considered and the pulse applied to the left end of the bar.

  Figure 2.11: Computed displacement, electric potential, magnetic potential and temperature in the bar middle point x 3 = l/2.

  Figure 2.13: Computed displacement, voltage and temperature in the bar middle point x 3 = l/2.

Figure 3

 3 Figure 3.1: Left, mesh used in this numerical example, each material represented in a different color: Bi 2 Te 3 in red, copper in blue, solder in green and alumina in yellow. Right, boundary conditions imposed in this problem.

  Figure3.2: Maximum von Mises stress for solder (red), copper (black) and Bi 2 Te 3 (blue) when an electric square pulse of gain P = 3.5 and duration ∆t p = 5 s is applied at t = 125 s. Full line, computed result with elastoplastic response; dashed line, computed results for elastic response.

Figure 3 . 3 :

 33 Figure 3.3: Von Mises stress [MPa] contour in deformed configuration for elastic (left) and plastic regime (right). All materials but the alumina. Detailed zoom at the bottom side of the Peltier cell.

Figure 4 . 1 :

 41 Figure 4.1: Tonti diagram for electrostatic problem with scalar potential.

  Figure 4.2: Locations of unknowns at nodes, edges and facets in a tetrahedral element and the corresponding positive orientations, with the reference frame placed at node 1.

  Figure 4.3: Schema of the location of the nodes depending on their location on a regular 15-node isoparametric tetrahedron.

  Figure 4.4: Tree set (continuous lines) for regular mesh. Every hexahedron represents six tetrahedra in cubic disposition. Every point is a vertex node.

Figure 4 . 5 :

 45 Figure 4.5: Voltage and electric displacement distributions calculated for the first validation example. First row, V-formulation; second row, Hellinger-Reissner.

Figure 4 . 6 :

 46 Figure 4.6: Voltage and electric displacement distributions calculated for the second validation example. First row, V-formulation; second row, Hellinger-Reissner.

  Figure 4.7: Voltage and electric displacement distributions calculated for the third validation example. First row, V-formulation; second row, Hellinger-Reissner.

  Figure 4.8: Energy convergence for both V-formulation and Hellinger-Reissner, and the exact solution.

  Figure 4.10: Voltage and norm of electric displacement distributions calculated for the fourth validation example, with both lateral faces with V boundary condition. First row, V-formulation; second row, Hellinger-Reissner; third row, complementary energy formulation.

Figure 4 .

 4 Figure 4.11: Left, energy convergence for V-formulation, Hellinger-Reissner and Complementary energy formulation compared with the exact solution. Right, computation time for all the formulations tested.

  Figure 4.14: Total electrostatic energy by number of elements of the mesh with two formulations: V-formulation (V-form) and Hellinger-Reissner (H-R).

Figure 5

 5 Figure 5.2: Local frame representation for a beam. The local direction 1 is oriented from local node 1 to 2. Local directions 2 and 3 are randomly chosen, bu remain perpendicular to the beam axis.

Figure 5 . 3 :

 53 Figure 5.3: Sketch of a tetragonal unit cell for PbTiO 3 below (tetragonal) and above (cubic) Curie temperature. Pb in black, O in white and Ti in gray.

  Figure 5.4: Flowchart to select the switch-state for the next time step s n+1 based on its previous value s n and the current value of electric field E 1,n+1

  Figure 5.5: Left, smoothed polarization multiplier m curves for negative switch-state s n = -1 with different values of parameter a = {10, 25, 60, 160, 400, 1000}. Right, detail of curves with chosen value a = 200 for |s n | = 1 and a = 244 for s n = 0.

  Figure 5.6: Left, smoothed polarization multiplier m curves for all three previous state of polarization s n with the chosen a = 200. Right, detail of the corresponding derivatives of m with respect to the electric field.

Figure 5

 5 Figure 5.7: Flowchart to select the switch for the next time step based on the previous state of switch s n , the current value of electric field E n+1 , and the previous value of axial force N 1,n .

1

  Figure 5.8: Smoothed polarization multiplier m curves for switching model 2 with different values of the previous switching state s n = 1 for top left, s n = -1 for top right, and s n = 0 for bottom figure.

Figure 5

 5 Figure 5.9: Diagram of decision for the current local regime depending on the previous crack switch-state.

  Permittivity 1 C/V m 56.3 × 10 -9 Spontaneous polarization P s C/m 2 0.25 Piezoelectric term e 11 C/m 2 1.18 × 10 -9 Coercive Electric Field E c V/m 0.36 × 10 6Table5.2: Material properties for the material used in the macro-model.

Figure 5 .

 5 Figure 5.11: Mesh used for the numerical examples. Left, Voronoi cell representation; right, lattice model, dual to Voronoi The top and bottom faces have prescriptions of the voltage of V = 0 and V = V t (t),

Table 3

 3

		Units	Alumina Copper Tin-Lead Bi 2 Te 3
	Hardening modulus K	GPa	5	5	5	1
	Yield stress σ y	MPa	15,400	100	180	55
	Mass density ρ	kg/m 3	3,570	8,960	7,310	7,530
	Specific heat c p	J/kg K	837	385	226	544
	Expansion coeff. α T First Lamé parameter λ	10 -6 / • C GPa	5.0 163.0	17.0 71.6	27.0 32.5	16.8 67.1
	Shear stress µ	GPa	151.0	43.9	16.8	16.8
	Thermal cond. κ	W/K m	35.3	386.0	48.0	κ(T )
	Electric cond. γ 10 6 A/V m	0	58.1	4.72	γ(T )
	Seebeck coeff. α	V/K	0	0	0	α(T )

.1: Material properties for every material used in the Peltier Cell

The circumcenter is defined as the point that is located at the same distance of each vertex of the polyhedron. It is the center of the circumscribed sphere.

compared against the Joule dissipative term or Peltier cooling. Thus, it would appear that if the main goal in this problem concerned the temperature evolution, the plasticity would not have decisive role for computing the sufficiently good results.

The boundary conditions for the equivalent previous numerical examples are:

Thus, the voltage expressions for each of the three cases are:

All these cases show good agreement with the analytical solution. As can be seen in Fig. 4.10 left, where the case for Laplace equation and imposed V at the external face is plotted, voltage distribution is no longer linear, but rather follows a logarithmic distribution starting from zero at the internal lateral face to 20 at the external lateral face. The problem being axisymmetric results with the only variation in the radial direction.

In the Fig. right, the norm of electric displacement is plotted. This variable follows an inverse distribution with the radius. In the figure, different plot bounds can be observed for both formulations as the exact solution cannot be reached with linear tetrahedral elements. The exact values for |D(r

The complementary energy formulation has also been used to solve the same problem as a comparison with the two formulations used through this section. As previously discussed, the voltage cannot be obtained directly from the calculated variables. On the other hand, the norm of electric displacement is plotted, obtaining similar or even more accurate results as with the previous formulations.

In Fig. 4.11 left, a small mesh convergence study for the case run above is shown. In this study, all formulations converge to the analytical solution using a few thousands of elements, being the complementary energy formulation the more accurate for every mesh configuration and V-formulation the more inaccurate.

However, as can be seen in the Fig right, the latter is unquestionably the fastest of them all. The slowest, in any case, is Hellinger-Reissner due to the higher amount of variables, which also make it the one that spends the most memory. Complementary energy shoots up as well due to the necessity to solve the Dirichlet boundary conditions, which makes it inefficient for being a single field formulation. equations.

With this hypothesis in mind, all the jump terms have to be cancelled. Thus, the rate of crack openings α i and the new hardening variable ζ i are functions of the plastic multipliers:

The softening criteria φ in every direction is defined at the discontinuity:

where Ni represents the force in direction i at the discontinuity x 1 = x 1 , N i,f is the ultimate force in direction i , and q i is the stress-like variable for exponential softening:

with G i,f as the fracture energy. This model of embedded discontinuity includes a crack switch-state s i for each direction i to describe if the element has not cracked yet s i = 0, or it has already reached the failure stress s i = 1.

Numerical implementation

In this section we present the details of the discrete approximation constructed by the finite element method. First of all, the displacements, rotations, and voltage are the variables acting as the degrees of freedom. They are approximated along the beam element using the standard linear interpolating functions and the value of the fields at the nodes:

(5.37)

Numerical examples

In this section, we present the results of several illustrative numerical simulations.

Single-Crystal simulation

In this example, a single beam of length L = 0.1 m is simulated to obtain the characteristic hysteresis loops. This beam is clamped on the left side and subjected to ground voltage. On the right side, a variable voltage is set in terms of a triangular force, with maximum and minimum values of 300 and -300 V. The material properties are gathered in Table 5.1

Table 5.1: Material properties for the material used in the single crystal polarization simulation

The results are displayed in Figure 5.10. On the left, the electric charge Q 1 has the characteristic loop for single crystal. The slope does not change since the non-polarized and the polarized value of the electric permittivity are the same. Also, the value of Q 1 at E = 0 is the value of P s A. On the right, the strain is plotted obtaining the typical butterfly loop. Notice that the strain takes negative values when the electric field is close to the coercive value; in other