Constantinos Kitsos

Gildas Besançon

Christophe Prieur

High-gain observer for 3 × 3 linear heterodirectional hyperbolic systems

Keywords: linear hyperbolic systems, high-gain observers, C 1 -convergence

The problem of High-Gain Observer Design is addressed for a class of 3 × 3 inhomogeneous linear hyperbolic systems with possibly distinct characteristic velocities and considering distributed measurement of part of the state. Applying an infinitedimensional state transformation, the system is mapped into a new set of partial differential equations, satisfying an appropriate form for observer design. The observer includes output correction terms and also spatial derivatives of the output, and ensures arbitrarily fast state estimation. The result is illustrated with a simulated example.

Introduction

High-gain observer design for nonlinear finite-dimensional systems has been considerably developed (see e.g. [START_REF] Khalil | High-gain observers in nonlinear feedback control[END_REF] and references therein). In short, it relies on a single tuning coefficient, to be chosen large enough so as to ensure exponential -and arbitrarily fast -convergence. High-gain observers apply to a large class of cases corresponding to uniformly observable systems Gauthier and Bornard (1981). In the recent papers [START_REF] Kitsos | High-gain observer design for a class of hyperbolic systems of balance laws[END_REF], [START_REF] Kitsos | High-gain observer design for a class of integro-differential quasilinear hyperbolic systems -Application to an epidemic model[END_REF], this approach was extended to infinite-dimensional systems, namely, quasi-linear hyperbolic systems made of n × n balance laws with one characteristic velocity and considering distributed measurement of a part of the state. In [START_REF] Kitsos | A high-gain observer for a class of 2x2 hyperbolic systems with C1 exponential convergence[END_REF] the more general case of two distinct characteristic velocities was considered, but for 2 × 2 systems written in a triangular form.

Other observer designs for hyperbolic systems can be found in the literature, but mainly considering the full state vector on the boundaries as measurement, and

The work of the third author has been partially supported by MIAI@Grenoble Alpes (ANR-19-P3IA-0003).

Email addresses: constantinos.kitsos@laas.fr (Constantinos Kitsos), gildas.besancon@grenoble-inp.fr (Gildas Besançon), christophe.prieur@grenoble-inp.fr (Christophe Prieur).

without the high gain features. Amongst others, one can mention Di Meglio et al (2013), [START_REF] Hasan | Boundary observer design for hyperbolic PDE-ODE cascade systems[END_REF], and [START_REF] Anfinsen | An adaptive observer design for n + 1 coupled linear hyperbolic PDEs based on swapping[END_REF] for backstepping designs, Besançon et al (2006) for a discretization approach (earlylumping techniques), Castillo et al (2013) for direct infinite-dimension-based Lyapunov techniques (see also Besançon et al (2013)), or [START_REF] Nguyen | State and parameter estimation in 1-d hyperbolic PDEs based on an adjoint method[END_REF] for optimization methods. For semigroup and spectral-based methods, see for instance Curtain (1982), Demetriou (2004) and [START_REF] Schaum | A simple observer scheme for a class of 1-D semilinear parabolic distributed parameter systems[END_REF]. In Christophides and Daoutidis (1996), there appear distributed observers relying on operator theoretic and geometric techniques.

The present paper aims at providing a solution to a highgain observer design problem (H-GODP) for a class of 3 × 3 inhomogeneous linear hyperbolic systems of balance laws with one measurement, written in a triangular form, as an extension of [START_REF] Kitsos | High-gain observer design for a class of integro-differential quasilinear hyperbolic systems -Application to an epidemic model[END_REF], [START_REF] Kitsos | High-gain observer design for a class of hyperbolic systems of balance laws[END_REF] for the case of systems with distinct velocities. This case presents technical difficulties, even for linear systems. Similar difficulties appear in the notable work of [START_REF] Alabau-Boussouira | Internal controllability of first order quasilinear hyperbolic systems with a reduced number of controls[END_REF], as well as in previous works of the same authors, about the internal controllability of hyperbolic systems in cascade form with reduced number of control laws. In these works, stronger regularity of the solutions is required. The approach here on observer design with reduced number of observations reveals some duality with this framework. The present paper also extends the results of Kitsos et Preprint submitted to Automatica al (2019) up to 3 × 3 strictly hyperbolic systems. Moreover, it can be noticed that, even though the system is linear, the use of high gain is instrumental, as a result of the considered domain of the hyperbolic operator. The class of systems that we study, written in an appropriate triangular form, can be found in various situations, like plug flow chemical reactors, where, by measuring the jacket temperature, we would expect to estimate the concentrations of masses of chemicals, or in linearized hyperbolic Lotka-Volterra systems with three species in cascade form (see for instance [START_REF] Zhang | Local Exponential Stabilization of Semi-Linear Hyperbolic Systems by Means of a Boundary Feedback Control[END_REF]). Stability and controllability for linear hyperbolic systems have been widely studied, see for instance Bastin and Coron (2016), [START_REF] Prieur | Boundary feedback control of linear hyperbolic systems: application to the Saint-Venant-Exner equations[END_REF], Espitia et al (2017) by employing Lyapunov and operator-theoretic strategies. In Coron and Nguyen (2019) and [START_REF] Auriol | Minimum time control of heterodirectional linear coupled hyperbolic PDEs[END_REF], the problem of minimum-time control of such systems arises, suggesting that a faster observer than a boundary one, which overcomes the transport phenomena restrictions, would be desirable.

The main contribution here is the proof of solvability of the H-GODP for 3 × 3 strictly hyperbolic triangular systems via an indirect approach. An important idea of this paper is to perform an appropriate infinite -dimensional state transformation to overcome the lack of a commutative property, yet needed in the Lyapunov stability analysis. Note that constraints on the source term can be found in some studies of stability problems as in Bastin and Coron (2016), Coron and Bastin (2015) and in [START_REF] Prieur | Stability of switched linear hyperbolic systems by Lyapunov techniques[END_REF] (Proposition 2.1), which allow a similar commutation, while this is not the case here. The transformation we use is related to the system's triangularity and maps the original hyperbolic system into a target system of partial differential equations (PDEs). This methodology requires stronger regularity for system's dynamics and additionally to the output correction terms, ouput's spatial derivatives up to order q -1 are injected in the high-gain observer dynamics, where q is an indicator of the number of different characteristic velocities. The convergence of the observer is then proven for the sup spatial norm.

The paper is organized as follows. The description of the H-GODP, the sufficient conditions and a solution to it are presented in Section 2, where Theorem 1 constitutes the main result. In Section 3 we prove Theorem 1 and in Section 4, we illustrate our methodology via a plug flow chemical reactor. Conclusions and perspectives are discussed in Section 5.

Notation: For a given w in R n , |w| denotes its usual Euclidean norm. For a given matrix A in R n×n , A denotes its transpose, |A| := sup {|Aw| , |w| = 1} is its induced norm and Sym(A) = A+A 2 stands for its symmetric part. By eig(A) we denote the minimum eigenvalue of a symmetric matrix A. By I n we denote the identity matrix of dimension n. By sgn(x) we denote the signum function sgn(x) = d dx |x|, when x = 0, with sgn(0) = 0. For given ξ : [0, +∞) × [0, L] → R n and time t ≥ 0 we use the notation ξ(t)(x) := ξ(t, x) to refer to the profile at time t. For a real-valued function f , by f we denote its first derivative. For a continuous (C

0 ) map [0, L] x → ξ(x) ∈ R n (or [0, L] x → A(x) ∈ R n×n ) we adopt the notation ξ ∞ := max{|ξ(x)| , x ∈ [0, L]} ( A ∞ := max{|A(x)| , x ∈ [0, L]} respectively). For a q -times continuously differentiable (C q ) map [0, L] x → ξ(x) ∈ R n we adopt the notation ξ q := q i=0 ∂ i x ξ ∞ .
2 Problem statement and main result

We are concerned with one-dimensional, first-order linear hyperbolic systems of balance laws, described by the following equations on a strip

Σ := [0, +∞) × [0, L] ξ t (t, x) + Λ(x)ξ x (t, x) = M (x)ξ(t, x), (1a)
where

ξ = ξ 1 ξ 2 ξ 3 : Σ → R 3 .
Consider also a distributed measurement, provided in the output, of the form

y(t, x) = Cξ(t, x), (1b)
where y is a mapping from Σ to R and the matrix of characteristic velocities

Λ(x) = diag {λ 1 (x), λ 2 (x), λ 3 (x)} ,
with λ 1 (x) > 0 for all x ∈ [0, L] (without loss of generality), contains m ∈ {1, 2, 3} positive characteristic velocities λ i (x) and 3 -m negative ones, implying that these are of constant sign in [0, L]. We have the following algebraic structures for the involved matrices

M (x) =     m 1,1 (x) m 1,2 (x) 0 m 2,1 (x) m 2,2 (x) m 2,3 (x) m 3,1 (x) m 3,2 (x) m 3,3 (x)     , C = 1 0 0 .
To give appropriate boundary conditions, let us first define a permutation matrix denoted by I π ∈ R 3×3 , which re-orders the two last components of ξ according to the signs of the related second and third characteristic velocities, namely, putting the m elements corresponding to positive velocities in ξ + and the 3 -m to negative velocities in ξ -. More explicitly,

I π ξ := ξ + ξ - .
Boundary conditions are written in the form

ξ + (t, 0) ξ -(t, L) = K ξ + (t, L) ξ -(t, 0) , ( 2a
)
where K is a matrix of the form K = K 00 K 01

K 10 K 11
, with 3-m) .

K 00 ∈ R m×m , K 01 ∈ R m×(3-m) , K 10 ∈ R (3-m)×m , K 11 ∈ R (3-m)×(
We also consider initial condition

ξ(0, x) =: ξ 0 (x), x ∈ [0, L]. (2b)
Notice that the considered system provides an internal measurement of the first state only, while not any information on the two other states is given, except for the general relationship between incoming and outgoing information on the boundaries, given by a general law (2).

Remark 1 We assumed, without loss of generality, that the first characteristic velocity λ 1 is positive and we perform all the calculations that follow for this particular case. Indeed, if λ 1 is negative, then applying the space variable transformation x → L -x, we get a hyperbolic system with opposite sign of all three velocities. Then, the observer design methodology remains unchanged for the new system. Therefore, in this paper, we may assume that the first characteristic velocity is positive.

At this point, let us define a, roughly speaking, "index of strict hyperbolicity" as follows q := min {i :

λ i ≡ λ j , ∀j = i, . . . , 3} ,
where we used the equivalence relation

λ i ≡ λ j ⇔ λ i (x) = λ j (x), for all x in [0, L]
. By this definition, we have q ∈ {1, 2, 3}, and in case of a strictly hyperbolic system, we have q = 3. The case where q = 1 (a single velocity), as a particular case of the general formulation here, has been already addressed in [START_REF] Kitsos | High-gain observer design for a class of hyperbolic systems of balance laws[END_REF], while the case where q = 1 is more complex and needs a different treatment. We further define q := max (1, q -1) .

We are now in a position to present the main assumptions.

Assumption 1 Functions λ i and m i,j belong to C q ([0, L]; R). Initial condition ξ 0 in C q [0, L]; R 3 satisfies compatibility conditions of order q (see (Bastin and Coron , 2016, Chapter 4.5.2) for definition of compatibility conditions of any order).

The following assumption concerns the structure of the matrix K linking the incoming with the outgoing information on the boundaries.

Assumption 2 If λ 3 (•) > 0, then K 11 is invertible and if λ 3 (•) < 0, then K 00 is invertible
The previous requirement is trivially met if all velocities λ 1 , λ 2 , and λ 3 are positive. Also, a case where this assumption is trivially satisfied is when K = I 3 . In that case, boundary conditions are x-periodic. Note that an x-periodicity assumption is made in Alabau-Boussouira et al ( 2017) for the somehow dual problem of controllability of underactuated systems.

The nature of the next assumption is revealed in [START_REF] Kitsos | High-gain observer design for a class of hyperbolic systems of balance laws[END_REF], where a triangular form is introduced for the case of quasi-linear hyperbolic systems. This assumption allows us to obtain a target system, that we introduce later, and is also a sufficient condition for the observer design.

Assumption 3 The following condition is satisfied for all x in [0, L]:

m 1,2 (x), m 2,3 (x) = 0.
Given the previous assumption, system (1a) satisfies some triangular structure, which presents an analogy to the finite-dimensional case (see [START_REF] Khalil | High-gain observers in nonlinear feedback control[END_REF]).

The following fact results from classical arguments borrowed from the theory of linear hyperbolic systems, combined with the manipulation of the extra regularity (see (Bastin and Coron , 2016, App. A, Chap. 4.5) and references therein).

Fact 1 Under Assumption 1, there exists a unique global solution to system (1a

)-(2) in C q [0, +∞) × [0, L]; R 3 .
The problem that we address in this work is stated in the following definition.

Definition 1 (H-GODP) The High-Gain Observer Design Problem is solvable for a system given by (1a)-( 2) with output (1b), while output's spatial derivatives of order q -1 might also be available, if there exists a wellposed observer system, which estimates the state of (1a) with a convergence speed that can be arbitrarily tuned via a single parameter (high-gain constant).

The observer design problem under consideration relies on a distributed measurement of a part of the state.

The main feature of it is the arbitrarily fast convergence rate, similarly as in the finite dimensions and to achieve this property, distributed measurement on the whole domain is assumed. We note here that, although boundary observers with the full-state measurement are preferred for practical reasons, in the present formulation, distributed measurement of part of the state might be available in many cases of distributed parameter systems. For instance, some setups include thermal cameras for chemical reactors or alternative methods, see for instance [START_REF] Zogg | Isothermal reaction calorimetry as a tool for kinetic analysis[END_REF], [START_REF] Pradere | Processing of temperature fields in chemical microreactors with infrared temperature fields[END_REF], providing the desired measurements. Also, approximations with distributed measurements within the domain would provide an approximated measurement on the whole domain. Furthermore, as indicated in the H-GODP definition, stronger regularity of the solutions to the initial systems is required for some classes of systems, coming from Assumption 1, since the observer dynamics may include higher-order spatial derivatives of the output. These spatial derivatives might be available, since their knowledge is causal, contrary to the time-derivatives which are strictly avoided. Approximations of such derivatives via kernel convolutions might be investigated more. In the following remark, we discuss the problem of solvability of the H-GODP, if instead of a distributed measurement, we had a boundary one.

Remark 2 H-GODP is not solvable in case of boundary measurement, instead of internal measurement as in (1b). First, arbitrary convergence condition would not be fulfilled, since a boundary observer for hyperbolic systems would experience a limitation with respect to convergence speed. The rate of convergence is limited by a minimal observation time which depends on the size of the domain L and the characteristic velocities λ i in that case (see [START_REF] Li | Exact boundary observability for quasilinear hyperbolic systems[END_REF] 2013) for the control design for underactuated system). In these cases, the considered boundary conditions allow a backstepping transformation of the observer error equations into a system with dissipative boundary conditions which is stable, but this is not the case here, since the present boundary conditions cannot allow us to obtain this dissipativity (see Coron and Bastin (2015) about linking dissipativity of boundary conditions with stability).

To proceed to the observer design, we introduce an indirect approach to deal with the generality of the considered hyperbolic operator. Although system (1a)-( 1b)-( 2) is written in an appropriate triangular form, as it was introduced in Kitsos et al ( 2018), [START_REF] Kitsos | High-gain observer design for a class of integro-differential quasilinear hyperbolic systems -Application to an epidemic model[END_REF], it seems that in presence of distinct characteristic velocities, we need to employ a different strategy, in order to find solutions for the H-GODP. The problem comes from the fact that the balance laws in (1a) do not allow the choice of a diagonal Lyapunov functional to be used in the stability analysis of the observer error equations.

A non-diagonal Lyapunov functional does not permit an integration by parts when taking its time-derivative, since the Lyapunov matrix and the matrix of velocities do not commute. To address this problem, we perform a transformation including spatial derivations of the state up to order q -2, in order to write the system in an appropriate form for which a Lyapunov approach is feasible. Then, for the obtained target system, we design the high-gain observer and, finally, returning to the initial coordinates, solvability of H-GODP is guaranteed.

Note that H-GODP for 2 × 2 quasilinear hyperbolic systems with two velocities has been already solved in Kitsos et al ( 2019). The latter is a simpler problem than the present one of the 3 × 3 case, since no state transformation is required. The increased difficulties with respect to the presence of distinct velocities appear in the somehow dual problems of internal controllability with reduced numbers of controls (see comments on algebraic solvability in Alabau-Boussouira et al ( 2017)) and the main technical difficulty of the present work is to tackle this problem.

Let us define a Banach space X by

X := C q ([0, L]; R) × C 1 ([0, L]; R 2 ),
equipped with the norm ξ X := ξ 1 q + ξ 2 1 + ξ 3 1 .

Assume that there exists a linear bounded and injective transformation T : (X , • X ) → (X , • X ), with bounded inverse, which maps system (1a)-( 2) into a target system ζ, as follows

ζ = T ξ; (3) with ζ 1 = ξ 1 .
The desired target system (T) of PDEs that we consider in this work satisfies the following equations on Σ, distinguishing two boundary cases, depending on the sign of λ 1 (•)

(T)            ζ t (t, x) + λ 3 (x)ζ x (t, x) = M (x)ζ(t, x) + Mζ 1 (t)(x), ζ + (t, 0) ζ -(t, L) = K ζ + (t, L) ζ -(t, 0) +K 1 ζ 1 (0) + K 2 ζ 1 (L), y ζ (t, x) = y(t, x) = Cζ(t, x), (4)
with initial condition ζ(0, x)

:= ζ 0 (x) = T ξ 0 (x),
where

M : C q-1 ([0, L]; R) → C 0 ([0, L]; R 3 ), K 1 , K 2 : C q-1 ([0, L]; R) → R n
are linear differential operators acting on ζ 1 , to be determined in the sequel, depending on the choice of T , M (x) is a matrix satisfying the same structure as M (x) and y ζ is target system's output, which remains equal to original system's output y. The existence of such a transformation T is shown in the following section.

The proposed high-gain observer for (4) satisfies the following equations on Σ ζt (t, x)

+ λ 3 (x) ζx (t, x) = M (x) ζ(t, x) -ΘN (x) y(t, x) -C ζ(t, x) + My(t)(x), (5a)
ζ+ (t, 0) ζ-(t, L) = K ζ+ (t, L) ζ-(t, 0) + K 1 y(0) + K 2 y(L), (5b)
with initial condition ζ0 (x) := ζ(0, x) (for a function ζ0 in X ), where Θ := diag θ, θ 2 , θ 3 , with θ ≥ 1 the candidate high-gain constant of the observer, to be selected precisely later. In the above equations, we considered, also, a vector gain N (•) in C q [0, L]; R 3 , selected in a way, such that, for P (•) in C q [0, L]; R 3×3 symmetric and positive definite, a Lyapunov equation of the following form is satisfied

2Sym (P (x) (M 1 (x) + N (x)C)) = -Q(x) (6)
for some symmetric and positive definite Q(x) of class C q , where M 1 (x) is derived by M (x) keeping only its 1 -diagonal, namely,

M 1 (x) =     0 m 1,2 (x) 0 0 0 m 2,3 (x) 0 0 0     .
Lyapunov equation ( 6) is solvable by a symmetric and positive definite P (x) for choice of N (x), such that M 1 (x)+N (x)C is Hurwitz for all x in [0, L]. The latter is feasible due to observability of the pair (M 1 (x), C) coming from Assumption 3. In addition, we note that the solution P (•) of ( 6), used in a Lyapunov functional that we introduce in the sequel, is never diagonal, meaning that it would not commute in general with the matrix of velocities Λ(•). This is a problem for Lyapunov-based stability analysis, which motivates for the proposed infinite-dimensional transformation.

We are now in a position to present our main result on the solvability of the H-GODP.

Theorem 1 Consider system (1a)-( 2), defined on Σ with output (1b) and suppose that Assumptions 1 -3 hold. Let also P in C q [0, L]; R 3×3 be symmetric and positive definite and let N in C q [0, L]; R 3 , both satisfying (6) for some Q symmetric and positive definite in C q [0, L]; R 3×3 . Then, the H-GODP is solvable by T -1 ζ (where ζ is the unique solution to (5)), for θ ≥ 1 as a high gain and initial condition T -1 ζ0 (x), with ζ0 satisfying zero and one-order compatibility conditions. More precisely, for every κ > 0, there exists θ 0 ≥ 1, such that for every θ ≥ θ 0 , the following holds for all t ≥ 0:

T -1 ζ(t, •) -ξ(t, •) ∞ ≤ e -κt T -1 ζ0 (•) -ξ 0 (•) X , (7)
with > 0 a polynomial in θ.

This observer convergence result is based on the existence of a transformation T as introduced in (3), that we show in the next section, along with the observer convergence of Theorem 1. We note also that in the study of internal controllability for underactuated systems, the phenomenon of loss of derivatives appears, as the regularity of the dynamics is stronger than the regularity of the control laws, whenever the velocities are distinct (see Theorem 3.1 in Alabau-Boussouira et al ( 2017)). In the present framework of the solutions to the H-GODP, the regularity of system's dynamics needs to be stronger than the regularity of the space in the norm of which the asymptotic convergence of the observer is exhibited (sup spatial norm).

Proof of Theorem 1

In this section, we prove Theorem 1. We first show the existence of the infinite-dimensional state transformation T and then we prove the observer convergence via Lyapunov analysis.

Transformation and target system

We show here the existence of T and provide the corresponding dynamics of the target system (T).

Let us choose T in (3) by

T := I 3 + T ; T :=     0 0 0 τ (x)∂ x 0 0 0 0 0     , (8)
where

τ (x) := λ 2 (x) -λ 3 (x) m 1,2 (x) .
Obviously, this transformation is bounded, invertible, with bounded inverse from X to X , independently of boundary conditions. Applying this transformation to system (1a)-( 2), we obtain target system (T) with

M (x) :=M (x) + τ (x)     0 0 0 m 1,1 (x) m 1,2 (x) 0 0 0 0     , M :=     a 1 (x)∂ x a 2 (x)∂ x + a 3 (x)∂ 2 x a 4 (x)∂ x     , K 1 := 1 + sgn(λ 1 (x)) 2 I π T C - 1 -sgn(λ 1 (x)) 2 KI π T C , K 2 :=- 1 + sgn(λ 1 (x)) 2 KI π T C + 1 -sgn(λ 1 (x)) 2 I π T C ,
where

a 1 (x) :=λ 3 (x) -λ 1 (x) -m 1,2 (x)τ (x), a 2 (x) :=τ (x)m 1,1 (x) -τ (x)λ 1 (x) -τ (x)(m 1,2 (x)τ (x)) -τ (x)m 2,2 (x) + τ (x)λ 2 (x), a 3 (x) := (λ 2 (x) -λ 1 (x) -τ (x)m 1,2 (x)) τ (x), a 4 (x) := -τ (x)m 3,2 (x).
By the above forms of mappings M and K 1 , K 2 , we see that the target system of PDEs (4) contains spatial derivatives of the first state ζ 1 up to order 2 in its dynamics and up to order 1 in its boundaries. Notice also that for the case of 3 positive velocities (m = 3) as a special case, we have I π = I 3 and mappings K 1 and K 2 take the simpler forms

K 1 :=     0 τ (0)∂ x 0     , K 2 := -K 00     0 τ (L)∂ x 0     .

Observer convergence

We are in a position to prove that the proposed observer is a high-gain observer for the target system, which is mapped from the original system via T . Observer's exponential convergence is proven for appropriate spatial norm. Injectivity of T and boundedness of its inverse, then, guarantees that T -1 ζ approaches exponentially the state ξ of the original system and, thus, solves the H-GODP.

We start from a prerequisite lemma for the wellposedness of the observer, which is a direct consequense of Assumption 1 and linearity of the system. The proof relies on the method of characteristics and a combination of classical arguments from proofs found in (Bastin and Coron , 2016, Chap. 4.5, App. A), therefore, it is omitted here.

Lemma 1 Under the regularity assumptions for the dynamics and for any y in C q ([0, +∞) × [0, L]; R), the problem described by (5) on domain Σ with initial condition ζ0 (x) satisfying compatibility conditions of order 1 admits a unique solution

ζ in C 1 [0, +∞) × [0, L]; R 3 ).
Consider, now, observer (5) for target system (T). We define a scaled observer error by ε := Θ -1 ζ -ζ , for which we derive the following hyperbolic equations on Σ

ε t (t, x) + λ 3 (x)ε x (t, x) = θ (M 1 (x) + N (x)C) ε(t, x) + Θ -1 ( M (x) -M 1 (x))Θε(t, x), (9a) ε(t, 0) = Θ -1 KΘε(t, L), when λ 3 (•) > 0, (9b)
ε(t, L) = Θ -1 KΘε(t, 0), when λ 3 (•) < 0, (9c)
where

K :=I π K 00 -K 01 K -1 11 K 10 K 01 K -1 11 -K -1 11 K 10 K -1 11 I π , K :=I π K -1 00 -K -1 00 K 01 K 10 K -1 00 K 11 -K 10 K -1 00 K 01 I π .
To prove exponential stability of with respect to its origin of the solutions to observer error system, we adopt a Lyapunov-based approach inspired by methodologies presented in Bastin and Coron (2016). Similar p-functionals have appeared in [START_REF] Kitsos | High-gain observer design for a class of hyperbolic systems of balance laws[END_REF], [START_REF] Kitsos | High-gain observer design for a class of integro-differential quasilinear hyperbolic systems -Application to an epidemic model[END_REF]. The stability is proven for the 1-spatial norm. We define a Lyapunov functional

W p : C 1 ([0, L]; R 3 ) → R by W p [ε] := L 0 π(x)e pµ θ x G p [ε](x)dx 1/p ; (10a) G p [ε](x) := ε (x)P (x)ε(x) + ε t (x)P (x)ε t (x) p (10b)
and p in N, P (•) is of class C 1 , symmetric and positive definite, satisfying (6), π : [0, L] → R is a function given by

π(x) := (π -1) x L + 1; π := max x∈[0,L] λ 3 (x) min x∈[0,L] λ 3 (x) (11)
and µ θ is a constant given by

µ θ := 1 + sgn(λ 3 (x)) 2L ln | K| 2 P (•) ∞ min x∈[0,L] eig(P (x)) θ 4 + 1 -sgn(λ 3 (x)) 2L ln min x∈[0,L] eig(P (x)) | K| 2 P (•) ∞ θ 4 . (12)
Note here that, by its definition, π is bounded as follows

1 ≤ π(x) ≤ π, ∀x ∈ [0, L]. (13)
By invoking Lemma 1 and Fact 1, we are in a position to define G p , W p : [0, +∞) → R + by

G p (t) := G p [ε](t), W p (t) := W p [ε](t), t ≥ 0. (14)
By temporarily assuming that ε has some extra regularity, i.e., it is C 2 , we obtain the following hyperbolic equations for ε t

ε tt (t, x) + λ 3 (x)ε xt (t, x) = θ (M 1 (x) + N (x)C) ε t (t, x) + Θ -1 ( M (x) -M 1 (x))Θε t (t, x), (15a) ε t (t, 0) = Θ -1 KΘε t (t, L), λ 3 (•) > 0, (15b)
ε t (t, L) = Θ -1 KΘε t (t, 0), λ 3 (•) < 0. (15c)
For the time-derivative Ẇp along the C 1 solutions of ( 15) we get the following

Ẇp = 1 p W 1-p p L 0 pπ(x)e pµ θ x G p-1 (x) × 1 i=0 ∂ i+1 t ε (x)P (x)∂ i t ε(x) + ∂ i t ε (x)P (x)∂ i+1 t ε(x) dx = W 1-p p 1 p T 1,p + T 2,p , (16)
where

T 1,p := - L 0 λ 3 (x)π(x)e pµ θ x [∂ x G p (x) -pG p-1 (x) 1 i=0 ∂ i t ε (x)P x (x)∂ i t ε(x) dx, T 2,p := L 0 (2θπ(x)e pµ θ x G p-1 (x) × 1 i=0 ∂ i t ε (x)Sym ((P (x)(M 1 (x) + N (x)C) ∂ i t ε(x) +2π(x)e pµ θ x G p-1 (x) 1 i=0 ∂ i t ε (x) ×Sym P (x)Θ -1 ( M (x) -M 1 (x))Θ ∂ i t ε(x) dx.
Using an integration by parts in term T 1,p and utilizing (13) and other trivial bounds, we get

T 1,p ≤ -λ 3 (L)π(L)e pµ θ L G p (L) + λ 3 (0)π(0)G p (0) + (ω 1 + pω 2 + p|µ θ | λ 3 (•) ∞ ) W p p ,
where

ω 1 := π-1 L λ 3 (•) ∞ + λ 3 (•) ∞ , ω 2 := λ3(•) ∞ Px(•) ∞
min x∈[0,L] eig(P (x)) . Substituting (15b), for each of the cases λ 3 (•) > 0 and λ 3 (•) < 0, the above yields

T 1,p ≤ -min x∈[0,L] λ 3 (x)π(L) min x∈[0,L] eig(P (x)) p e pµ θ L +π(0) max x∈[0,L] λ 3 (x) θ 2 | K| 2p P (•) p ∞ × 1 i=0 |∂ i t ε(L)| 2 p + (ω 1 + pω 2 +p|µ θ | λ 3 (•) ∞ ) W p p , when λ 3 (x) > 0, T 1,p ≤ -min x∈[0,L] λ 3 (x)π(L) P (•) p ∞ θ 2 | K| 2p ×e pµ θ L + π(0) max x∈[0,L] λ 3 (x) min x∈[0,L] eig(P (x)) p 1 i=0 |∂ i t ε(0)| 2 p + (ω 1 +pω 2 + p|µ θ | λ 3 (•) ∞ ) W p p , when λ 3 (x) < 0,
which, by virtue of ( 11), ( 12), is bounded as follows

T 1,p ≤ (ω 1 + pω 2 + p|µ θ | λ 3 (•) ∞ ) W p p . (17)
Next, for θ ≥ 1 and invoking (6), we obtain for T 2,p

T 2,p ≤ (-θω 3 + ω 4 ) W p p , (18)
where ω 3 := (x)) . By ( 16), in conjunction with ( 17), (18) we get

min x∈[0,L] eig(P (x)) min x∈[0,L] eig(Q(x)) P (•) ∞ , ω 4 := 2 P (•) ∞ M (•)-M1(•) ∞ min x∈[0,L] eig(P
Ẇp ≤ (-θω 3 + |µ θ | λ 3 (•) ∞ + ω 1 + ω 2 + ω 4 ) W p . (19)
We obtained the above estimate of Ẇp for ε of class C 2 , By invoking density arguments, the results remain valid with ε only of class C 1 (see Coron and Bastin (2015) for further details). Now, one can select the high gain θ, such that θ ≥ θ 0 , where θ 0 ≥ 1 is such that

-θω 3 + |µ θ | λ 3 (•) ∞ + ω 1 + ω 2 + ω 4 ≤ -2κ θ , ∀θ ≥ θ 0 ,
for some κ θ > 0. One can easily check that for any κ θ > 0, there always exists a θ 0 ≥ 1, dependent on the involved constants, such that the previous inequality is satisfied. Subsequently, (19) yields

W p (t) ≤ e -2κ θ t W p (0). ( 20
)
Taking also into account (13), we get the following prop-erty holding for continuous ε, ε t

lim p→∞ W p = lim p→∞ e µ θ • π(•) 1 p G 1/p p (•) L p (0,L) (21) = 1 i=0 e µ θ • ∂ i t ε (•)P (•)∂ i t ε(•) ∞ . (22)
By ( 20), in conjunction with (21), we derive

1 i=0 e µ θ • ∂ i t ε (•)P (•)∂ i t ε(•) ∞ ≤ e -2κ θ t × 1 i=0 e µ θ • ∂ i t ε 0 (•)P (•)∂ i t ε 0 (•) ∞ , (23)
where ε 0 (x) := ε(0, x). Now, from error dynamics (15) in conjunction with continuous differentiability of the dynamics, we can perform simple differentiations, so as to calculate constants ρ θ,i , σ θ,i , depending polynomially on θ, such that

1 i=0 ρ θ,i ∂ i x ε(•) ∞ ≤ 1 i=0 ∂ i t ε(•) ∞ ≤ 1 i=0 σ θ,i ∂ i x ε(•) ∞ .
Combining the above estimates with the following inequality

e µ θ -|µ θ | 2 L min x∈[0,L]
eig(P (x))

1 i=0 ∂ i t ε(•) ∞ 2 ≤ 1 i=0 e µ θ • ∂ i t ε (•)P (•)∂ i t ε(•) ∞ ≤ e µ θ +|µ θ | 2 L P (•) ∞ 1 i=0 ∂ i t ε(•) ∞ 2 ,
we obtain an exponential stability result as follows

ζ(t, •) -ζ(t, •) 1 ≤ ¯ θ e -κ θ t ζ0 -ζ 0 1 , (24)
where ¯ θ is a polynomial in θ (as in high-gain observers in finite dimensions). To return to the original coordinates, we notice that T is bounded from

X to X , X is continuously embedded in C 1 ([0, L]; R 3 ). Also the extension of T -1 on C 0 [0, L]; R 3 ) is bounded in C 0 [0, L]; R 3 ) and C 1 ([0, L]; R 3 ) is continuously embed- ded in C 0 ([0, L]; R 3 ).
Thereby, by (24), we can calculate constant , polynomial again in θ, such that ( 7) is satisfied with κ = κ θ .

The proof of Theorem 1 is complete.

Remark 3 Notice that although the considered system is linear, the high-gain technique is of special interest to dominate "extra terms" in the derivative of the Lyapunov function, similarly to the nonlinear terms in finitedimensional high-gain observers. In the present case, indeed, there appears a term µ θ in the Lyapunov derivative (see (17)), coming from the boundary conditions, and having an effect similar to nonlinearities in finitedimensional systems. In an abstract sense, passing from the finite dimensions to the infinite dimensions, the hyperbolic (differential) operator imposes extra difficulties, even for linear systems.

Simulation

In this section, we apply the high-gain observer design to a system of an exothermic plug flow chemical reactor. Control and observer designs for chemical reactors have been widely investigated (see for instance Boskovic and Krstic (2002) and Christophides and Daoutidis (1996)).

Here we consider a linearized model, where system's states ξ 1 , ξ 2 , ξ 3 represent the deviation with respect to their steady values, i.e., ξ

1 = T c -T * c , ξ 2 = T r -T * r , ξ 3 = C A -C *
A , where T c is the coolant temperature, T r is the reactor temperature and C A is the concentration of the chemical components (see (Bastin and Coron , 2016, 5.1.1) for more details). The hyperbolic dynamics satisfy (1a) with

Λ =     V c 0 0 0 V r 0 0 0 V r +     , M (x) =     k 0 -k 0 0 -k 0 k 0 + k 1 φ 1 (x) k 1 φ 0 (x) 0 -φ 1 (x) -φ 0 (x)    
for positive V c , V r , , k 0 , k 1 , and boundary conditions

ξ 1 (t, 0) = 2ξ 1 (t, L), ξ 2 (t, 0) = ξ 2 (t, L), ξ 3 (t, 0) = ξ 3 (t, L)
(chosen to be fictitious in order to get unstable trajectories for the sake of illustration), while φ 0

(x) = (a + b)exp -E RT * r (x) , φ 1 (x) = C * A (x) -b a+b C in A E R(T * r (x)) 2 φ 0 (x), for con- stants a, b, E, R, C in A .
The steady states satisfy the following differential equations over [0, L]:

V r d dx T * r = -k 0 (T * c -T * r ) + k 1 r(T * r , C * A ), V r d dx C * A = -r(T * r , C * A ), V c d dx T * c =k 0 (T * c -T * r ), with r(T r , C A ) = (a + b)C A -bC in A exp -E
RTr the reaction rate and boundary conditions T

* r (0) = T in r , C * A (0) = C in A , T * c (0) = T in c .
For simulation, numerical values are as follows:

T in r = 340K, C in A =
Fig. 1. Time and space evolution of system's output 0.02mol • L -1 , T in c = 293K. The length of the reactor is L = 1m, the reactive fluid velocity in the reactor is V r = 0.025m • s -1 , the coolant velocity in the jacket is V c = 1.13m • s -1 , the activation energy is E = 11250cal • mol -1 , a = 0.56s -1 and b = 0.12s -1 are rate constants, and R = 1.986cal • mol -1 • K -1 is the Boltzmann constant. We also add an artificial constant = 0.005m • s -1 in the third characteristic velocity to make system strictly hyperbolic and deal with the problem in its full generality (with = 0, we have q = 3 instead of q = 2, for = 0). Assume that measured output is the cooling temperature, i.e., y(t, x) = ξ 1 (t, x). We now follow the steps of the H-GODP described in the previous sections. In particular, we apply the transformation (8), introduced in Subsection 3.1, with τ = /k 0 . We choose some initial conditions, such that Assumption 1 is satisfied with q = 2, and now system meets all sufficient conditions for solvability of the H-GODP. We apply Theorem 1, with θ = 4 and N (x) = -1 5 10 . As expected, the convergence of the inversely transformed observer state to the unknown state ξ is guaranteed by Theorem 1. In Fig. 1 we represent the output ξ 1 . In figures 2 -4 we see the observation errors for each of the original states ξ 1 , ξ 2 and ξ 3 , after choosing arbitrary observer initial conditions satisfying appropriate compatibility conditions. In Fig. 5 we represent the time evolution of the sup spatial norm of the second observer error for different values of high-gain constant θ. For θ = 1, the observer converges very slowly, while for increasing values of θ above 4, the convergence becomes faster, but with a trade off with respect to transients (and overshoot) to be handled in practice (as expected from ( 7)). In Fig. 6 we see the second observer error when we consider a smaller gain (θ = 0.1), which is divergent, as we expected.

Conclusion

A high-gain observer for a class of 3×3 linear hyperbolic systems in triangular form with possibly distinct char- acteristic velocities has been presented, considering distributed measurement of a part of the state. This result constitutes an extension of the high-gain observer design for finite-dimensional systems to a class of hyperbolic systems and, also, an extension of previous works to the case of systems with distinct velocities. To that end, the hyperbolic system is first mapped into a target system of PDEs and an observer for this system is designed, utilizing output correction terms and injection of output spatial derivatives. The extension of this methodology to wider classes of infinite-dimensional systems will be subject to our future work.

Fig. 2 .

 2 Fig. 2. First observer error for θ = 4

Fig. 4 .

 4 Fig. 4. Third observer error for θ = 4

Fig. 5 .

 5 Fig. 5. sup spatial norm of the second observer error for various high-gain constants

G. Bastin and J.-M. Coron, Stability and boundary stabilization of 1-D hyperbolic systems, Progress in Nonlinear Differential Equations and Their Applications. Springer International Publishing, 2016. G. Besançon, D. Georges, and Z. Benayache, Sur la commande en dimension finie d'une classe de systèmes non linéaires de dimension infinie, Conf. Internationale Francophone d'Automatique, Bordeaux, France, 2006. G. Besançon, B. Pham, and D. Georges, Robust state estimation for a class of convection-diffusion-reaction systems, IFAC Workshop on Control of Systems Governed by Partial Differential Equations, Paris, France, 2013. D. Boskovic and M. Krstic, Backstepping control of chemical tubular reactors, Computers and Chemical Engineering, vol. 26, pp. 1077-1085, 2002. F. Castillo, E. Witrant, C. Prieur, and L. Dugard, Boundary observers for linear and quasi-linear hyperbolic systems with application to flow control, Automatica, vol. 49, no. 11, pp. 3180-3188, 2013. P. D. Christofides and P. Daoutidis, Feedback control of hyperbolic pde systems, AIChE Journal, vol. 42, no. 11, pp. 3063-3086, 1996. J.-M. Coron and G. Bastin, Dissipative boundary conditions for one-dimensional quasi-linear hyperbolic systems: Lyapunov stability for the C1-norm, SIAM J. Control Optim., vol. 53, no. 3, pp. 1464-1483, 2015. J.-M. Coron, R. Vazquez, M. Krstić, and G. Bastin, Local exponential stabilization of a 2x2 quasilinear hyperbolic system using backstepping, SIAM J. Control Optim., vol. 51, no. 3, pp. 2005-2035, 2013. J.-M. Coron and H. M. Nguyen, Optimal time for the controllability of linear hyperbolic systems in one dimensional space, SIAM Journal on Control and Optimization, vol. 57, no. 2, pp. 1127-1156, 2019. R. F. Curtain, Finite-dimensional compensator design for parabolic distributed systems with point sensors and boundary input, IEEE Transactions on Automatic Control, vol. 27, no. 1, pp. 98-104, 1982. M. Demetriou, Second order observers for second order distributed parameter systems, Syst. Control Lett., vol. 51, pp. 225-234, 2004. A. Deutschmann, L. Jadachowski, and A. Kugi, Backstepping-based boundary observer for a class of time-varying linear hyperbolic PIDEs, Automatica, vol. 68, pp. 369-377, 2016. F. Di Meglio, R. Vazquez, and M. Krstic, Stabilization of a system of n + 1 coupled first-order hyperbolic linear PDEs with a single boundary input, IEEE Transactions on Automatic Control, vol. 58, pp. 3097-3111, 2013. N. Espitia, A. Tanwani, and S. Tarbouriech, Stabilization of boundary controlled hyperbolic PDEs via Lyapunov-based event triggered sampling and quantization, IEEE Conference on Decision and Control, pp. 3112-3117, Melbourne, Australia, 2017. J. P. Gauthier and G. Bornard, Observability for any u(t) of a class of nonlinear systems, IEEE Transactions on Automatic Control, vol. 26, no. 4, pp. 922-926,


Wenqiong Tu

Qiang Chen
email: qiang.chen@ensam.eu

Orcid Id

Electromechanical response of multilayered piezoelectric BaTiO 3 / PZT-7A composites with wavy architecture

Keywords: Multilayered wavy composites, homogenization, piezoelectric effect, finite-volume theory, finite-element theory

Electromechanical laminated composites with piezoelectric phases are increasingly being explored as multifunctional materials providing energy conversion between electric and mechanical energies. The current work explores thus-far undocumented combined microstructural effects of amplitude-to-wavelength ratio, volume fraction, poling direction of piezoelectric phases on both the homogenized properties and localized stress/electric field distributions in multilayered configurations under fully coupled electro-mechanical loading. In particular, the Multiphysics Finite-Volume Direct Averaging Micromechanics (FVDAM) and its counterpart, an in-house micromechanical multiphysics finite-element model, are utilized to investigate the homogenized and localized responses of wavy multilayered piezoelectric BaTiO 3 / PZT-7A architectures. These two methods generate highly agreeable results. Moreover, we critically examine the convergence of the finite-volume and finite element-based approaches via the Average Stress Theorem and Average Electric Displacement Theorem. The comparison shows the finite volume-based approach possesses a better numerical convergence. This study illustrates the FVDAM's ability toward the analysis and design of engineered multilayered piezoelectric materials with wavy architecture.

Introduction

As an essential component of sensors, actuators, transducers, and energy harvesting devices under intensive research and development [START_REF] Li | Stepping piezoelectric actuators with large working stroke for nano-positioning systems: A review[END_REF][START_REF] Liu | A comprehensive review on piezoelectric energy harvesting technology: Materials, mechanisms, and applications[END_REF][START_REF] Baek | A review of the piezoelectric electromechanical impedance based structural health monitoring technique for engineering structures[END_REF], the piezoelectric material possesses the capability to convert mechanical energy to electric energy through direct piezoelectric effect or vis-a`-vis through the inverse piezoelectric effect. While monolithic piezoelectric materials have several limitations, such as brittleness, limited range of coupled properties, and pronounced directionality (Kar-Gupta and [START_REF] Venkatesh | Electromechanical response of 1-3 piezoelectric composites: Effect of poling characteristics[END_REF][START_REF] Venkatesh | Electromechanical response of 1-3 piezoelectric composites: An analytical model[END_REF], the composite piezoelectric materials offer a promising solution toward the above issues [START_REF] Venkatesh | Electromechanical response of (2-2) layered piezoelectric composites[END_REF][START_REF] Kar-Gupta | Electromechanical response of 1-3 piezoelectric composites: Effect of fiber shape[END_REF][START_REF] Sun | Layout design optimization for magneto-electro-elastic laminate composites for maximized energy conversion under mechanical loading[END_REF].

To use and design the composite piezoelectric materials effectively, it's indispensable to have certain micromechanics techniques which can take into account the piezoelectric effect under thermal, electric, and mechanical loads in an accurate and efficient manner. In past decades, the classical micromechanics models, such as Mori-Tanaka, composite cylinder/sphere assemblage, and three-phase models, which are based on single fiber embedded in another media, have been extensively developed to calculate the homogenized properties of piezoelectric composites [START_REF] Chatzigeorgiou | Piezoelectric-piezomagnetic behaviour of coated long fiber composites accounting for eigenfields[END_REF][START_REF] Jiang | An exact solution for the three-phase piezoelectric cylinder model under antiplane shear and its applications to piezoelectric composites[END_REF][START_REF] Li | Micromechanics of magnetoelectroelastic composite materials: Average fields and effective behavior[END_REF][START_REF] Xiao | Generalized selfconsistent electroelastic estimation of piezoelectric nanocomposites accounting for fiber section shape under antiplane shear[END_REF][START_REF] Li | Variational bounds for the effective moduli of heterogeneous piezoelectric solids[END_REF]. However, those models are microstructural detail-free models which only estimate the homogenized properties and do not predict the accurate internal stress/strain, as well as electric displacement distributions. To fully understand the performance of piezoelectric materials and facilitate their applications in the real world, it's important to develop effective computational micromechanics models for piezoelectric composites which can both predict accurate homogenized and local responses. The finite element approach is the major approach to create micromechanics models. However, in order to obtain the homogenized and localized responses of unidirectional piezoelectric composites, a three-dimensional unit cell analysis needs to be conducted in the commercial software, such as Abaqus, and the periodic boundary conditions of the nodal displacements and electric potential need to be enforced explicitly during the model creation process (Chen et al., 2018a). To overcome the above computational inefficiency and applicational cumbrousness, a fully anisotropic twodimensional Q8-type multiphysics quadratic element has been constructed under the generalized plane strain constraint [START_REF] Wang | Computationally-efficient homogenization and localization of unidirectional piezoelectric composites with partially cracked interface[END_REF]. Hence, a full set of homogenized moduli can be generated directly without using the time-consuming 3D analysis of unit cell problems. Meanwhile, the parametric multiphysics Finite Volume Direct Averaging Micromechanics has proven to be an effective tool to obtain an accurate homogenized and local responses of unidirectional piezoelectric composites with the added advantages of stability and efficiency [START_REF] Wang | Computationally-efficient homogenization and localization of unidirectional piezoelectric composites with partially cracked interface[END_REF].

As an important type of piezoelectric composites, the flat multilayered piezoelectric composites have been studied extensively. For instance, the multilayered piezocomposites are utilized to yield high hydrostatic piezoelectric response, to satisfy the demanding requirements of clinic ultrasound imaging [START_REF] Cannata | Development of a 35-MHz piezo-composite ultrasound array for medical imaging[END_REF], and to improve the shear piezoelectric actuation of beams and annular plates via an off-axis lamination design [START_REF] Dubey | Shear actuation mechanism and shear-based actuation capability of an obliquely reinforced piezoelectric fibre composite in active control of annular plates[END_REF][START_REF] Panda | A balanced laminate of piezoelectric fiber composite for improved shear piezoelectric actuation of beams[END_REF], as well as to create effective energy harvesting devices via parallel links [START_REF] Asano | Energy harvester for safety shoes using parallel piezoelectric links[END_REF] or a stacking design [START_REF] Kim | Improved multilayered (Bi, Sc) O3-(Pb, Ti) O3 piezoelectric energy harvesters based on impedance matching technique[END_REF]. Further, to better understand the relation between the constituent material properties and the overall response of the layered system which would enable the tailoring of composite system's performance, the electromechanical response and homogenized material properties have been extensively investigated via both analytical and numerical approaches by various researchers [START_REF] Camarena | Improved analytical homogenization of the piezoelectric macro-fiber composite: Active layer embedded among passive layers[END_REF][START_REF] Challagulla | Electromechanical response of 2-2 layered piezoelectric composites: A micromechanical model based on the asymptotic homogenization method[END_REF][START_REF] Venkatesh | Electromechanical response of (2-2) layered piezoelectric composites[END_REF][START_REF] Kim | Effective properties of multi-layered multi-functional composites[END_REF][START_REF] Wang | Effective properties of multilayered functionally graded multiferroic composites[END_REF][START_REF] Zhang | Exact impact response of multi-layered cement-based piezoelectric composite considering electrode effect[END_REF].

In contrast with the well-studied regular flat multilayered piezoelectric composites, research work on the piezoelectric multilayers with wavy architectures is limited but shows great potential. Most recently, the piezoelectric multilayers with wavy architectures are employed to develop an advanced intelligent acoustic emission (AE) and mechanical impedance hybrid sensor [START_REF] Feng | A smart acoustic emission and mechanical impedance hybrid sensor with static force detecting and dynamic measuring capabilities[END_REF][START_REF] Feng | Development of structure enhanced micromachined acoustic emission sensors with wide-bandwidth and improved sensitivity[END_REF], an efficient wearable energy harvesting system with high stretchability [START_REF] Feng | Stretchable ferroelectric nanoribbons with wavy configurations on elastomeric substrates[END_REF][START_REF] Qi | Enhanced piezoelectricity and stretchability in energy harvesting devices fabricated from buckled PZT ribbons[END_REF] and a sandwich-type piezoelectric energy harvester powered by human walking [START_REF] Han | A self-powered insole for human motion recognition[END_REF][START_REF] Singh | Piezoelectric Based Energy Harvester Embedded in Shoe for Wearable Electronics[END_REF]. It is worth mentioning that the mechanical behaviors of multilayers with wavy architectures have been extensively studied, including thermal-elastic behavior (Khatam and Pindera, 2009b), plasticity [START_REF] Katz | Plastic strain localization in periodic materials with wavy brick-and-mortar architectures and its effect on the homogenized response[END_REF]Khatam and Pindera, 2009a), finite deformation (Tu and [START_REF] Pindera | Targeting the finitedeformation response of wavy biological tissues with bioinspired material architectures[END_REF]. The above work illustrates wavy multilayers' capability in achieving excellent toughness and mimics the biomaterial's performance by designing the wavy microstructures and constituent properties. To fully exploit the potential of the wavy multilayered piezoelectric composites, it is essential to understand the relation between the constituent's material properties, material microstructures, and the macro electromechanical behavior of the wavy layered system. However, to the best knowledge of the authors, the homogenized and local responses of piezoelectric multilayers with wavy architecture have never been investigated.

The main objective of the current paper is to understand the effects of amplitude-to-wavelength ratio, volume fraction, and the poling direction of constituent on both the homogenized properties and local stress/electric field distributions of the wavy multilayered piezoelectric BaTiO 3 /PZT-7A architectures. Secondly, the paper critically examines the accuracy of finite volume and finite element-based micromechanics approaches. The remainder of the paper is organized as follows: Section 2 describes the theoretical framework of the multiphysics FVDAM, including the homogenization framework, generalized local stiffness construction, generalized stiffness matrix assembly, and homogenization. Section 3 validates the multiphysics FVDAM against an analytical solution found in the literature and then focuses on the numerical results of the wavy multilayered piezoelectric BaTiO 3 /PZT-7A architectures, demonstrating the effects of amplitude-towavelength ratio, volume fraction, the poling direction of constituent on both the homogenized properties and local stress/electric field distributions. In Section 4, the accuracy of FVDAM and finite-element based multiphysics micromechanics are accessed based on the average stress and electric displacement theory. Discussions and future work are described in Section 5. The summary and conclusions are presented in Section 6.

Multiphysics FVDAM

The details of the mechanical version of the finitevolume direct averaging micromechanics approach have been presented in a sequence of papers [START_REF] Cavalcante | Generalized FVDAM theory for elastic -plastic periodic materials[END_REF][START_REF] Cavalcante | Finitevolume micromechanics of periodic materials: Past, present and future[END_REF][START_REF] Chen | Homogenization and localization of elastic-plastic nanoporous materials with Gurtin-Murdoch interfaces: An assessment of computational approaches[END_REF]Chen et al., 2018b). The FVDAM with coupled thermo-electro-elastic capability was developed recently by Chen et al. (2018a) for periodic materials undergoing infinitesimal deformations in order to generate higher fidelity of stress and electric fields relative to the classical micromechanics theories and finiteelement techniques. For a piezoelectric material, the generalized Hooke's Law reads:

s ij = C ijkl e kl À e kij E k D i = e ikl e kl + k ik E k ð1Þ
where s ij and D i are the stress and electric displacement components, respectively, e ij and E k are the strain and electric field components, respectively. C ijkl , e ikl , and k ik denote the stiffness, piezoelectric coefficients, and dielectric coefficients, respectively. Equation ( 1) can be expressed in matrix and vector notations for an easier implementation into the FVDAM analysis:

s D ! = Z e ÀE ! ð2Þ
where s = ½s 11 s 22 s 33 s 23 s 13 s 12 T , e = ½e 11 e 22 e 33 2e 23 2e 13 2e 12

T , D = ½D 1 D 2 D 3 T , E = ½E 1 E 2 E 3 T ,

Homogenization framework

The multiphysics FVDAM framework is a simplified version of the mathematical homogenization theory. The latter is based on a systematic asymptotic analysis of periodic media whose response is characterized by governing differential equations with periodically varying coefficients that reflect the spatial variation of the microstructures [START_REF] Charalambakis | Homogenization techniques and micromechanics. A survey and perspectives[END_REF][START_REF] Yang | Determination of metamaterial parameters by means of a homogenization approach based on asymptotic analysis[END_REF][START_REF] Zhu | Numerical investigation of necking in perforated sheets using the periodic homogenization approach[END_REF]. The mathematical homogenization theory provides a consistent framework for taking into account the effect of macroscopic strain and electric field variations in a periodic material whose microstructural scale is characterized by the small parameter §. Within this framework, the macroscopic and local variations of displacements and electric potentials are described using the global and local scales x and y, respectively, connected through the small parameter § = y=x. The operator ∂=∂x i in the governing differential equations is then replaced by

∂ ∂x ! ∂ ∂x + 1 § ∂
∂y which separates the governing differential equations into different orders of § or scales.

Consistent with the concept of scale separation, a multiscale displacement

u § ð Þ i x, y ð Þ and electric potential a § ð Þ x, y ð Þ field representation based on the small para- meter §, u § ð Þ i x, y ð Þ= u 0 ð Þ i x, y ð Þ+ §u 1 ð Þ i x, y ð Þ+ O § 2 À Á a § ð Þ x, y ð Þ= a 0 ð Þ x, y ð Þ+ §a 1 ð Þ x, y ð Þ+ O § 2 À Á ð3Þ
is used to separate the governing differential equations into equations of different scale orders.

In the finite-volume homogenization theory, the displacement and electric potential field in each phase of the microstructure is partitioned into average and fluctuating contributions dependent on the global and local coordinates, x = x 1 , x 2 , x 3 ð Þand y = y 1 , y 2 , y 3 ð Þ , respectively, without considering the higher-order terms:

u q ð Þ i x, y ð Þ= e ij x j + u 0 q ð Þ i y ð Þ a q ð Þ x, y ð Þ= À E j x j + a 0 q ð Þ y ð Þ ð4Þ
where e ij and E j are macroscopic strain and electric field. u

0 q ð Þ i i = 1, 2, 3 ð
Þ and a 0 q ð Þ are displacement and electric potential fluctuations.

Herein, we limit our analysis to multilayered wavy composites with continuous reinforcements whose response may be characterized by periodically repeating material microstructure with two-dimensional periodicity, Figure 1(a). The unit cell is discretized into quadrilateral subvolumes designated by the index q ð Þ, whose location is specified by the subvolume vertices

y i, q ð Þ 2 , y i, q ð Þ 3
referred to the fixed RUC coordinate system. The qth subvolume is generated by mapping the square-sided subvolume in the reference coordinates h À j ð Þ onto a quadrilateral subvolume resident in the actual microstructure, Figure 1(b), using the Q4-type mapping function:

y q ð Þ i (x, y) = X 4 j = 1 N j (h, j)y j, q ð Þ i , i = 2, 3 ð5Þ
where N j h, j ð Þ is the shape function given by:

N 1 h, j ð Þ= 1 4 1 À h ð Þ 1 À j ð Þ, N 2 h, j ð Þ= 1 4 1 + h ð Þ 1 À j ð Þ N 3 h, j ð Þ= 1 4 1 + h ð Þ 1 + j ð Þ, N 4 h, j ð Þ= 1 4 1 À h ð Þ 1 + j ð Þ ð6Þ

Generalized local stiffness construction

The fluctuating displacements u

0 q ð Þ i i = 1, 2, 3 ð
Þand electric potentials a 0 q ð Þ in the qth subvolume are approximated using the following polynomial expansion in the reference coordinates h À j ð Þ: q) i(00) + hW (q) i(10) + jW (q) i( 01)

u 0 i (q) = W (
+ 1 2 (3h 2 À 1)W (q) i(20) + 1 2 (3j 2 À 1)W (q) i(02) a 0 k ð Þ = W q ð Þ 4 00 ð Þ + hW q ð Þ 4 10 ð Þ + jW q ð Þ 4 01 ð Þ + 1 2 3h 2 À 1 À Á W q ð Þ 4 20 ð Þ + 1 2 3j 2 À 1 À Á W q ð Þ 4 02 ð Þ ð7Þ
where W (q) i(::) are the unknown microvariables associated with each subvolume, which ensures that e 11 = e

q ð Þ
11 and

E 1 = E q ð Þ 1 due to the continuous reinforcement in x 1 direction.
The employed two-scale expansion of displacements and electric displacements produces local strain and electric field in the form:

e (q) ij = e ij + e 0 q ð Þ ij = e ij + 1 2 ( ∂u 0 i ∂y j + ∂u 0 j ∂y i ) (q) E q ð Þ i = E i + E 0 q ð Þ i = E i À ∂a 0 q ð Þ ∂y i ð8Þ
Evaluating surface-averaged interfacial displacements and electric potentials for the qth subvolume in the reference coordinates produces the surface-averaged displacements and electric potentials (the subscript q is omitted for convenience), which provides the fundamental unknowns for the parametric multiphysics FVDAM:

û0(1, 3) i = 1 2 ð + 1 À1 u 0 i (h, 71)dh û0(2, 4) i = 1 2 ð + 1 À1 u 0 i (61, j)dj â0 1, 3 ð Þ = 1 2 ð + 1 À1 a 0 (h, 71)dh â0 2, 4 ð Þ = 1 2 ð + 1 À1 a 0 (61, j)dj ð9Þ
The first-(W

q ð Þ i 10 ð Þ , W q ð Þ i 01 ð Þ ) and second-(W q ð Þ i 20 ð Þ , W q ð Þ i 02 ð Þ )
order coefficients are explicitly expressed in terms of surface-averaged displacements and electric potentials as well as the zeroth-order coefficients (W

q ð Þ i 00 ð Þ
) upon using the definition in equation ( 9). The four remaining equations for the determination of the zeroth-order coefficients are introduced by satisfying the stress equilibrium equation and conservation equations in an average sense within each subvolume:

1 V q ð V q ∂s ij ∂x j dV = 0, 1 V q ð V q ∂D i ∂x i dV = 0 ð10Þ
Evaluating the surface-averaged interfacial tractions and normal electric displacements in terms of surfaceaveraged mechanical displacements and electric potentials produces a system of generalized local stiffness equations which can be symbolically expressed as:

Ŷ q ð Þ = K q ð Þ Û q ð Þ + N q ð Þ Z q ð Þ X ð11Þ
where

Ŷ q ð Þ = p 1 ð Þ p 2 ð Þ p 3 ð Þ p 4 ð Þ Â Ã q ð Þ T , p p, q ð Þ = t1 t2 t3 f  à p, q ð ÞT , t i = s ij n j , f = D i n i ,
and n 1 n 2 n 3 ½ are the unit normal vectors that define the orientation of each face of a quadrilateral subvolume.

Û q ð Þ = û0 1 ð Þ û0 2 ð Þ û0 3 ð Þ û0 4 ð Þ h i q ð ÞT , û0 (p, q) = û0 1 û0 2 û0 3 â0 ½ (p, q) T , X = e 11 e 22 e 33 2 e 23 ½ 2 e 13 2 e 12 À E 1 À E 2 À E 3 T , N q ð Þ = n 1 ð Þ n 2 ð Þ Â n 3 ð Þ n 4 ð Þ q ð ÞT , n p, q ð Þ = 0 0 0 n 3 n 2 0 0 n 2 0 n 3 0 0 0 0 0 n 3 n 2 0 0 0 0 0 0 0 0 0 n 2 n 3 2 6 6 4 3 7 7 5 p, q ð Þ . Z q ð Þ is
the generalized local stiffness matrix of the q th subvolume:

Z q ð Þ = C e T e Àk ! q ð Þ
. The generalized local stiffness matrices K q ð Þ given explicitly by the subvolume geometry and material properties can be found in reference (Chen et al., 2018a).

Generalized stiffness matrix assembly

The unknown interfacial surface-averaged displacements and electric potentials are determined by solving a global system of equations generated by first imposing traction and electric displacement continuity conditions at common faces of adjacent subvolumes, followed by direct enforcement of fluctuating mechanical displacement and electric potential continuity conditions. The resulting global system of equations can be symbolically represented by:

K Û0 = DZ X ð12Þ
where K denotes the generalized global stiffness matrix, the vector Û0 is comprised of all the unknown interfacial mechanical displacements and electric potentials. The matrix DZ represents the differences in the generalized local stiffness matrices of adjacent subvolumes.

Homogenization

The solution of the global system of equations, equation ( 12), yields the fluctuating displacement components and electric potentials, hence the local strains and electric fields in equation ( 8). The generalized electroelastic concentration matrices relating averaged subvolume strains and electric fields and their applied counterparts [START_REF] Hill | Elastic properties of reinforced solids: Some theoretical principles[END_REF], also noted as localization relations, are established as:

e E ! q ð Þ = A q ð Þ e E ! ð13Þ
where the matrices A q ð Þ can be determined by successive imposition of unit macroscopic strain or electric field at a time and then solving the corresponding boundary value problems to obtain average subvolume strains and electric fields at the specified loading condition.

The homogenized stresses and electric displacements for the unit cell of a volume V are expressed as a weighted sum of the averaged stresses and electric displacements over all elements:

s = 1 V ð s x ð ÞdV = 1 V X Nq q = 1 ð Vq s q ð Þ x ð ÞdV q ð Þ = X Nq q = 1 v q ð Þ s q ð Þ D = 1 V ð D x ð ÞdV = 1 V X Nq q = 1 ð Vq D q ð Þ x ð ÞdV q ð Þ = X Nq q = 1 v q ð Þ D q ð Þ ð14Þ
where V q ð Þ is the volume of the qth subvolume, and

v q ð Þ = V q ð Þ
V is the corresponding volume fraction. Taking the localization relations in the expression for the averaged composite stresses and electric displacements, in conjunction with the element-averaged constitutive relation of qth element, we obtain the global (or effective) constitutive equation for the multiphase piezoelectric composite in the form:

s D ! = C à e ÃT e à Àk à ! e À E ! ð15Þ
where C à , e à , and k à are the effective elastic, piezoelectric and dielectric matrices, respectively, given by:

C à e ÃT e à Àk à ! = X N q q = 1 v q ð Þ C e T e Àk ! q ð Þ A q ð Þ ð16Þ
3. Numerical results

Homogenized responses

In this section, the multiphysics FVDAM is employed to investigate the homogenized stiffness, piezoelectric and dielectric coefficients of multilayered piezoelectric composites. The effectiveness of the multiphysics FVDAM will be validated extensively, first against the theoretical values found in the literature, and then an in-house finite-element homogenization technique.

A flat multilayer can be considered as a special case of the wavy multilayers when the amplitude-towavelength ratio is zero. The multiphysics FVDAM is first validated against a theoretical solution based on the asymptotic homogenization model (AHM) developed by Kar-Gupta and Venkatesh (2013). Figure 2 illustrates the unit cell of a flat piezoelectric composite, consisting of two different phases, that is, the PVDF (white layer) and PZT-7A (gray layer). The elastic, piezoelectric, and dielectric coefficients which can be found in reference (Kar-Gupta and Venkatesh, 2013) are given in Table 1. The large property contrast between the constituent phases yields large stress and deformation gradients in the vicinity of the interface, providing a critical assessment of the multiphysics FVDAM. Figure 3 describes the comparison of the homogenized piezoelectric and dielectric coefficients generated by the multiphysics FVDAM and asymptotic homogenization model as a function of PZT-7A volume fraction. The results show an excellent level of accordance between the two approaches, providing a strong support for the multiphysics FVDAM approach.

We proceed to investigate the homogenized response of a multilayered composite consisting of PZT-7A and BaTiO 3 piezoelectric phases of different volume fractions and poling directions. Figure 4 illustrates the geometric representation of a repeating unit cell of a sinusoidal shape multilayered composite. The gray layer is occupied by the PZT-7A while the two white layers are both occupied by BaTiO 3 . The volume fraction of the PZT-7A phase is denoted as y f . Both of them are transversely isotropic with the axis of symmetry oriented in the x 1 direction, perpendicular to x 2 À x 3 plane. According to the relatively poling direction of each ply, the piezoelectric multilayered composites can be either ''transverse'' or ''longitudinal.'' In this work, the longitudinally poled multilayers are defined when all plies are poled along the x 1 direction, whereas the transversely poled multilayers are defined when the PZT-7A ply is poled in the x 1 direction while the BaTiO 3 plies are poled in the x 2 direction.

Figure 5 presents the variation of selected homogenized elastic moduli of the PZT-7A/BaTiO 3 multilayered composite system as a function of PZT-7A volume content over a wide range of volume fractions for three different amplitude-to-wavelength ratios, namely, 0% (flat), 10%, and 20%. It should be noted that the poling directions of PZT-7A and BaTiO 3 phases are both in x 1 direction. For validation purposes, the finite-element predictions are enclosed in the figure for comparison. We note that the imposition of the periodic boundary conditions is a difficult task for multilayered wavy architectures with curved boundaries in the commercial finite-element tools. Therefore, an in-house finite-element method (FEM) developed by the present authors (shown in the Appendix) is employed to generate reference results for comparison with the multiphysics FVDAM predictions. The inhouse multiphysics finite-element method mimics the unit cell solution methodology adopted by the multiphysics FVDAM theory. Hence the two sets of results may be compared on an equal footing. It is important to mention that both multiphysics FVDAM and inhouse FEM results are generated using MATLAB codes. The displacements, electric potentials, local stress, and electric fields are calculated directly without interpolation across element or subvolume interfaces. This is in contrast to the commercial finite-element codes that typically use averaging techniques in calculating local responses.

As expected, the homogenized moduli

C Ã 22 , C Ã 33 , C Ã 44 , C Ã
55 , and C Ã 66 decrease with increasing PZT-7A volume fractions while C Ã 23 exhibits an opposite trend. This is due to the fact that the material stiffness coefficients of PZT-7A phase are smaller than those of BaTiO 3 phase, except for C Ã 23 . The ply waviness has a pronounced effect on the homogenized moduli C Ã 33 , C Ã 55 , and C Ã 66 . C Ã 33 and C Ã 66 decrease with the increase of amplitude-towavelength ratio for a given PZT-7A volume fraction. In contrast, C Ã 55 increases with the increase of amplitude-to-wavelength ratio for a given volume fraction. Increasing the PZT-7A volume fraction first accentuates the differences due to the ply waviness. After reaching a maximum just below the volume fraction of 0.5, the ply waviness effect becomes less remarkable with increasing PZT-7A volume fraction. The multiphysics FVDAM predictions coincide with the finite-element results in the entire volume fraction range and for all the amplitude-to-wavelength ratios, providing good support for the developed FVDAM technique.

Comparisons of homogenized piezoelectric and dielectric coefficients generated by the multiphysics FVDAM and finite-element approaches as a function of PZT-volume fraction for three different ply amplitude-to-wavelength ratios are provided in Figures 6 and7, respectively. Similar to the previous observations, no visible differences are observed between the finite-element and finite-volume predictions over a wide range of volume fractions. While the ply waviness has a negligible effect on the axial piezoelectric and dielectric coefficients, namely, e In order to demonstrate the effect of the poling direction on the homogenized piezoelectric and dielectric coefficients of the piezoelectric multilayered composites, a transverse PZT-7A/BaTiO 3 multilayered composite is considered. Table 2 presents a detailed comparison of the complete set of piezoelectric and dielectric coefficients of 40% PZT-7A/BaTiO 3 multilayered composites with different poling directions. The ply amplitude-to-wavelength ratios are prescribed as 10%. It is noted that different results may be obtained by changing the relative poling direction of the PZT-7A and BaTiO 3 plies. The transverse isotropy ceases to hold in the case of transversely poled composites. The most interesting behavior is the moduli e à 16 , e à 21 , e à 22 , e à 23 , e à 34 and k à 12 , k à 21 that are present in the transversely poled composites but not in either of the individual phase or the longitudinal composites. Therefore, those coefficients are new material properties that are exhibited by composites with designated poling directions.

Besides the above-mentioned new material properties that do not exist in the longitudinal material system, certain homogenized piezoelectric properties of the transversely poled multilayered PZT-7A/BaTiO 3 composites vary much more significantly as the change of volume fraction of PZT-7A and amplitude-towavelength ratio. For instance, in the flat longitudinally poled composites as shown in Figure 6,e à 26 varies within the range of [9.7-11.4] C/m 2 and the difference between lower bound and higher bound is 1.7 C/m 2 , while in the flat transversely poled multilayers as shown in Figure 8,e à 26 varies within the range of [0-8.125] C/ m 2 with a much larger difference between the lower bound and higher bound, 8.125 C/m 2 . Meanwhile, in the longitudinally poled multilayers, the most significant change of e à 26 due to amplitude-to-wavelength ratio variation is 0.2 C/m 2 , occurring at the volume fraction of 0.35. In contrast, in the transversely poled composites, the most significant change of e à 26 due to the amplitude-to-wavelength ratio is much larger, that is 1.25 C/m 2 , occurring at the volume fraction of 0.5. e à 35 also has a similar trend. In the flat longitudinally poled multilayers, e à 35 varies in the range of [9.8-11.5] C/m 2 with 1.7 C/m 2 difference between the lower and upper bound. In contrast, for the flat transversely poled multilayers, e à 35 varies in the range of [0-9.3] C/ m 2 with 7.5 times higher difference. Meanwhile, in the longitudinally poled multilayers, the highest difference of e à 35 due to amplitude-to-wavelength ratio variation is 0.18 C/m 2 , occurring at volume fraction of 0.25. In contrast, in the transversely poled multilayers, the highest difference of e à 35 due to amplitude-towavelength ratio variation is nine times higher, occurring at volume fraction of 0.45. As we can see, in transversely poled multilayers, the volume fraction and amplitude-to-wavelength ratio have a significant effect on the homogenized piezoelectric properties, e à 26 , e à 35 and offer a large design space to tailor the homogenized piezoelectric properties. Figures 7 and9 show the variations of selected homogenized dielectric properties of longitudinally and transversely poled multilayers at different volume fractions for three different amplitude-to-wavelength ratios: 0%, 10%, and 20%. Dielectric properties, k à 22 and k à 33 of longitudinally poled multilayers are more sensitive to the change of amplitude-to-wavelength ratio than the transversely poled multilayers. However, the biggest difference between the longitudinally and transversely poled multilayers is the introduction of a coupling dielectric property, k à 12 , which reaches a maximal value at the volume fraction of 0.38. In particular, the nonlinearity of k à 12 increases significantly as the increase of amplitude-to-wavelength ratios. The new findings open new opportunities to design new applications of piezoelectric devices.

Table 3 shows comparison of the selected homogenized properties of transversely-poled PZT-7A/BaTiO 3 multilayered composites generated by the multiphysics FVDAM and FEM approaches at two different amplitude-to-wavelength ratios, that is, 10% and 20%. The relative errors defined by: ( :

ð Þ FVDAMÀ : ð Þ FEM )= :
ð Þ FEM j 3 100% are also included in Table 3. The PZT-7A phase volume fraction is fixed at 40%. As observed, the multiphysics FVDAM and FEM agree to two or three significant digits for all cases. The relative percentage errors of the multiphysics FVDAM predictions relative to the FEM results are negligible. The above comparison provides additional evidence of the accuracy of the multiphysics FVDAM for simulating the multilayered piezoelectric architectures. It should be mentioned that the finite-element homogenized properties are generated based on the same mesh discretization as the multiphysics FVDAM technique, namely, 120 3 40. We have verified that for the finiteelement homogenization, this mesh refinement produces results with sufficient accuracy, hence the finiteelement homogenized properties are reliable as reference solutions.

Local stress/electric fields

The ability to correctly recover the full-field localized stress and electric displacement distributions within the composite microstructures helps to understand the underlying deformation mechanisms of the composite materials and structures. To the best knowledge of the present authors, there have been very few reliable tools that yield accurate localized field distribution within the heterogeneous microstructures with only two exceptions, namely, the finite-element and finite-volume techniques that explicitly impose periodic boundary conditions. cf. Chen et al. (2018a) and Chen andWang (2018, 2020). Figure 10 displacements which are nonuniform within the wavy microstructures is significant. The correlation between the two computational approaches is remarkable, providing good support for the developed multiphysics FVDAM technique.

Figure 11 illustrates the comparison of the local axial shear stress s 12 and s 13 distributions during the unidirectional electric displacement loading by D 2 = 1C=m 2 , reflecting the coupling effect between the mechanical response and applied electric field. Once again, the differences between the multiphysics FVDAM and finite element predictions are small.

The effect of the ply amplitude-to-wavelength ratio on the localized stress and electric displacement distributions is illustrated in Figure 12. We consider a unidirectional multilayered composite comprising 40% PZT-7A for three different amplitude-to-wavelength ratios. The composites are subjected to unidirectional electric displacement loading D 2 = 1C=m 2 . In the case of a flat multilayer, the shear stress s 12 and electric displacement D 2 are uniform within each phase, whereas the shear stress s 13 and electric displacement D 3 are zero, indicating no coupling in s 13 and D 3 . In contrast, in the case of a wavy multilayer, the shear stress s 13 and electric displacement D 3 arise and the local stress and electric displacement become nonuniform under the same loading condition. As the ply amplitude-towavelength ratio increases, larger stress and electric field concentration are observed between the interface of different plies. Hence the location and conditions for the onset of mechanical failure in the piezoelectric composites could be predicted.

Assessment of accuracy

Multiphysics FVDAM and FEM are fundamentally different in the way of satisfying equilibrium and conservation equations in the discretized domain of the unit cell. In FVDAM, equilibrium and conservation equations are directly satisfied in an integral sense for every subvolume, hence the automatic satisfaction of global equilibrium and conservation. In contrast, in FEM, equilibrium and conservation equations are achieved via the minimization of total potential energy, which requires sufficient mesh refinement to achieve Table 3. Quantitative comparison of the homogenized properties of a transversely-poled PZT-7A/BaTiO 3 multilayered composite with 40% PZT-7A volume content for two different amplitude-to-wavelength ratios. global equilibrium and conservation, and then equilibrium and conservation at the element level. To gauge the satisfaction of equilibrium and conservation of an element or a subvolume for current piezoelectric analysis, unbalanced average stress and unbalanced average electric displacement can be defined via the Average Stress Theorem and the Average Electric Displacement Theorem. According to the Average Stress Theorem, for an arbitrary domain in equilibrium status, the average stress ŝij can be calculated by integrating the surface tractions as follows:

C à 22 (GPa) C à 44 (GPa) e à 12 (C/m 2 ) e à 22 (C/m 2 ) e à 35 (C/m 2 ) k à 11 (nC/Vm) k à 12 (nC/Vm) k à 33 (nC/Vm)
ŝij = 1 V ðð S T i x j dS = 1 V ðð S s ik n k x j dS ð17Þ
Converting the surface integral to volume integral using Gauss's Theorem, we obtain

ŝij = 1 V ðð S s ik n k x j dS = 1 V ððð V ∂ ∂x k s ik x j À Á dV = 1 V ððð V ∂s ik ∂x k x j + s ij ! dV ð18Þ
In the absence of body forces, when equilibrium is satisfied in a pointwise sense, ∂s ik ∂x k = 0, the above integral reduces to the standard definition for the volume-averaged stress, as follows:

1 V ððð V ∂s ik ∂x k x j + ∂x j ∂x k s ik ! dV = 1 V ððð V s ij dV = s ij ð19Þ
It should be noted that in FEM the equilibrium at the subdomain is not automatically satisfied unless the mesh discretization is sufficiently refined. Therefore, we can define unbalanced average stress to gauge the numerical convergence as follows:

D s ij = s ij À ŝij = 1 V ððð V s ij dV À 1 V ðð S T i x j dS ð20Þ
For each element or subvolume, the unbalanced average stress can be expressed as: D s

q ð Þ ij = 1 V q ððð V q s ij dV À 1 V ðð S q T i x j dS ð21Þ
While for the entire unit cell, a global unbalanced average stress can be calculated by summing up all the absolute value of unbalanced average stress for each element or subvolume as follows:

D s ij = X N q q = 1 V q V D s q ð Þ ij ð22Þ
Similarly, the unbalanced average electric displacement can be defined and utilized to gauge the convergence of the electric field. According to the Average Electric Displacement Theorem, for an arbitrary domain in the absence of electric charge, the average electric displacement Di can be calculated by integrating electric displacements along the domain surface as follows:

Di = 1 V ðð S D i n i x j dS ð23Þ
Converting the surface integral to volume integral using Gauss's Theorem, we obtain

Di = 1 V ðð S D i x j À Á n i dS = 1 V ððð V ∂D i ∂x i x j + D i ! dV ð24Þ
In the absence of electric charge, when the conservation equation is satisfied in a pointwise sense, ÐÐÐ V ∂D i ∂x i = 0, the above integral reduces to the standard definition for the volume-averaged electric displacement, as follows:

1 V ððð V ∂D i ∂x i x j + D i ! dV = 1 V ððð V D i dV = D i ð25Þ
In FEM the conservation at the subdomain is not automatically satisfied unless the mesh discretization is sufficiently refined. Therefore, we can define an unbalanced average electric displacement to gauge the numerical convergence as follows:

D D i = D i À Di = 1 V ððð V D i dV À 1 V ðð S D i n i x j dS ð26Þ
For each subvolume/element, the unbalanced electric displacement can be expressed as:

D D i = D i À D q ð Þ i = 1 V q ððð V q D i dV À 1 V q ðð S q D i n i x j dS ð27Þ
While for the entire unit cell, a global unbalanced average electric displacement can be calculated by summing up all the absolute value of unbalanced average electric displacement for each element or subvolume as follows:

D D i = X N q q = 1 V q V DD q ð Þ i ð28Þ
To illustrate the degree of equilibrium and local conservation is satisfied as a function of mesh refinement in the FEM and FVDAM approaches, the piezoelectric multilayers are analyzed under electric displacement loading D 2 of 1 C/m 2 . Figure 13 D 2 changes from the value of 18.6% and 12% to the value of 1.71% and 1.14%, respectively. It is important for researchers and engineers to notice that a mesh refinement study needs to be conducted toward the analysis of Multiphysics Piezoelectric Composites if the FEMbased Multiphysics micromechanics homogenization technique is utilized. In contrast, in multiphysics FVDAM, both the unbalanced average stress and unbalanced average electric displacement remain zero and they do not change as the mesh refinement. This is because in the FVDAM framework the equilibrium and local conservation equations are enforced during the process of analytical evaluation of the subvolume's generalized local stiffness matrices (Chen et al. 2018a).

Figures 14 and 15 present the distributions of normalized unbalanced average stress Dŝ 13 and electric displacement D D2 at element/subvolume level for three types of mesh refinement, 30 3 8, 60 3 14, and 120 3 40, in FEM and FVDAM. And Dŝ 13 is normalized by the maximal stress, s max 13 in the analysis domain and D D2 is normalized by the applied electric displacement, D 2 . In FVDAM, the unbalanced average stress and electric displacement for the entire domain are exactly zeros due to the exact satisfaction of the equilibrium equation and conservation equation at the subvolume level. In contrast, in FEM, the value of normalized unbalanced average stress and electric displacement varies significantly, in the range of [0%-30%]. The unbalanced average stress and electric displacement are more widespread in the unit cell with coarser mesh than the unit cell with refined mesh. It is interesting to observe the unbalanced average stress and electric displacement concentrates along the interface of two layers. This is consistent with the fact that local stress and electric displacement concentration occurs along the interface of different constituents with distinctive material properties, and equilibrium and conservation are more difficult to be satisfied in the FEM framework due to the large gradient of stress and electric displacement fields. Therefore, the multiphysics FVDAM can serve as a more reliable tool to predict the onset of failure due to large mechanical stress in the piezoelectric composites without the concern of mesh refinements.

Discussions and future work

The increasing applications of the piezoelectric multilayered composites in emerging technologies require an accurate yet efficient micromechanics model to study their homogenized and localized responses such that laborious experimental test can be avoided or reduced. Relative to the numerical technique such as the finiteelement method, the strengths of the present multiphysics FVDAM are several folds. First, the FVDAM satisfies the local (hence global) governing differential equations (stress equilibrium equations and conservation equation) in each volume at each level of mesh refinement. The finite-element method, however, is based on the minimization of the global potential energy of the analysis domain. Therefore, sufficient mesh refinement may be required to ensure local equilibrium and conservation conditions. Secondly, continuity conditions are applied on both interfacial tractions and displacements, as well as electric displacements and potentials, across all common subvolume interfaces, whereas only continuity of displacements and electric potentials are enforced at the adjacent nodes in the finite-element techniques. In the case of multiphased and heterogeneous materials, the satisfaction of continuity of both tractions and displacements, as well as electric displacements and potentials, is both demanding and essential, especially along the interface separating different phases. Thirdly, as demonstrated recently by the present authors, cf. Tu and Chen (2020), the FVDAM is naturally well-suited for simulating the damage initiation and evolution in composite materials without experiencing the convergence and numerical issues common to the finite-element approaches. The ability to simulate evolutionary damage sets the FVDAM apart from other micromechanics models, as damage simulation used to largely rely on the finiteelement techniques. The present work provides a foundation for the future extension for simulating damage initiation and evolution in piezoelectric composites.

Further, the usability of the FVDAM is at least twofold. The FVDAM can be easily used in a standalone manner in the parametric studies aimed at understanding the microstructure-homogenized structure relationship, as conducted in the present work. Secondly, to make full use of the potential offered by the proposed method, the FVDAM can be incorporated into a structure analysis code for multiscale applications. For instance, the structure version of the finite-volume theory is incorporated into topology optimization code for compliance minimization of continuum linear elastic structures, cf. [START_REF] Araujo | Checkerboard free topology optimization for compliance minimization applying the finite-volume theory[END_REF]. The multiphysics FVDAM with some modification is also suitable for incorporation into a structural analysis tool, such as Abaqus via the user material subroutine. This work is in progress and will be reported elsewhere in our future publications.

Summary and conclusion

With the aid of multiphysics FVDAM and in-house micromechanical multiphysics finite-element model, a comprehensive study is conducted to understand the effects of poling direction, volume fraction, and amplitude-to-wavelength ratio on the homogenized and localized response of wavy multilayered piezoelectric BaTiO3/PZT-7A architectures. Both approaches generated highly agreeable results, while multiphysics FVDAM proves a better convergence in satisfying the conservation equations. The principal conclusions of this study are as follows:

1. In the longitudinally poled multilayers, the increase of PZT-7A volume fraction first accentuates the variation of homogenized mechanical properties of longitudinally poled multilayers due to the amplitude-to-wavelength ratio, while 21 . Meanwhile, for the transversely poled multilayers, the amplitude-to-wavelength ratio has a significant effect on the homogenized piezoelectric properties, e à 26 , e à 35 and the dielectric properties, k à 12 , k à 21 , and offers a higher degree of tunability for the material system. 3. The concentration of stress and electric displacement field typically occurs at the interface of PZT-7A/BaTiO 3 layers. As the increase of amplitude-to-wavelength ratio, the stress and electric displacement concentration also increase. And the multiphysics FVDAM and FEM approaches enable us to predict the onset of mechanical failure. 4. FEM-based Multiphysics mechanics is more sensitive to mesh refinement. Multiphysics FVDAM exhibits a better convergence due to its direct satisfaction of the conservation equation at the subvolume level in an integral sense.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article. x, y ð Þ= e ij x j + u 0

i q ð Þ x, y ð Þ a q ð Þ x, y ð Þ= À E j x j + a 0 q ð Þ x, y ð Þ ðA1Þ
The fluctuating parts of the displacement and electric potential in each element are approximated using the Q8-type interpolation function and nodal quantities:

u 0 i q ð Þ = X 8 p = 1 N p h, j ð Þu o p, q ð Þ i , a 0 q ð Þ = X 8 p = 1 N p h, j ð Þa o p, q ð Þ ðA2Þ
where u o p, q ð Þ i and a o p, q ð Þ are the nodal displacements and electric potential. The index p indicates the node number at the four corners and four midpoints of the actual element starting from the lower-left corner and processing counterclockwise.

The employed two-scale expansion of displacements and electric displacements produces local strain and electric field in the form:

e q ð Þ x, y ð Þ= e + B q ð Þ u u 0 q ð Þ o , E q ð Þ x, y ð Þ= E À B q ð Þ u ϕ 0 q ð Þ o ðA3Þ
where u 0 q ð Þ o and ϕ 0 q ð Þ o are vectors including the fluctuating nodal displacements and electric potentials, respectively. B q ð Þ u is the standard strain-displacement matrix of the q th element while B q ð Þ u represents the gradient vector relating the electric-field and electric potentials. For a piezoelectric solid within a domain V , the total potential energy of the system is the summation of bulk strain and electric energies and the work done by the external force and electric charge:

p = U b + U e À W = 1 2 ð V s T e À D T E À Á dV À s T e À D T E ðA4Þ
Using the generalized Hooks' law, equation (A4) can be further expressed as:

p = U b + U e À W = 1 2 ð V Ce À eE ð Þ T e À e T e + kE À Á T E h i dV À s T e À D T E ðA5Þ
Substituting equation (A3) into equation (A5) and making the first variation yields:

dp = du 0 T o K uu u 0 o + K uu ϕ 0 o À f u À Á + dϕ 0 T o K T uu u 0 o + K uu ϕ 0 o À f u ðA6Þ
where

K uu = Ð V B T u CB u dV , K uu = Ð V B T u CB u dV , f u = Ð V B T u e T EÀB T u Ce dV and f u = Ð V ÀB T u eeÀB T u k E dV .
The global system of equations is obtained by dp=0:

K uu K uu K T uu K uu ! u 0 o ϕ 0 o ! = f u f u ! ðA7Þ

Figure 1 .

 1 Figure 1. (a) Discretization of repeating unit cell of multilayered piezoelectric composites with sinusoidal wavy architecture and (b) mapping transformation in parametric FVDAM.

Figure 2 .

 2 Figure 2. Unit cell of a bilayered piezoelectric composite used for comparison with the asymptotic homogenization model (Kar-Gupta and Venkatesh, 2013).

Figure 3 .

 3 Figure 3. Comparison of selected homogenized piezoelectric and dielectric coefficients as a function of PZT-7A volume fraction generated by the multiphysics FVDAM and asymptotic homogenization model (Kar-Gupta and Venkatesh, 2013).

Figure 4 .

 4 Figure 4. Geometric representation of a repeating unit cell representative of multilayered composites with wavy architectures.

  presents a comparison of the local electric displacements D 2 and D 3 distributions with the imposition of unit homogenized electric displacement D 2 generated by the multiphysics FVDAM and an in-house finite-element code for the 40% PZT-7A/BaTiO 3 composites with an amplitudeto-wavelength ratio of 20%. Both the multiphysics FVDAM and finite-element are generated using the same mesh discretization, namely, 120 3 40. It is observed that the concentration of electric

Figure 5 .

 5 Figure 5. Comparison of homogenized mechanical properties of multilayered PZT-7A/BaTiO 3 architectures as a function of PZT-7A volume faction for three different amplitude-to-wavelength ratios: 0% (flat), 10% and 20%.

Figure 6 .

 6 Figure 6. Comparison of selected homogenized piezoelectric properties of multilayered PZT-7A/BaTiO 3 architectures as a function of PZT-7A volume faction for three different amplitude-to-wavelength ratios: 0% (flat), 10%, and 20%.

Figure 7 .Figure 8 .

 78 Figure 7. Comparison of selected homogenized dielectric properties of multilayered PZT-7A/BaTiO 3 architectures as a function of PZT-7A volume faction for three different amplitude-to-wavelength ratios: 0% (flat), 10%, and 20%.

Figure 9 .

 9 Figure 9. Comparison of homogenized dielectric properties of transversely poled multilayered PZT-7A/BaTiO 3 architectures as a function of PZT-7A volume faction.

Figure 10 .

 10 Figure 10. Comparison of local electric displacement distributions of longitudinally poled multilayered PZT-7A/BaTiO 3 architectures with 40% PZT-7A volume content during unidirectional electric displacement loading by D 2 = 1C m 2 generated by the FVDAM (left column) and FEM (right column): (a) D 2 and (b) D 3 .

Figure 11 .

 11 Figure 11. Comparison of local stress distributions of longitudinally poled multilayered PZT-7A/BaTiO 3 architectures with 40% PZT-7A volume content during unidirectional electric displacement loading by D 2 = 1C m 2 generated by the FVDAM (left column) and FEM (right column): (a) s 12 and (b) s 13 .

Figure 12 .

 12 Figure 12. Effect of amplitude-to-wavelength ratio on the local stress and electric displacement distributions of longitudinally poled multilayered PZT-7A/BaTiO 3 architectures with 40% PZT-7A volume content during unidirectional electric displacement loading by D 2 = 1C m 2 : flat (left column); 10% (middle column); 20% (right column): (a) s 12 and (b) s 13 (c) D 2 , and (d) D 3 .

  illustrates the variation of overall unbalanced shear stresses D s 13 , D s 12 and normalized unbalanced electric displacements D D 2 = D 2 , D D 3 = D 2 over a wide range of mesh refinement from 12 3 4 to 120 3 40 and it shows the unbalanced quantity decreases with the increase of element number in FEM. Particularly the magnitude of D s 13 , D s 12 decreases from the value of 40 MPa and 63 MPa to the value of 4 MPa and 6 MPa with the coarsest 12 3 4 and the finest mesh 120 3 40 discretizations, respectively. And the percentage of D D 2 = D 2 and D D 3 =

Figure 13 .

 13 Figure 13. Comparison of unbalanced stresses and non-conservation of electric displacements as a function of mesh refinement generated by the FVDAM and FEM approaches.

Figure 14 .

 14 Figure 14. Comparison of subvolume unbalanced stress Dŝ 13 s max 13 y ð Þ as a function of mesh refinement generated by the FVDAM (left column) and FEM (right column) approaches: (a) 30 3 8, (b) 60 3 14, and (c) 120 3 40.

Figure 15 .

 15 Figure 15. Comparison of subvolume non-conservation of electric displacements D D2 D 2 as a function of mesh refinement generated by the FVDAM (left column) and FEM (right column) approaches: (a) 30 3 8, (b) 60 3 14, and (c) 120 3 40.

12 e 13 e 14 e 15 e 16 e 21 e 22 e 23 e 24 e 25 e 26 e 31 e 32 e 33 e 34 e 35 e 36

  22 C 23 C 24 C 25 C 26 C 13 C 23 C 33 C 34 C 35 C 36 C 14 C 24 C 34 C 44 C 45 C 46 C 15 C 25 C 35 C 45 C 55 C 56 C 16 C 26 C 36 C 46 C 56 C 66

			2		3
			6 6 6 6 6 6 6 6 4		7 7 7 7 7 7 7 5 7	,
	e =	e 11 e 2 6 4	3 7 5,
	k =	2 6 4	k 11 k 12 k 13 k 12 k 22 k 23	3 7 5:
				k 13 k 23 k 33

and Z = C e T e Àk ! . To accommodate the general anisotropic properties of poled piezoceramics, the submatrices of the generalized stiffness matrix are defined in this work as: C = C 11 C 12 C 13 C 14 C 15 C 16 C 12 C

Table 1 .

 1 Electromechanical constants of constituent phases[START_REF] Venkatesh | Electromechanical response of (2-2) layered piezoelectric composites[END_REF].

		PZT-7A	BaTiO 3	PVDF
	C 11 (GPa)	131	145.5	4.84
	C 12 (GPa)	74.2	65.94	2.22
	C 13 (GPa)	74.2	65.94	2.22
	C 22 (GPa)	148	150.4	4.84
	C 23 (GPa)	76.2	65.63	2.72
	C 33 (GPa)	148	150.4	4.84
	C 44 (GPa)	35.9	42.37	1.06
	C 55 (GPa)	25.3	43.86	5.26
	C 66 (GPa) e 11 (Cm 22 ) e 12 (Cm 22 ) e 26 (Cm 22 ) k 11 (nCV 21 m 21 ) k 22 (nCV 21 m 21 )	25.3 10.99 22.324 9.31 2.081 3.984	43.86 17.36 24.322 11.4 15.1 12.8	5.26 20.11 0.004 20.002 0.07083 0.06641

à 11 , e à 12 , e à 13 , and k à 11 , it plays an important role in affecting the transverse coefficients, namely, e à 26 , e à 35 , k à 22 , and k à 33 .

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The first author W. Tu acknowledges the support of Jiangsu University Research Startup Fund for Senior Talent.

Appendix

The multiphysics finite-element method developed for comparison with the FVDAM predictions follows the same zeroth-order homogenization framework, cf. [START_REF] Wang | Computationally-efficient homogenization and localization of unidirectional piezoelectric composites with partially cracked interface[END_REF]. That is, the displacement and electric potential in the qth element are represented by a two-scale expansion involving macroscopic and fluctuating contributions: