Adrien Morel
email: adrien.morel@univ-smb.fr

Gaël Pillonnet

David Gibus

Adrien Badel

Discrete Optimal Couple Tracking: a control scheme for electrical-tuning of strong coupling piezoelectric energy harvesters

Keywords: Piezoelectricity, Vibration energy harvesting, Nonlinear electronics, Discrete Maximum Power Point Tracking, Predictive model, Gradient algorithm, Frequency tuning, Broadband

In order to track the maximum power point of piezoelectric energy harvesters, an effective approach consists in designing interfaces able to electrically tune the harvester dynamics. Such interface should exhibit -at least -two tunable parameters in order to independently optimize the harvester damping and resonant frequency, i.e., an electromechanical impedance matching of the real and imaginary parts of the electrical load. In order to optimize and control these two parameters, it is necessary to implement two control loops that simultaneously impact the harvester dynamics. Furthermore, these loops should combine quick convergence time with low-power consumption in order to react sufficiently fast to compensate for any shift in the vibration frequency and to consume only a small portion of the harvested power. In this paper, we propose a control-law based on the successive evaluation of couples of parameters in order to select the optimal one. First, we determine analytically the expression of a set of optimal couples valid for any two-parameters electrical interfaces. Thereafter, we prove that only a few well-chosen couples of parameters are sufficient to maximize the power on a large frequency band. The proposed methodology has been experimentally verified on a tunable interface, the short-circuit synchronous electric charge extraction. Combined with a strongly coupled piezoelectric energy harvester, we have been able to maintain the harvested power close to the maximal achievable power (i.e., >70% of the maximal achievable power) on a 20 Hz frequency band, with only five couples of parameters. Such approach allows to substitute a multi-parameters convergence algorithm by a singlecouple of parameters convergence algorithm, greatly simplifying the algorithm and reducing its convergence time.

Introduction

Energy harvesting has been widely investigated as a way to replace batteries for powering small electronic devices such as sensors nodes and wake-up radios [START_REF] Penella | Powering wireless sensor nodes: Primary batteries versus energy harvesting[END_REF]. In close environments, where solar radiations and thermal gradients are weak, vibration energy harvesters constitute a promising solution [START_REF] Siang | Review of vibration-based energy harvesting technology: Mechanism and architectural approach[END_REF]. Among the various mechanical-to-electrical transduction mechanisms, piezoelectricity offers a good trade-off between compacity, power, and simplicity [START_REF] Roundy | Improving power output for vibration-based energy scavengers[END_REF]. Most piezoelectric energy harvesters (PEH) are made of resonant structures that maximize the harvested power when the vibration frequency matches the PEH resonant frequency [START_REF] Williams | Analysis Of A Micro-electric Generator For Microsystems[END_REF]. However, in many environments, the frequency of the vibrations varies with time, leading to poor performances when it shifts away from the PEH resonant frequency [START_REF] Yildirim | A review on performance enhancement techniques for ambient vibration energy harvesters[END_REF]. In order to enlarge the frequency bandwidth of PEH, several approaches and techniques have been proposed in the literature [START_REF] Twiefel | Survey on broadband techniques for vibration energy harvesting[END_REF]. For instance, some researchers proposed to purposely bring nonlinearities into the behaviors of PEH in order to enhance their frequency bandwidth [START_REF] Cottone | Nonlinear Energy Harvesting[END_REF], but such nonlinear approach may lead to complicated behaviors that are amplitude-and orbit-dependent [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF][START_REF] Saint-Martin | Statistical performance assessment of a bistable vibration energy harvester[END_REF]. Another approach consists in adapting in real-time the dynamics of the PEH and to track the shifts of vibration frequency in order to always operate at resonance. Such approach can be implemented with many frequency tuning mechanisms, such as tunable magnets [START_REF] Hoffmann | A self-adaptive energy harvesting system[END_REF], screws or springs [START_REF] Niri | A passively tunable mechanism for a dual bimorph energy harvester with variable tip stiffness and axial load[END_REF]. One of them consists in the electrical adjustment of the PEH resonant frequency and damping using tunable interface circuits (which can be seen as an electromechanical impedance matching) [START_REF] Brenes | Maximum power point of piezoelectric energy harvesters[END_REF]. Such electrical tuning allows the elaboration of compact self-powered harvesting solutions that can adapt their own dynamics to the environment. However, the control of the electrical interface remains an important challenge that need to be tackled. Such approach requires the optimization of, at least, two electrical parameters in order to simultaneously tune the resonant frequency and damping and reach the PEH maximum power point (MPP) 1 . Furthermore, these parameters have to converge toward their optimal values in a minimal time and with a small power consumption.

Multi-parameters electrical interfaces that can tune the PEH resonant frequency and enlarge the harvesting bandwidth have been widely investigated in the last decade, as illustrated in Table 1. The first electrical tuning interfaces have been proposed by Wu et al. [START_REF] Wu | Tunable resonant frequency power harvesting devices[END_REF] and Liao et al. [START_REF] Liao | Optimal parameters and power characteristics of piezoelectric energy harvesters with an RC circuit[END_REF] and rely on the adjustment of a capacitive banks connected to the PEH. Such approach has been implemented by Bouhadjar et al. in 2014 with a capacitively-tuned full-bridge rectifier [START_REF] Ahmed-Seddik | Self-powered resonant frequency tuning for piezoelectric vibration energy harvesers[END_REF]. In 2014, Badel and Lefeuvre proposed the frequencytuning synchronized electric charge extraction (FTSECE) [START_REF] Badel | Wideband piezoelectric energy harvester tuned through its electronic interface circuit[END_REF]. FTSECE is a non-linear interface that allows the tuning of the resonant frequency and electrical damping by mean of a fine control of the phase-shift of the energy harvesting events and of the voltage inversion ratio. Such technique has been experimentally validated by Brenes et al. [START_REF] Brenes | Large-bandwidth piezoelectric energy harvesting with frequency-tuning synchronized electric charge extraction[END_REF] and allows a significant gain in term of power and bandwidth, compared to a passive approach. Many other tunable interfaces such as the short-circuit SECE (SCSECE) that we proposed [START_REF] Morel | Frequency tuning of piezoelectric energy harvesters thanks to a short-circuit synchronous electric charge extraction[END_REF] and the phase-shift synchronized switch harvesting on inductor (PS-SSHI) proposed by Hsieh et al. [START_REF] Hsieh | Improving the scavenged power of nonlinear piezoelectric energy harvesting interface at off-resonance by introducing switching delay[END_REF] have been developed in the last few years as interesting alternatives to the FTSECE.

In order to maximize the harvesting power on a large frequency band, the aforementioned interfaces rely on the adjustment of two electrical parameters (e.g., a capacitive bank, the phase-shift of the energy harvesting event, or the duration of a short-circuit sequence). In most papers, these parameters have been tuned manually, in order to experimentally validate the basic principle of the tunable interface without any application-ready implementation [START_REF] Brenes | Large-bandwidth piezoelectric energy harvesting with frequency-tuning synchronized electric charge extraction[END_REF][START_REF] Morel | Frequency tuning of piezoelectric energy harvesters thanks to a short-circuit synchronous electric charge extraction[END_REF]. However, the design of an automated control of these parameters is not easy: it has to combine low-power consumption (i.e., in the µW range) with a fast convergence (mss convergence, depending on the application) of multiple electrical parameters (usually two parameters to optimize both the damping and the resonant frequency of the harvester). Among the literature, only a few papers tackle the need of an automated control. [START_REF] Morel | 2 self-tunable phaseshifted SECE piezoelectric energy-harvesting ic with a 30nW mppt achieving 446% energy-bandwidth improvement and 94% efficiency[END_REF] proposes to realize such control with a low-power (≈1 µW) perturb & observe algorithm (for a single parameter interface). However, the convergence time of the algorithm (around 15-20s) might be too slow for many applications, and this algorithm might get trapped in local power maxima, hindering the potential of the tunable interface (Fig. 1.a). In [START_REF] Cai | A piezoelectric energy-harvesting interface circuit with fully autonomous conjugate impedance matching, 156% extended bandwidth, and 0.38 µW power consumption[END_REF], Cai and Manoli proposed to realize such control using a look-up table, with pre-defined couples of parameters (depending on the detected vibration frequency). Such control, because of its relative simplicity, allows to combine lowpower consumption (≈1 µW) with quasi-instantaneous convergence time (only depending on the electronics sensing the vibration frequency). However, the choice of the pre-defined couples of parameters is not explained in the paper, and is crucial in order to design such simple low-power fast-convergence algorithm.

In this paper, we propose a generic methodology in order to determine a minimal set of couples to maintain a harvested power larger than a targeted percentage of the achievable power on the largest possible frequency bandwidth, for any two-tuning electrical interfaces (Fig. 1.b). The first section reminds the electromechanical model of linear piezoelectric energy harvesters as well as the impact of the electrical interface on the harvester dynamics, and proves the need of multi-parameter control of tunable electrical interfaces. The second section presents the generic analytical derivation of the optimal couples for any tunable interface. The third section proposes to analyze the impact of the number of couples on the harvester power-frequency response and on the convergence time of the tunable parameters. Finally, the last section provides an experimental validation of the proposed methodology on the SCSECE, a two-tuning electrical interface. This experimental validation, while prototyped on a laptop and not being self-powered, validates the proposed methodology and paves the way toward low-power algorithms for PEH frequency tuning.

! ! and ! " ! ! , ! " ! ! ! , ! " " ! ! , ! " # ! ! , ! " $
Choice of an optimal couple in a table of a few pre-defined couples

! ! , ! " ! ! ! , ! " # ! ! , ! " $ ! ! , ! " " b.

Electromechanical model of the energy harvester

In the literature, most linear piezoelectric energy harvesters (PEH) are made of piezoelectric plates bonded on a mechanical resonator (e.g., a cantilever beam). Under a sinusoidal vibration, the resonator starts to oscillate, and the piezoelectric material under strain generates electrical charges that can be collected by an electrical interface. Such single degree of freedom (SDoF) system is illustrated by Fig. 2.a, with 𝑀 being the inertial mass of the resonator, 𝐾 being its stiffness, and 𝐷 being the damping factor [START_REF] Williams | Analysis Of A Micro-electric Generator For Microsystems[END_REF]. The electromechanical equations modeling the dynamic mass displacement 𝑥, the current flowing in the electrical interface 𝑖 # and the voltage across the piezoelectric element 𝑣 # can be written as [START_REF] Penella | Powering wireless sensor nodes: Primary batteries versus energy harvesting[END_REF].

with 𝑦̈= 𝛾 the ambient acceleration, and 𝐵 $ the forcing term [START_REF] Erturk | Issues in mathematical modeling of piezoelectric energy harvesters[END_REF]. 𝐶 # is the piezoelectric material capacitance, and 𝛼 is the force factor between the mechanical part and the electrical part of the PEH. This electromechanical model remains accurate as long as the strain in the piezoelectric material remains in the linear region, that the vibration frequency is close to the resonant frequency of the harvester, and that the dielectric losses are negligible [START_REF] Morel | A comparative study of electrical interfaces for tunable piezoelectric vibration energy harvesting[END_REF]. We will consider that the external vibration is a monochromatic acceleration of amplitude 𝛾 ! whose pulsation 𝜔 might slowly vary with time. In this case, the mass displacement 𝑥(𝑡) can be approximated as a cosine wave of amplitude 𝑋 ! , such that 𝛾 = 𝛾 ! cos(𝜔𝑡 + 𝜓) and 𝑥(𝑡) = 𝑋 ! cos (𝜔𝑡), with 𝜓 the phase difference between 𝛾(𝑡) and 𝑥(𝑡). The piezoelectric voltage 𝑣 # (𝑡) is a periodic function and can be expressed as a Fourier series

𝑣 # (𝑡) = 𝑎 % + ∑ [𝑎 & cos(𝑛𝜔𝑡) + 𝑏 & sin(𝑛𝜔𝑡)] '( &)*
, with 𝑎 & and 𝑏 & the 𝑛 +, cosine and sine parameters of 𝑣 # (𝑡) Fourier series. Due to the filtering effect of the mechanical resonator (first equation of (1)), the higher frequency components of 𝑣 # (𝑡) share a negligible impact on the harvester dynamics (first harmonic assumption [START_REF] Morel | A comparative study of electrical interfaces for tunable piezoelectric vibration energy harvesting[END_REF]). Hence, in order to compute the influence of the electrical interface on the PEH dynamics, only the fundamental of the piezoelectric voltage 𝑣 # | * has to be determined. As detailed in [START_REF] Brenes | Large-bandwidth piezoelectric energy harvesting with frequency-tuning synchronized electric charge extraction[END_REF] and [START_REF] Morel | A comparative study of electrical interfaces for tunable piezoelectric vibration energy harvesting[END_REF], 𝑣 # | * can be expressed as an in-phase and out-of-phase components with the mechanical displacement:

(𝑡) = -𝐵 $ 𝑦̈(𝑡) 𝑖 # (𝑡) = 𝛼𝑥(𝑡) -𝐶 # 𝑣 # ̇(𝑡) (1)
with

𝜀 -(𝜓 * , 𝜓 " ) = - . ! / " 01 # and 𝜀 2 (𝜓 * , 𝜓 " ) = 3 ! / " 01 #
the dimensionless electrically-induced damping and stiffness, respectively, and 𝜓 * and 𝜓 " two electrically-adjustable parameters (i.e., the frequency of the harvesting events [START_REF] Morel | A unified N-SECE strategy for highly coupled piezoelectric energy scavengers[END_REF], the angular phase of the harvesting events [START_REF] Lefeuvre | Power and frequency bandwidth improvement of piezoelectric energy harvesting devices using phase-shifted synchronous electric charge extraction interface circuit[END_REF], the voltage inversion ratio [START_REF] Richter | Tunable interface for piezoelectric energy harvesting[END_REF], etc.) that can be tuned in order to control the electrical damping and stiffness. Combining (2) with the first equation of ( 1) yields [START_REF] Roundy | Improving power output for vibration-based energy scavengers[END_REF].

with 𝜔 % = I𝑀/𝐾, 𝑄 = √𝐾𝑀/𝐷 and 𝑘 ! " = 𝛼 " /L𝐾𝐶 # M being the short-circuit resonance pulsation, the mechanical quality factor and the expedient electromechanical coupling of the PEH, respectively. (3) proves that the electrical interface has both an impact on the damping of the PEH (𝜀 -) and on its resonance frequency (𝜀 2 ). The analysis of (3) in the frequency domain leads to the expressions of the mechanical displacement 𝑥 and the harvested power 𝑃 45+ .

with 𝑃 67! = L𝐵 $ 𝛾 " 𝑄M/(8 𝜔 % 𝑀) the maximum power that can be harvested with a linear PEH, as proved in [START_REF] Morel | A comparative study of electrical interfaces for tunable piezoelectric vibration energy harvesting[END_REF], and Ω = 𝜔/𝜔 % the normalized vibration frequency. [START_REF] Yildirim | A review on performance enhancement techniques for ambient vibration energy harvesters[END_REF] demonstrates that the power-frequency response of any linear PEH can be adjusted using an electrical interface. The normalized power 𝑃 &89 ∈ [0,1] is expressed by [START_REF] Twiefel | Survey on broadband techniques for vibration energy harvesting[END_REF].

It is possible to harvest the power limit 𝑃 67! (𝑃 &89 = 1) for any vibration frequency and for any PEH characteristics (𝑘 ! " , 𝑄) as long as the electrical interface respects the impedance matching conditions given by ( 7) [START_REF] Erturk | Issues in mathematical modeling of piezoelectric energy harvesters[END_REF].

𝜀 - : $%# and 𝜀 2 : $%# are the values of 𝜀 -and 𝜀 2 that allow to harvest the power limit, 𝑃 67! . Tunable electrical interfaces have been recently developed in order to realize such electromechanical impedance matching, and enlarge the harvested power of PEH on large frequency ranges [START_REF] Badel | Wideband piezoelectric energy harvester tuned through its electronic interface circuit[END_REF], [START_REF] Morel | Frequency tuning of piezoelectric energy harvesters thanks to a short-circuit synchronous electric charge extraction[END_REF]. Such harvesting interfaces require the use of multidimensional control to adjust the values of (𝜓 * , 𝜓 " ) and (𝜀 -, 𝜀 2 ) to respect the matching conditions given by [START_REF] Cottone | Nonlinear Energy Harvesting[END_REF]. In this article, we propose a control methodology in order to replace the multidimensional control by a monodimensional control that can ease the convergence of (𝜓 * , 𝜓 " ), save power, and reduce the overall complexity of the electrical interface.

𝑣 # | * (𝑡) = 𝛼 𝐶 # T𝜀 2 (𝜓 * , 𝜓 " ) 𝑥(𝑡) + 𝜀 -(𝜓 * , 𝜓 " ) 𝜔 𝑑𝑥 𝑑𝑡 (𝑡)V (2) 1 𝜔 % " 𝑑 " 𝑥 𝑑𝑡 " (𝑡) + T 1 𝑄𝜔 % + 𝑘 ! " 𝜀 -(𝜓 * , 𝜓 " ) 𝜔 V 𝑑𝑥 𝑑𝑡 (𝑡) + [1 + 𝑘 ! " 𝜀 2 (𝜓 * , 𝜓 " )] 𝑥(𝑡) = - 𝛾(𝑡) 𝜔 % " 𝐵 $ 𝑀 (3) 𝑥 = -𝛾/𝜔 % " [1 -Ω " + 𝑘 ! " 𝜀 2 (𝜓 * , 𝜓 " )] + 𝑗[Ω/𝑄 + 𝑘 ! " 𝜀 -(𝜓 * , 𝜓 " )] (4)
𝑃 45+ = 𝑃 67! 4 𝑘 ! " 𝜀 -Ω/𝑄 [1 -Ω " + 𝑘 ! " 𝜀 2 (𝜓 * , 𝜓 " )] " + [Ω/𝑄 + 𝑘 ! " 𝜀 -(𝜓 * , 𝜓 " )] " (5)
𝑃 &89 = 𝑃 45+ 𝑃 67! = 4 𝑘 ! " 𝜀 -Ω/𝑄 [1 -Ω " + 𝑘 ! " 𝜀 2 (𝜓 * , 𝜓 " )] " + [Ω/𝑄 + 𝑘 ! " 𝜀 -(𝜓 * , 𝜓 " )] " (6)
⎩ ⎪ ⎨ ⎪ ⎧ 𝜀 - : $%# = Ω 𝑄 𝑘 ! " 𝜀 2 : $%# = 1 -Ω " 𝑘 ! " (7)

Principle of the discrete optimal couple tracking

The proposed scheme, named discrete optimal couple tracking (DOCT), consists in pre-determining the couples of electrically-controllable parameters (𝜓 * , 𝜓 " ) that maximize the robustness of the harvested power on a predefined frequency band [𝛺 !7& , 𝛺 !35 ], in order to implement a monodimensional gradient algorithm that chooses the best parameters couple among the pre-determined ones. First, we propose to determine the optimal couples (𝜀 -, 𝜀 2 ) maximizing the robustness of the harvested power on [𝛺 !7& , 𝛺 !35 ]. Then, the optimal couples (𝜓 * , 𝜓 " ) can be found using the interface-dependent relation between (𝜓 * , 𝜓 " ) and (𝜀 -, 𝜀 2 ).

From ( 6), the resonant frequency of the system (i.e., the vibration frequency that maximizes the PEH powerfrequency response) can be adjusted with 𝜀 2 , while the damping of the system can be adjusted with 𝜀 -. Thus, if we define an electrical interface that can emulate a number

𝑁 ; ∈ ℕ * of couples (𝜀 -,& & , 𝜀 2,& & ) (with 𝑛 ; ∈ [1, 𝑁 ; ]
representing the index of the considered couple of electrical influences), this interface can be associated with a number 𝑁 ; of power-frequency responses. For example, the power-frequency responses of an electrical interface capable of emulating three unbalanced 2 couples L𝜀 -,* , 𝜀 2,* M, L𝜀 -," , 𝜀 2," M and (𝜀 -,> , 𝜀 2,> ) is illustrated by Fig. 3.a. The maximum power achievable power for any frequencies with such interface will be called 𝑃 4&G , the power envelope, and is represented by a black dashed line in Fig. 3. This envelope is defined mathematically by [START_REF] Huguet | Drastic bandwidth enhancement of bistable energy harvesters: Study of subharmonic behaviors and their stability robustness[END_REF].

In order to maximize the robustness of the harvested power on a given frequency band [𝛺 !7& , 𝛺 !35 ], we have to maximize the minimum of the power envelope on [𝛺 !7& , 𝛺 !35 ]. Such minimum is defined by [START_REF] Saint-Martin | Statistical performance assessment of a bistable vibration energy harvester[END_REF].

In Fig. 3.a., the minimum of the power envelope 𝑃 4&G is relatively low (𝑃 4&G,!7& = 0.42). In order to increase 𝑃 4&G,!7& , 𝜀 2," can be increased, shifting the red power-frequency response on the right (higher resonance frequency), as shown in Fig. 3.b. In this case, we can observe that 𝑃 4&G,!7& = 0.6.

0 D,1 =0.24, 0 K,1 =0.1 0 D,2 =0.1, 0 K,2 =0.5 0 D,3 =0.4, 0 K,3 =0
0 D,1 =0.24, 0 K,1 =0.1 0 D,2 =0.1, 0 K,2 =0.2 0 D,3 =0.4, 0 K,3 =0.9 P env a.
b.

! !"#,%&" = 0.42 ! !"#,%&" = 0.6 𝑃 4&G (Ω) = maxp𝑃 &89 LΩ, (𝜀 -,& & , 𝜀 2,& & )M; 𝑛 ; ∈ [1, 𝑁 ; ]r (8)
𝑃 4&G,!7& = min{𝑃 4&G (Ω); Ω ∈ [Ω !7& , Ω !35 ]} (9)
Considering the linear behavior of the PEH, the minimum points 𝑃 4&G,!7& always correspond to the intersections of two power-frequency responses. For example, in Fig. 3.a., 𝑃 4&G,!7& corresponds to the intersection of the blue and yellow power-frequency responses. In Fig. 2 b), 𝑃 4&G,!7& becomes the intersection of the red and yellow powerfrequency responses. The optimal case that maximizes 𝑃 4&G,!7& can be attained when the intersections of all the power-frequency responses share the same (maximized) value and are all fairly distributed on the frequency interval, as illustrated by Fig. 4 with 𝑁 ; = 3 (balanced couples case). Hence, in order to maximize 𝑃 4&G,!7& , the first condition is that the 𝑁 ; resonant frequencies corresponding to the 𝑁 ; power-frequency responses must be evenly spaced in [𝛺 !7& , 𝛺 !35 ]. This condition is expressed analytically by [START_REF] Hoffmann | A self-adaptive energy harvesting system[END_REF].

Combining [START_REF] Hoffmann | A self-adaptive energy harvesting system[END_REF] with the matching condition on 𝜀 2 expressed by [START_REF] Cottone | Nonlinear Energy Harvesting[END_REF] leads to [START_REF] Niri | A passively tunable mechanism for a dual bimorph energy harvester with variable tip stiffness and axial load[END_REF].

If [START_REF] Niri | A passively tunable mechanism for a dual bimorph energy harvester with variable tip stiffness and axial load[END_REF] is respected, the normalized frequencies Ω % that correspond to the minima of the power envelope are given by [START_REF] Brenes | Maximum power point of piezoelectric energy harvesters[END_REF].

To maximize 𝑃 4&G,!7& , the second condition is to maximize the value of the power envelope when

Ω = Ω %,& & ∀𝑛 ; ∈ [1, 𝑁 ; ].
This condition is expressed analytically by [START_REF] Wu | Tunable resonant frequency power harvesting devices[END_REF].

Solving [START_REF] Wu | Tunable resonant frequency power harvesting devices[END_REF], we obtain the following condition on the 𝑁 ; electrical damping 𝜀 -,& & :

∀𝑛 ; ∈ [1, 𝑁 ; ], Ω 94H,& & = Ω !7& + (Ω !35 -Ω !7& ) 2𝑁 ; (2𝑛 ; -1) (10)
∀𝑛 ; ∈ [1, 𝑁 ; ], 𝜀 2,& & = Ω 94H,& & " -1 𝑘 ! " = wΩ !7& + Ω !35 -Ω !7& 2𝑁 ; (2𝑛 ; -1)x " -1 𝑘 ! " (11)
Ω %,& & = Ω !7& + (Ω !35 -Ω !7& ) 𝑁 ; 𝑛 ; (12)
∀𝑛 ; ∈ [1, 𝑁 ; ], 𝜕𝑃 &89 𝜕𝜀 - | J)J ',) & ,K * )K *,) & = 0 (13)
Hence, the optimality conditions on the 𝑁 ; discrete couples L𝜀 -,& & , 𝜀 2,& & M have been determined and are given by ( 11) and ( 14). The only parameters that are needed in order to find such optimal couples are the piezoelectric electromechanical coupling 𝑘 ! " and 𝑄. These two parameters can easily be found with an impedance analysis of the PEH, as shown in [START_REF] Badel | Wideband piezoelectric energy harvester tuned through its electronic interface circuit[END_REF]. Note that none of the 𝑁 ; optimal discrete couples L𝜀 -,& & , 𝜀 2,& & M are equal to (𝜀 - : $%# , 𝜀 2 : $%# ) because the aim of DOCT is not to harvest the power limit on the whole frequency range of interest, but to harvest a power close to this power limit with a few optimal couples.

Critical coupling and [𝛀 𝒎𝒊𝒏 , 𝛀 𝒎𝒂𝒙 ] for implementing the DOCT

Practically, the values of 𝜀 -and 𝜀 2 are bounded and depend on the choice of the electrical interface. This mean that it might not be possible to fulfill the optimality conditions on 𝜀 2,& & and 𝜀 -,& & given by ( 11) and ( 14) if the values of max(𝜀 -) and max(𝜀 2 ) (with a given interface) are too small, or if 𝑘 ! " is too small (meaning that the impact of the electrical interface does not allow to induce the optimal damping or stiffness). In order to ensure that ( 14) can be fulfilled for all 𝑁 ; , max(𝜀 -) should be greater or equal to a limit value given by [START_REF] Ahmed-Seddik | Self-powered resonant frequency tuning for piezoelectric vibration energy harvesers[END_REF].

If the case of a weakly coupled PEH, Ω !35 ≈ Ω !7& ≈ 1. Therefore, in this case, ( 15) becomes ( 16). This means that the PEH electromechanical coupling 𝑘 ! " should be greater than (𝑄 max(𝜀 -)) L* in order to implement the DOCT. 𝑘 !,; " = (𝑄 max(𝜀 -)) L* is the critical coupling which separates weakly coupled PEH (𝑘 ! " < 𝑘 !,;

" ) from strongly coupled PEH (𝑘 ! " > 𝑘 !,; " ), as defined in [START_REF] Liao | Maximum power, optimal load, and impedance analysis of piezoelectric vibration energy harvesters[END_REF]. For instance, in the case of the SCSECE interface, max(𝜀 -) = 4/𝜋 [START_REF] Morel | Frequency tuning of piezoelectric energy harvesters thanks to a short-circuit synchronous electric charge extraction[END_REF]. To implement the DOCT with the SCSECE interface, the electromechanical coupling of the PEH should be greater than 𝜋/(4𝑄). ( 16) is consistent with the definition of critical coupling and coupling state given in [START_REF] Liao | Maximum power, optimal load, and impedance analysis of piezoelectric vibration energy harvesters[END_REF] and [START_REF] Liao | Generalized modeling and analysis of piezoelectric vibration energy harvesters[END_REF]. Indeed, [START_REF] Liao | Maximum power, optimal load, and impedance analysis of piezoelectric vibration energy harvesters[END_REF] proves that such critical coupling depends on the choice of electrical interface and on the quality factor of the PEH.

In order to ensure that ( 11) is fulfilled, the pre-defined frequency band [Ω !7& , Ω !35 ] where the DOCT algorithm maximizes the harvested power should be chosen such that (17) is respected. [START_REF] Brenes | Large-bandwidth piezoelectric energy harvesting with frequency-tuning synchronized electric charge extraction[END_REF] shows that the stronger the electromechanical coupling 𝑘 ! " , the larger the frequency band where the DOCT algorithm can be used.

For instance, with a RC load (min(𝜀 2 ) = 0, max(𝜀 2 ) = 1 [START_REF] Morel | Frequency tuning of piezoelectric energy harvesters thanks to a short-circuit synchronous electric charge extraction[END_REF]), the maximum pre-defined frequency band is given by [1, I1 + 𝑘 ! " ]. With a SCSECE interface (min(𝜀 2 ) = 1/2 -2/𝜋, max(𝜀 2 ) = 1.6366 [START_REF] Morel | Resistive and reactive loads' influences on highly coupled piezoelectric generators for wideband vibrations energy harvesting[END_REF]), the maximum pre-defined frequency band is given by

[I1 + 𝑘 ! " (1/2 -2/𝜋), I1 + 1.6366 𝑘 ! " ]. ∀𝑛 ; ∈ [1, 𝑁 ; ], 𝜀 -,& & = ƒ L1 -Ω %,&; " + 𝑘 ! " 𝜀 2,& & M " + Ω %,& & " 𝑄 𝑘 ! " (14) max(𝜀 -) ≥ … (Ω !35 -Ω !7& ) " 4 + Ω !35 " 𝑄 𝑘 ! " (15) 𝑘 ! " ≥ 1 𝑄 max(𝜀 -) = 𝑘 !,; " (16)
D Ω !7& ≥ I1 + 𝑘 ! " min (𝜀 2 ) Ω !35 ≤ I1 + 𝑘 ! " max (𝜀 2 ) (17)

Number of couples: a trade-off between power and convergence time

Figure 5 shows the values of these optimal couples in the plane (𝜀 -, 𝜀 2 ) for various number of couples 𝑁 ; , as well as the associated power-frequency responses. The black dotted line (𝑁 ; → ∞) corresponds to the case where the number of couples is infinite, and matches the electromechanical impedance matching conditions given by [START_REF] Cottone | Nonlinear Energy Harvesting[END_REF] (which is, indeed, the solution to the continuous limit-case of the optimality problem studied in this paper). From Fig. 5, we can verify that a greater number of couples 𝑁 ; leads to a greater 𝑃 4&G,!7& . If 𝑁 ; is sufficiently large, 𝑃 4&G,!7& gets close to 1 (corresponding to an extraction of the power limit 𝑃 67! on the whole frequency band). Interestingly, the values of the optimal couples L𝜀 -,& & , 𝜀 2,& & M greatly vary with the number of these couples 𝑁 ;3 , meaning that the optimal couples cannot be directly determined from the impedance matching theory and are not equal to (𝜀 - : $%# , 𝜀 2 : $%# ) (7).

Figure 6.a shows the evolution of 𝑃 4&G,!7& as a function of the number of couples 𝑁 ; , for various piezoelectric energy harvesters (𝑘 ! " 𝑄), with a fixed 𝑄 = 30. In in each case, Ω !7& has been fixed to 1 and Ω !35 has been fixed to I1 + 𝑘 ! " = √1.3, which correspond to the normalized short-circuit and open-circuit resonant frequencies of the PEH [START_REF] Erturk | Issues in mathematical modeling of piezoelectric energy harvesters[END_REF][START_REF] Liao | Maximum power, optimal load, and impedance analysis of piezoelectric vibration energy harvesters[END_REF]. Note that the results shown in Fig. 6.a remain sensibly the same with other values of 𝑄, as long as the products 𝑘 ! " 𝑄 remain unchanged. With a greater 𝑘 ! " 𝑄, it becomes necessary to increase the number of couples 𝑁 ; in order to maintain a large 𝑃 4&G,!7& . Indeed, if 𝑘 ! " is increased, the frequency range

[Ω !7& , Ω !35 ]
where the PEH dynamics can be tuned becomes larger [START_REF] Morel | Resistive and reactive loads' influences on highly coupled piezoelectric generators for wideband vibrations energy harvesting[END_REF], leading to a decrease of 𝑃 4&G,!7& for a given 𝑁 ; . If 𝑄 is increased, the maximum harvested power increases (because 𝑃 67! grows linearly with 𝑄 [START_REF] Morel | A comparative study of electrical interfaces for tunable piezoelectric vibration energy harvesting[END_REF]), but the bandwidth associated with each power-frequency response become smaller (∝ 1/𝑄) which also leads to a decrease of 𝑃 4&G,!7& for a given 𝑁 ;4 . However, a large 𝑁 ; leads to a slower algorithm convergence. The algorithm convergence time can be majored by considering an algorithm that tests every couple, maintain each of them during a specific duration, then select the optimal one (such as the one illustrated in Fig. 7). In this case, the total time 𝑡 ;8&G to find the optimal couple is given by the product of the number of couples with the time taken to converge to steadystate for a single couple. By considering that the PEH behaves as a second-order system, its settling time periodic regime is given by 𝑡 H4+ = -

0 K [-] N c =1 N c =2 N c =3 N c =4 N c =6 N c =8 N c =10 N c ! 1
" P QRS+86 T*L ! ,-. U V '
, with 𝑡𝑜𝑙 being the tolerance fraction that is the maximum percentage of error between the final value and the value obtained after a time 𝑡 H4+ , as depicted in Fig. 7. Therefore, 𝑡 ;8&G can be evaluated as [START_REF] Morel | Frequency tuning of piezoelectric energy harvesters thanks to a short-circuit synchronous electric charge extraction[END_REF]:

The number of periods before convergence, 𝑁 #4978W , can thus be evaluated as: . As shown in Fig. 6.b, a larger number of couples as well as a larger 𝑄 tend to increase the time taken by the algorithm before its convergence. Therefore, there exists a tradeoff between the time taken by the algorithm to converge and the minimum value of the power envelope, 𝑃 4&G,!7& . Depending on the application requirements and on the PEH characteristics, 𝑁 ; should be chosen in order to obtain a sufficient power in a minimum convergence time. 𝑡 ;8&G = -

2 𝑄 𝑁 ; ln "𝑡𝑜𝑙…1 - 1 4𝑄 " • 𝜔 % ( 18
)
𝑁 #4978W = - 2 𝑄 𝑁 ; Ω ln "𝑡𝑜𝑙 …1 - 1 4𝑄 " • 2𝜋 (19)

Implementing the DOCT with a tunable electrical interface

Equations ( 11) and ( 14) give the optimal values of the electrical damping 𝜀 -and stiffness 𝜀 2 for a given number of couples 𝑁 ; . To find the optimal couples of electrical parameters (𝜓 * , 𝜓 " ), the interface-dependent relation between (𝜀 -, 𝜀 2 ) and (𝜓 * , 𝜓 " ) has to be known. Such relations are given in Table 2 for some tunable interfaces that have been proposed in the recent literature. Note that for each interface, the electrically-controllable parameters can be of different nature: for the SCSECE, 𝜓 * represents the phase-shift between the displacement extremum and the harvesting event, and 𝜓 " is the angular duration of a short-circuit sequence. For the FTSECE and the FT-OSECE, 𝜓 * represents the voltage inversion ratio, and 𝜓 " represents the phase-shift between the displacement extremum and the harvesting event. In the case of the RC load, 𝜓 * represents the normalized resistive load 𝜓 * = 𝑅 683W 𝐶 # 𝜔 % and 𝜓 " represents the normalized capacitive load 𝜓 " = 𝐶 # /(𝐶 # + 𝐶 683W ).

Table 2 -Relations between the electrically-controllable parameters (𝝍 𝟏 , 𝝍 𝟐 ) and the electrical damping and stiffness (𝜺 𝑫 , 𝜺 𝑲 ) for some two-tuning interfaces in the literature.

Strategy 𝜀 - 𝜀 2 Ref.
Tunable RC load

(𝝍 𝟏 Ω) " 1 -𝝍 𝟐 '1 + ) 𝝍 $ Ω 1 -𝝍 𝟐 * " + %$ 𝝍 𝟏 Ω '1 + ) 𝝍 $ Ω 1 -𝝍 𝟐 * " + %$ [13, 14] FTSECE 4 𝜋 1 -𝝍 𝟏 1 + 𝝍 𝟏 cos " (𝝍 𝟐 ) 1 + 2 𝜋 1 -𝝍 𝟏 1 + 𝝍 𝟏 sin(2𝝍 𝟐 ) [17] SCSECE [cos(𝝍 𝟏 ) + cos(𝝍 𝟏 + 𝝍 𝟐 )] " 𝜋 1 - 𝝍 𝟐 𝜋 + sin(2𝝍 𝟏 + 2𝝍 𝟐 ) 2𝜋 + sin (2𝝍 𝟏 ) 2𝜋 + 2 cos(𝝍 𝟏 + 𝝍 𝟐 )sin (𝝍 𝟏 ) 𝜋 [18] FT-OSECE 4 𝜋 1 -𝝍 𝟏 1 + 𝝍 𝟏 cos " (𝝍 𝟐 ) 1 + 2 𝜋 1 -𝝍 𝟏 1 + 𝝍 𝟏 sin(2𝝍 𝟐 ) [33]
Some other tunable interfaces in the literature such as the PS-SSHI [START_REF] Hsieh | Improving the scavenged power of nonlinear piezoelectric energy harvesting interface at off-resonance by introducing switching delay[END_REF][START_REF] Zhao | Phase-variable control of parallel synchronized triple bias-flips interface circuit towards broadband piezoelectric energy harvesting[END_REF] or the SC-SEH [START_REF] Morel | Active AC/DC control for wideband piezoelectric energy harvesting[END_REF] have not been included in Table 2 because their expressions of 𝜀 -and 𝜀 2 are too long. However, such expressions can be found in the associated references and can be used to apply the DOCT methodology to these interfaces as well.

From the expressions given in Table 2 and equations ( 11) and ( 14), it becomes possible to compute the optimal couples of electrically-controllable parameters (𝜓 * , 𝜓 " ) for any tunable interfaces.

Time Instantaneous Power

Time during which the optimal couple is held

! !"#$ ! !"#$ Time Instantaneous Power ! %&' ! %&' ! %&' ! %&' (" !,# , " $,# ) (" !,% , " $,% ) (" !,& , " $,& ) (" !,' , " $,' )
Power measurement with each couple

! !"#$ " ( , " ) = (" (,+ , " ),+ ) Optimal couple (Re-)Start the algortihm (' = 1)
Set the ' ', couple " (,-, " ),- Measure the power with the ' ', couple " (,-, " ),- Wait for ! !"#$ ' = ' + 1

Hold the optimal couple

Every couple has been tested?

Yes No

Experimental validation

In order to validate the proposed methodology, we designed a custom strongly-coupled PEH which is shown in Fig. 8. This harvester is made of a steel cantilever beam (total volume: 60 mm 3 ), covered by two plates of PZN-PT single crystals. Two long profile masses (total volume: 375 mm 3 ) in steel are fixed to the harvester in order to decrease its resonance frequency and increase its power density [START_REF] Morel | Frequency tuning of piezoelectric energy harvesters thanks to a short-circuit synchronous electric charge extraction[END_REF]. The dimensions of this prototype are indicated in Fig. 8 and its measured characteristics are summarized in Table 3.

This custom harvester is fixed on an electromagnetic shaker powered by a power amplifier and vibrating at an acceleration of amplitude 𝛾 ! = 0.3𝑚. 𝑠 L" .The displacement and speed of the inertial mass are sensed with a laser Doppler vibrometer, and the base acceleration is sensed using an accelerometer. The electrodes of the PEH are connected to a circuit on a breadboard that emulates the tunable SCSECE strategy. The whole experimental setup, illustrated in Fig. 9, is connected to a dSPACE control system that allows the real-time adjustment of the vibration frequency and of the tunable parameters (𝜓 * , 𝜓 " ) of the SCSECE interface, and store the voltage, displacement and acceleration waveforms in order to compute, for each set of parameters, the extracted power. Note that the proposed experimental setup allows to verify the validity of the DOCT methodology, but does not constitute a self-powered low-power implementation of a DOCT algorithm. However, we think that such low-power implementation is possible with a dedicated ASIC, as proved by recent realizations in the literature [START_REF] Morel | 2 self-tunable phaseshifted SECE piezoelectric energy-harvesting ic with a 30nW mppt achieving 446% energy-bandwidth improvement and 94% efficiency[END_REF][START_REF] Cai | A piezoelectric energy-harvesting interface circuit with fully autonomous conjugate impedance matching, 156% extended bandwidth, and 0.38 µW power consumption[END_REF]. The optimal (𝜓 * , 𝜓 " ) are numerically found from equations ( 11) and ( 14) and from the electrical damping and stiffness expressions given in Table 2. These optimal (𝜓 * , 𝜓 " ) (corresponding to the phase-shift between the displacement extremum and the harvesting event, and the angular duration of a short-circuit sequence, respectively) are shown in Table 3. We have experimentally fixed the values of 𝜓 * and 𝜓 " to the values indicated in Table 3. Then, for 38 vibration frequencies chosen between 95 Hz (Ω = 1.00) and 116 Hz (Ω = 1.22), we have acquired the voltage waveforms. Thereafter, the extracted power from the piezoelectric energy harvester is determined by calculating 𝑃 45+ = 𝐶 # 𝑓 𝑣 ,39G " [START_REF] Morel | Frequency tuning of piezoelectric energy harvesters thanks to a short-circuit synchronous electric charge extraction[END_REF] with 𝑓 the vibration frequency and 𝑣 ,39G the piezoelectric voltage right before the energy extraction. The measured extracted powers as a function of the vibration frequency (normalized by the limit power 𝑃 67! ) for a single couple of parameters (𝑁 ; = 1), three couples of parameters (𝑁 ; = 3), and five couples of parameters (𝑁 ; = 5) are shown in Fig. 10.

Discussions

As shown in Fig. 10, the theoretical predictions are in good agreement with the experimental data that have been measured with the automated testbench. The differences between experimental data and theoretical predictions may come from some mechanical and piezoelectric nonlinearities that have not been considered in the linear model given by (1), and from the influence of the higher resonant modes of the cantilever beam that have been neglected in our analysis. As shown in Fig. 10, the minimum power envelope 𝑃 !7&,4&G is equal to 0.35 for a single couple of parameters, 0.63 for three parameters, and 0.73 for five parameters. This shows, for instance, that more than 63% of the power limit 𝑃 67! can be harvested on a frequency band as large as 22% of the PEH short-circuit resonant frequency with only three well-chosen couples of electrical parameters (that have been determined with the DOCT methodology). While the proposed approach requires the characterization of the PEH (because it is necessary to know its 𝑘 ! " and 𝑄 to find the optimal couples), Fig. 11 shows that an error on these values (or a derive of these values because of aging or temperature variations) does not impact much the power envelope obtained with the optimal couples. Indeed, an error of 20% on 𝑄 decreases the minimum of the power envelope by less than 0.01𝑃 67! . An error of 5% on 𝑘 ! " leads to a decrease of about 0.05𝑃 67! .

Control system (dSPACE)

Theoretical results

Experimental results

! ! ! " ! ! ! " ! ! ! " ! ! ! " ! ! ! " ! ! ! " ! ! ! " ! ! ! " ! ! ! "
Note that the proposed methodology, DOCT algorithm, as well as the sensitivity analysis (Fig. 11) could be applied to any electrical interface that exhibits, at least, two tunable parameters. The only condition on the electrical interface is that its electrically-induced damping and stiffness (with a given set of tunable parameters) should meet the optimality criteria given by ( 11) and [START_REF] Liao | Optimal parameters and power characteristics of piezoelectric energy harvesters with an RC circuit[END_REF]. The convergence times of any electrical interface associated with the DOCT algorithm, are, in theory, identical and can be estimated with [START_REF] Hsieh | Improving the scavenged power of nonlinear piezoelectric energy harvesting interface at off-resonance by introducing switching delay[END_REF]. As a matter of example, a DOCTbased algorithm (with the five couples of optimal parameters shown in Table 3) has been simulated on Simulink, with a model of the SCSECE interface and a model of the PEH shown in Fig. 8. Figure 12 shows the piezoelectric voltage waveforms when the normalized vibration frequency shifts from Ω = 1.16 to Ω = 1.05. As shown in Fig. 12, because of this frequency shift, the voltage (and the power) drastically decreases when 𝑡 = 8𝑠. Around 𝑡 = 10𝑠, all the five optimal couples (𝜓 * , 𝜓 " ) are successively tested (in approximately 1 second, as estimated by [START_REF] Hsieh | Improving the scavenged power of nonlinear piezoelectric energy harvesting interface at off-resonance by introducing switching delay[END_REF]) and the optimal couple is kept and stored in the system.

) % = 32°) & = 67°) % = 176°) & = 112°) & ) % ) & ) %
Such approach enables the design of low-power and fast algorithms that could determine, depending on the vibration frequency, which couple among the three optimal ones should be chosen. The convergence time of such algorithms could be much smaller than traditional approaches based on multidimensional gradients between a set of hundreds of parameters values. For instance, the time taken to test every couple of parameters and find the optimal one with the PEH shown in Fig. 8 and for a single couple of parameters, three couples of parameters and five couples of parameters, are 0.2, 0.6 and 1 second, respectively. In comparison, the time taken by a gradient algorithm (with only a single tuning parameter) [START_REF] Morel | 2 self-tunable phaseshifted SECE piezoelectric energy-harvesting ic with a 30nW mppt achieving 446% energy-bandwidth improvement and 94% efficiency[END_REF] is around 15 to 20 seconds. Such quick convergence allows the fast adjustment of the resonant frequency, enabling the compensation of aging and temperature effects that might induce resonant frequency shifts. Another notable advantage of DOCT-based algorithms is their intrinsic stability: the optimal couple can be found with an exhaustive comparison of the power obtained with each couple, or with a derivative-free method such as a dichotomous search.

Conclusion

In this paper, we introduce a methodology to choose a few optimal couples of tunable parameters that maximize the harvested power of a given harvester on a large frequency band. The generic expression of the optimal couples has been derived analytically, and is valid for any two-tuning electrical interfaces. We showed that in most case, only three to five couples are enough to scavenge a power relatively close to the the available power on the largest possible frequency band. The proposed methodology has been verified experimentally on a tunable interface, the SCSECE. Combined with a strongly coupled piezoelectric energy harvester, we have been able to harvest more than 70% of the available power on a 20Hz frequency band (for a resonant frequency of 95Hz), with only five couples of parameters. Such methodology paves the way toward simple algorithms combining low-power consumption with fast-convergence.

Figure 2 .

 2 Figure 2. (a) Representation of a linear piezoelectric energy harvester and (b) power chain of the harvesting system.

  ) + 𝐷𝑥(𝑡) + 𝐾𝑥(𝑡) + 𝛼𝑣 #

Figure 3 .

 3 Figure 3. Power envelopes of a PEH with unbalanced couples: (a) 𝑷 𝒆𝒏𝒗,𝒎𝒊𝒏 = 𝟎. 𝟒𝟐 and (b) 𝑷 𝒆𝒏𝒗,𝒎𝒊𝒏 = 𝟎. 𝟔, with 𝛀 𝒎𝒊𝒏 = 𝟏, 𝛀 𝒎𝒂𝒙 = 𝟏. 𝟏𝟒 and a PEH exhibiting 𝒌 𝒎 𝟐 = 𝟎. 𝟑 and 𝑸 = 𝟑𝟎.

Figure 4 .

 4 Figure 4. Optimization of the power envelope in order to maximize 𝑷 𝒆𝒏𝒗,𝒎𝒊𝒏 with three balanced couples, with 𝛀 𝒎𝒊𝒏 = 𝟏, 𝛀 𝒎𝒂𝒙 = 𝟏. 𝟏𝟒 and a PEH exhibiting 𝒌 𝒎 𝟐 = 𝟎. 𝟑 and 𝑸 = 𝟑𝟎.

Figure 5 .

 5 Figure 5. Values of the optimal couples L𝜺 𝑫,𝒏 𝒄 , 𝜺 𝑲,𝒏 𝒄 M in the plane (𝜺 𝑫 , 𝜺 𝑲 ) with various number of couples 𝑵 𝒄 with 𝛀 𝒎𝒊𝒏 = 𝟏, 𝛀 𝒎𝒂𝒙 = 𝟏. 𝟏𝟒 and a PEH exhibiting 𝒌 𝒎 𝟐 = 𝟎. 𝟑 and 𝑸 = 𝟑𝟎.

Figure 6 .

 6 Figure 6. (a) 𝑷 𝒆𝒏𝒗,𝒎𝒊𝒏 as a function of the number of couples (𝜺 𝑫 , 𝜺 𝑲 ) and of the product of the PEH electromechanical coupling and quality factor 𝒌 𝒎 𝟐 𝑸 (with 𝑸 = 𝟑𝟎), and (b) Number of periods before convergence as a function of the number of couples (𝜺 𝑫 , 𝜺 𝑲 ) and of the quality factor of the PEH 𝑸, with 𝒕𝒐𝒍 = 𝟓%.

Figure 6 .

 6 Figure 6.b shows the evolution of the number of periods before convergence, 𝑁 #4978W , as a function of the quality factor of the PEH, 𝑄, and the number of couples 𝑁 ; . As shown in Fig.6.b, a larger number of couples as well as a larger 𝑄 tend to increase the time taken by the algorithm before its convergence. Therefore, there exists a tradeoff between the time taken by the algorithm to converge and the minimum value of the power envelope, 𝑃 4&G,!7& . Depending on the application requirements and on the PEH characteristics, 𝑁 ; should be chosen in order to obtain a sufficient power in a minimum convergence time.

  Number of periods before convergence [-]Normalized minimum power of the enveloppe[-]

Figure 7 .

 7 Figure 7. Illustration of the algorithm convergence after evaluating the harvested power for each couple of parameters (with 𝑵 𝒄 = 𝟒)

Figure 8 .

 8 Figure 8. Custom strongly coupled piezoelectric energy harvester

Beam, 40×5×0.3 mm 3 3 SupportFigure 9 .

 339 Figure 9. Experimental testbench that has been used in order to evaluate the proposed approach

Figure 10 .

 10 Figure 10. Optimal couples for 𝛀 ∈ [𝟏, 𝟏. 𝟐𝟐] for the SCSECE strategy and with (a) 𝑵 𝒄 = 𝟏, (b) 𝑵 𝒄 = 𝟑 and (c) 𝑵 𝒄 = 𝟓. The lines are the predicted power from (6) and the expressions of the electrical damping and stiffness given in Table 2. The diamonds correspond to the measured power on the experimental PEH shown in Fig.8 and using the testbench shown in Fig.9.

Figure 11 .

 11 Figure 11. Simulation of the variation of the power envelopes with an error of ±𝟐𝟎% on 𝑸 and an error of ±𝟓% on 𝒌 𝒎 𝟐 , with 𝒌 𝒎 𝟐 = 𝟎. 𝟑, 𝑸 = 𝟑𝟎 and 𝑵 𝒄 = 𝟑.

Figure 12 .

 12 Figure 12. Simulation of a DOCT-based algorithm convergence with the PEH shown in Fig.8 and the SCSECE interface (𝒌 𝒎 𝟐 = 𝟎. 𝟒𝟖, 𝑸 = 𝟐𝟎 and 𝑵 𝒄 = 𝟓).

1

 1

  all optimal couples ) % , ) &

Table 1 -Adding tunability to the electrical interface to maximize the harvested power

 1

	Number	Examples	Type of matching	Reach power	Power-frequency response	References
	of tuning	of		limit?
		electrical
		interfaces
						[22]
						[23]
	No	SECE,	No matching	No (except
	tuning	OSECE,		for a given
		etc.		𝑘 ! " 𝑄)
						[24]
		Tunable				[25]
	Single	SECE,	Resistive-Reactive	Yes (for one		[26]
	tuning	Phase-shift	dependent	or two
		SECE, N-	matching	vibration
		SECE, etc.		frequencies)
		Tunable				[15]
	Multiple	cap. bank,	Electromechanical	Yes (for a full		[17]
	tuning	FTSECE,	impedance	range of		[18]
		SCSECE,	matching	vibration		[19]
		PS-SSHI		frequencies if
				𝑘 ! " 𝑄 is large
				enough)

Figure 1. (a) Multi-parameter optimization by mean of a multi-dimensional algorithm and (b) proposed approach based on couples of parameter optimization by mean of a single mono-dimensional algorithm.

				Power	Power limit
					Frequency
				Power	Power limit
					Max. power
					(single tuning)
					Frequency
				Power	Power limit
					Max. power
					(multi tuning)
					Frequency
				First parameter ! !
	a.	Harvested power	Multidimensional algorithm
		information
	Piezoelectric harvester	Two-tuning			MPP
		interfaces	! ! adjustment
		Tuning parameters: ! ! and ! "	! " adjustment		Second parameter ! "
		Harvested power information	Monodimensional algorithm	Power	MPP
	Piezoelectric harvester	Two-tuning	Optimal couple (! ! , ! " )
		interfaces	adjustment
		Tuning parameters:
					Parameter couple ! ! , ! "

  5 in pseud-

	1
	0.9
	0.8
	0.7
	0.6
	0.5
	0.4
	0.3
	0.2
	0.1
	0.1 0	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5	0.55
					0 D [-]

Table 3 -Characteristics of the custom PMN-PT piezoelectric energy harvester

 3

	Parameter	Identified value	Unit
	𝑘 ! "	0.48	-
	𝜔 %	595	rad/s
	𝑄	20	-
	𝐶 #	1.22	nF
	𝑃 QYZ (for 𝛾 ! = 0.3𝑚. 𝑠 L" )	6	µW

Table 3 -

 3 Optimal couples of electrical parameters (𝝍 𝟏 , 𝝍 𝟐 ) for the SCSECE

		𝑁 ; = 1		𝑁 ; = 3				𝑁 ; = 5
		𝑛 ; = 1	𝑛 ; = 1 𝑛 ; = 2	𝑛 ; = 3	𝑛 ; = 1	𝑛 ; = 2	𝑛 ; = 3	𝑛 ; = 4	𝑛 ; = 5
	𝜀 -,& &	0.540	0.194	0.207	0.220	0.143	0.149	0.155	0.161	0.167
	𝜀 2,& &	0.476	0.153	0.476	0.820	0.091	0.276	0.476	0.680	0.891
	𝜓 *,8#+	164°	94.6°	112°	33.3°	139°	177°	18.0°	31.5°	42.3°
	𝜓 ",8#+	86.4°	39.6°	4.50°	57.6°	136°	112°	87.3°	66.6°	47.7°

  Variation of the power envelope with an error of ±20% on & Variation of the power envelope with an error of ±5% on ( !

			1
	nor [-]	0.8
	P
		0.6		1.02	1.04	1.06	1.08	1.1	1.12	1.14
			1							"
	nor [-]	0.8
	P
			0.6	1	1.02	1.04	1.06	1.08	1.1	1.12	1.14
							[-]
					0	5	10	15	20
							Time [s]

Indeed, as shown in Table1, a single tuning interface does not allow independent adjustment of the PEH resonant frequency and of the optimal damping, meaning that the electromechanical impedance matching conditions cannot be met for a whole range of vibration frequencies. In the one hand, multiple tuning interfaces allow to adjust independently the real and imaginary part of the load, meaning that the PEH resonant frequency and PEH damping can be optimized independently.

In other terms, !𝜀 #,% ! , 𝜀 &,% ! $ with 𝑁 ' = 2 is not a subset of !𝜀 #,% ! , 𝜀 &,% ! $ with 𝑁 ' =

[START_REF] Williams | Analysis Of A Micro-electric Generator For Microsystems[END_REF] Note that the absolute value of the minimum of the power envelope with a given 𝑁 ' does not change with an increase of 𝑄, but its normalized value 𝑃 (%),*+% relative to the achievable power, 𝑃 ,+* , decreases because 𝑃 ,+* grows linearly with 𝑄.[START_REF] Yildirim | A review on performance enhancement techniques for ambient vibration energy harvesters[END_REF] The settling time is "the time required for the response curve to reach and stay within a range of a certain percentage of the final value", as defined in[START_REF] Teng-Tiow | High performance control[END_REF].


Keywords: Homogenization, Multiscale, Domain decomposition, Radial approximation

A new multiscale computational strategy for the analysis of structures (such as composite structures) described in detail both in space and in time was introduced recently. This strategy is iterative and involves an automatic homogenization procedure in space as well as in time. At each iteration, this procedure requires the resolution of a large number of linear evolution equations, called the "micro" problems, on the micro scale. In this paper, we present a robust approximate resolution technique for these micro problems based on the concept of radial approximation. This very general technique, which leads to the construction of a relevant reduced basis of space functions, is particularly suitable for the analysis of composite structures.

Multiscale computational strategy with time and space homogenization: a radial-type approximation technique for solving micro problems

Anthony Nouy * , Pierre Ladevèze *

INTRODUCTION

The behavior of heterogeneous structures, such as composite structures, often results from complex phenomena taking place on the microscale. Attempts to predict the behavior of such structures have led to the development of increasingly sophisticated models describing the material on a very fine scale compared to the scale of the structure. Therefore, the numerical simulation of these problems requires the development of efficient computational strategies taking into account these multiscale aspects. For the description of local phenomena with short lengths of variation in space as well as in time to be accurate, the calculations must be performed on a fine discretization of the structure and of the time interval being considered. This condition leads to problems with very large numbers of degrees of freedom, whose computation costs with standard calculation codes are generally prohibitive. One of the main objectives in the last few decades has consisted in developing efficient and robust computational strategies suitable for this type of problem. One of these efficient strategies uses the theory of the homogenization of periodic media [START_REF] Sanchez-Palencia | Comportement local et macroscopique d'un type de mileux physiques hétérogènes[END_REF][START_REF] Duvaut | Matériaux élastiques composites à structure périodique, homogénéisation[END_REF][START_REF] Benssoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Sanchez-Palencia | Non homogeneous media and vibration theory[END_REF]. Other developments and the associated computational approaches can be found in [START_REF] Devries | Homogenization and damage for composite structures[END_REF][START_REF] Ghosh | Multiple scale analysis of heterogeneous elastic structures using homogenisation theory and Voronoï cell finite element method[END_REF][START_REF] Zohdi | Hierarchical modeling of heterogeneous bodies[END_REF][START_REF] Fish | Computational plasticity for composite structures based on mathematical homogenization : theory and practice[END_REF][START_REF] Smit | Prediction of the mechanical behaviour of non-linear heterogeneous systems by multilevel finite element modelling[END_REF][START_REF] Feyel | F E 2 multiscale approach for modelling the elasto-viscoplastic behaviour of long fiber SiC/Ti composite materials[END_REF][START_REF] Kouznetsova | Multiscale constitutive modelling of heterogeneous materials with a gradientenhanced computational homogenization scheme[END_REF][START_REF] Feyel | A multilevel finite element method (F E 2 ) to describe the response of highly non-linear structures using generalized continua[END_REF]. These strategies, however, have some limitations. Besides periodicity, these strategies rely on the fundamental assumption that the ratio between the small-scale length and the large-scale length is small. Moreover, the boundary zones require specific treatment because in these zones the material cannot be homogenized.

A new multiscale computational strategy for nonlinear evolution problems was introduced recently. This strategy involves an automatic homogenization technique in space as well as in time [START_REF] Ladevèze | A multiscale computational method with time and space homogenization[END_REF][START_REF] Ladevèze | Multiscale modelling and computational strategies for composites[END_REF][START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF] which is an extension of previous works limited to space [START_REF] Ladevèze | A new micro macro computational strategy for structural analysis[END_REF][START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF]. This strategy, developed in a general framework, makes no a priori assumption on the form of the solution and, therefore, does not have the limitations of standard homogenization techniques. Until now, it has been developed in the framework of small disturbances of (visco)plastic structures under possible contact conditions with or without friction.

The first characteristic of this strategy consists in partitioning the structure into substructures and interfaces, each of these entities having its own variables and equations. The change of scale takes place only at the interfaces, where forces and displacements are split into "macro" and "micro" contributions. The macro quantities are some mean values of the forces and displacements over space and time. The second characteristic of this strategy is the use of the LATIN method [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF], an iterative resolution technique taking into consideration the whole time interval being studied. At each iteration, one must solve a homogenized macro problem defined over the whole space-time domain, and a set of independent micro problems, which are linear evolution problems defined within each substructure, at their boundaries, and over the time interval being studied.

The resolution of a large number of micro problems at each iteration can lead to prohibitive calculation costs if classical incremental techniques are used. Moreover, when dealing with composite structures containing a large number of similar substructures, the operators and right-hand sides of the micro problems may vary only slightly from one resolution to the next, but classical techniques are unable to take advantage of this characteristic.

A possible answer to this problem is an alternative resolution technique based on the concept of generalized radial approximation [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]. This technique, which is often used in the context of the LATIN method, consists in approximating the solution defined over the space-time domain by a sum of radial functions, each of which is the product of a scalar function of the time variable by a function of the space variable. First of all, this approach reduces the computation cost drastically by transforming an evolution problem defined over the spacetime domain into a series of uncoupled problems defined either over the space domain or over the time domain. It also reduces the storage cost for the variables defined over the space-time domain. Last but not least, it enables the construction of a relevant reduced basis of space functions which turns out to be particularly suitable for composite structures containing families of many identical substructures.

Of course, the radial approximation is not defined uniquely. The classical approach is to define this approximation using an orthogonality criterion associated with a primal or dual formulation of the evolution equation [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]. This definition leads to the "best" approximation only for particular cases, and it provides no simple and systematic way to evaluate the quality of this approximation. Still, this classical approach is efficient in the single-scale framework. However, in some special situations, it may lead to non-convergence of the LATIN method. This lack of robustness is particularly notable for the multiscale computational strategy.

In this paper, we introduce a new and extremely robust approach which maintains a completely mixed vision of the evolution problem. This approach is based on a minimization problem and it now provides a simple and relevant error criterion to evaluate the quality of the radial approximation. This new approximation technique leads to significant savings in terms of computational and storage costs without affecting the efficiency and the robustness of the multiscale computational strategy.

THE REFERENCE PROBLEM ON THE MICRO SCALE

Sections 2 to 4 present a brief review of the main aspects of the multiscale computational strategy. Further details can be found in [START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF][START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF].

Description of the problem

Let us consider, assuming small perturbations and a quasi-static isothermal state, the equilibrium of a structure defined in the spatial domain Ω. At each time t of the interval (0, T ) being studied, this structure is subjected to volume forces f d and surface forces F d on a part ∂ 2 Ω of the boundary. On the complementary part ∂ 1 Ω, the displacement u d is prescribed. All the quantities with subscript "d" are known. The displacements, strains and stresses are subject to initial conditions at t = 0.

Constitutive relation model with internal state variables

The inelastic strain ε p is treated independently of the other internal variables, which are denoted X. The conjugate variable of X is Y; thus, the dissipation rate is:

σ : εp -Y • Ẋ def = σ Y • εp -Ẋ Preprint 4
A. NOUY AND P. LADEVÈZE where • designates a scalar product between dual quantities. From the free energy ρΨ(ε e , X), under the usual decoupling assumptions, one obtains the following "normal" formulation of the state equations (see [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]):

σ = ρ ∂Ψ ∂ε e = Kε e and Y = ρ ∂Ψ ∂X = ΛX (2.1)
where Λ and the Hooke's tensor K are linear, constant, symmetric and positive definite operators characterizing the material. The advantage of such a model lies in the simplicity of the corresponding constitutive relation which, except for (2.1), can be written as a differential equation:

εp -Ẋ = B σ Y (2.2)
where B is a monotonous operator characterizing the material.

Decomposition of the structure into substructures and interfaces

The structure is viewed as an assembly of simple components: substructures and interfaces [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]. Each component has its own variables and equations.

A substructure Ω E , E ∈ E, is subjected to the action of its environment (the neighboring interfaces) defined by a force distribution F E and a displacement distribution W E at its boundary ∂Ω E . An interface Γ EE ′ between two substructures E and E ′ transfers the distributions of both the displacement and the force: W E , W E ′ and F E , F E ′ (Figure 1). h,E respectively. We will designate by s E = (ε pE , X E , W E , σ E , Y E , F E ) the set of the variables associated with substructure Ω E and its boundary. The corresponding space is denoted S (0,T ) h,E . In practice, the interfaces and substructures are discretized in space, classically, using finite elements. On the time scale, one uses the discontinuous Galerkin MULTISCALE COMPUTATIONAL STRATEGY WITH TIME AND SPACE HOMOGENIZATION: A RADIAL-TYPE APPROXIMATION TECHNIQUE FOR SOLVING MICRO PROBLEMS 5 method [START_REF] Eriksson | Time discretization of parabolic problems by the discontinuous galerkin method[END_REF]. The time interval (0, T ) being considered is decomposed into T h = {I i = (t i , t i+1 )} N -1 i=0 , t 0 = 0 < t 1 < ... < t N = T . Let us denote P r (I; S) the space of the polynomial functions of degree r defined over I with coefficients in S. Then, let us introduce the space Z r (T h ; S), r ∈ N N , defined by:

Z r (T h ; S) = {v | v |Ii ∈ P ri (I i ; S), ∀i ∈ {0, ..., N -1}}
To simplify the notations, let Z (0,T ) r = Z r (T h ; R). A space S (0,T ) is said to be of type Z (0,T ) r if S (0,T ) = Z r (T h ; S). On the time scale, all spaces will be considered to be of type Z (0,T ) r .

Reformulation of the reference problem

E (0,T ) h,E,ad ⊂ E (0,T ) h,E × W (0,T )
h,E designates the space of the kinematically admissible variables

(ε E , W E ). F (0,T ) h,E,ad ⊂ S (0,T ) h,E × F (0,T )
h,E designates the space of the statically admissible variables (σ E , F E ) (i.e. the variables which verify the equilibrium equations) and F (0,T ) h,E,0 designates the associated vector space.

Definition 2.1. s E = (ε pE , X E , W E , σ E , Y E , F E ) ∈ S (0,T )
h,E is said to be Eadmissible if it verifies:

the static admissibility of (σ E , F E ), -the kinematic admissibility of (ε E , W E ), -the state equations (2.1),

The corresponding space is designated by S (0,T ) h,E,ad . S (0,T ) h,E,0 is the associated vector space.

Then, the reference problem consists in finding s = {s E } E∈E which verifies:

s E E-admissible, ∀E ∈ E, -the constitutive relation describing the state evolution (2.2), -the behavior at the interfaces, of which the boundary conditions on ∂ 1 Ω and ∂ 2 Ω are particular cases, -the initial conditions.

DESCRIPTION OF THE VARIABLES ON THE MACRO SCALE

Definition of macro quantities

The distinction between the micro level and the macro level is made only at the interfaces. For space, the macro scale is defined by the characteristic length of the interfaces. For time, the macro scale is associated with a coarse subdivision

T h ′ = {I M i = (t ′ i , t ′ i+1 )} N ′ -1 i=0 , t ′ 0 = 0 < t ′ 1 < ... < t ′ N ′ = T , of the time interval being studied.
Let us consider the interface Γ EE ′ between two substructures E and E' and choose the spaces in which the macro displacements and forces will be sought, i.e. W (0,T ),M h,EE ′ and F (0,T ),M h,EE ′ . Classically, for space, one takes affine functions over Γ EE ′ . For displacements, this corresponds to translations, rotations and expansions of the interface (4 functions in 2D). For forces, this corresponds to Preprint 6 A. NOUY AND P. LADEVÈZE the resultants, moments and expansions at the interface. Regarding time, one chooses for the grid T h ′ functions of the type Z (0,T ) p , p ∈ N N ′ . The macro parts W M , F M of W , F are defined as the best approximations in the sense of the work's bilinear form (see [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF]). The micro forces and displacements are simply defined by:

F m = F -F M and W m = W -W M
Thus, one ends up with uncoupled macro and micro works at the interface:

Γ EE ′ ×[0,T ] F • Ẇ ds dt = Γ EE ′ ×[0,T ] F M • Ẇ M + F m • Ẇ m ds dt (3.1)
Extended to the set of the interfaces, this decomposition leads to the spaces

W (0,T ),M h and F (0,T ),M h .

Admissibility of the macro quantities

Another important aspect of the multiscale computational strategy consists in the choice of the admissibility conditions for the macro quantities: the macro forces are required to verify the transmission conditions, including the boundary conditions. These conditions remain unchanged under unilateral contact. The corresponding space is designated by F (0,T ),M h,ad .

MULTISCALE COMPUTATIONAL STRATEGY WITH TIME AND SPACE HOMOGENIZATION

Principle

The partial verification a priori of the transmission conditions at the interfaces leads to the following reformulation of the reference problem:

find s = {s E } E∈E , with s E ∈ S (0,T )
h,E , which verifies:

A d -the E-admissibility of s E , E ∈ E -the admissibility of F M -the initial conditions

Γ

the constitutive relation describing the state evolution (2.2) -the behavior at the interfaces the initial conditions Γ constitutes a set of (possibly nonlinear) equations which are local in space and in time. A d is a set of global linear equations. With this decomposition, it is possible to apply the LATIN method [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF], a general computational strategy for time-dependent nonlinear problems which operates globally over the entire time-space domain. Figure 2 shows the scheme for one iteration, which consists of two stages, called the local stage and the linear stage.

Preprint

MULTISCALE COMPUTATIONAL STRATEGY WITH TIME AND SPACE HOMOGENIZATION: A RADIAL-TYPE APPROXIMATION TECHNIQUE FOR

SOLVING MICRO PROBLEMS 7 (I,Y,F) (Ap,X,W) s n+1/2 -- -+ s n+1 s n s ex ) d /

The local stage at iteration n

The problem consists in constructing s n+1/2 ∈ Γ knowing s n ∈ A d . ∆ s = s n+1/2s n must follow a search direction E + ; more precisely, for each substructure Ω E , ∆ s must verify:

H + ∆ σ E ∆ Y E + ∆ ˙ ε pE -∆ ˙ X E = 0 (4.1)
h + ∆ F E -∆ ˙ W E = 0 (4.2)
H + is a positive definite operator which is a classical parameter of the method. h + is another parameter of the method which can be interpreted as a "viscosity" effect at the interface. The problems in the substructures and at the interfaces are independent of one another. Moreover, these problems are local in space and in time and, therefore, they lend themselves to the highest degree of parallel computing.

The linear stage at iteration n

The problem consists in constructing

s n+1 ∈ A d knowing s n+1/2 ∈ Γ. s n+1 - s n+1/2 must follow a search direction E -.
Because of the admissibility of the macro forces, this search direction is global. Let us introduce a Lagrange multiplier W M to guarantee the admissibility of the macro forces in a weak sense (see [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF]).

W M ∈ W (0,T ),M h,0
, the space of the macro displacements which are continuous at the interfaces and zero on ∂ 1 Ω. Let us denote ∆s = s n+1s n . Then, the search direction can be written as follows: for each substructure Ω E ,

H -∆σ E -∆ σ E ∆Y E -∆ Y E - ∆ εpE -∆ ˙ ε pE -(∆ ẊE -∆ ˙ X E ) = 0 (4.3) h -∆F E -∆ F E + ∆ Ẇ E -∆ ˙ W E = ˙ W M E (4.4)
The resolution of the linear stage can be split into two parts: the resolution of a set of micro problems defined in each substructure, and the resolution of a global macro problem defined over the entire space-time domain. The micro problem associated with substructure E is a linear evolution equation which can be expressed as: find ∆s E ∈ S (0,T )

h,E which verifies:

the E-admissibility to zero:

∆s E ∈ S (0,T ) h,E,0
the search directions (4.3) and (4.4).

Since Relation (4.4) is local at boundary ∂Ω E , the micro problems in each substructure are independent of one another. If H -and h -are positive definite operators, the micro problem defined in substructure Ω E and at its boundary has a unique solution such that:

∆s E = s (1) E + s (2) E ( W M E ) (4.5)
where s

E depends on ∆ s E and s

E depends linearly on W M E . In particular, one has:

∆F M E|(0,T ) = F M (1) E|(0,T ) + L E W M E|(0,T ) (4.6)
L E is a linear operator which can be interpreted as a homogenized behavior operator over the space-time substructure Ω E × (0, T ). This operator can be calculated by solving a set of micro problems over Ω E × (0, T ) in which one takes successively for W The macro problem consists in finding ∆F M and W M which verify:

-∆F M ∈ F (0,T ),M h,0
(static admissibility to zero),

-W M ∈ W (0,T ),M h,0
(kinematic admissibility to zero), -the homogenized behavior relation (4.6).

This problem can be formulated in terms of the macro displacements W M .

If the number of macro space-time substructures is large, an approximation technique based on the introduction of a third scale can be used [START_REF] Loiseau | On a multiscale computational strategy for the analysis of composite structures: discretization and performances[END_REF][START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF].

Resolution of the linear stage

To solve the linear stage, one proceeds as follows: first, one solves a series of micro problems defined over Ω E × (0, T ), E ∈ E, in which one considers only the data ∆ s E of the previous stage. This leads to s

E . Then, one solves the macro problem to obtain W M . Finally, in order to calculate s

E , one solves a second series of micro problems with the Lagrange multiplier as the only data.

Remark 4.1. If the search direction operators H -and h -are updated in substructure E, one must also update the homogenized operator L E . The calculation of this operator requires the resolution of a set of micro problems. The approximation technique presented further on can also be applied to this calculation.

Preprint

MULTISCALE COMPUTATIONAL STRATEGY WITH TIME AND SPACE HOMOGENIZATION: A RADIAL-TYPE APPROXIMATION TECHNIQUE FOR SOLVING MICRO PROBLEMS 9

RESOLUTION OF THE MICRO PROBLEMS: THE CLASSICAL RADIAL APPROXIMATION TECHNIQUE

At each iteration of the computational strategy and for each substructure, one must solve two micro problems which represent evolution equations defined over Ω E × (0, T ), E ∈ E. The resolution of these micro problems with standard incremental methods can lead to prohibitive computation costs, especially if the search direction operators H -and h -are time-dependent. Moreover, in the case of composite structures containing large numbers of identical substructures, the operators and right-hand sides of the micro problems vary only slightly from one resolution to another. The classical approaches are unable to take advantage of this particularity. The radial loading approximation, commonly used in the framework of the LATIN method [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF], is a possible answer to this problem. In this section, we present the classical radial approximation technique [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF] based on the dual variational formulation of the micro problem.

Dual variational formulation of the micro problem

Let us write a micro problem defined over Ω E × (0, T ) formally as follows:

Problem 5.1. Find s E ∈ S (0,T )
h,E,0 which verifies

H -σ E -σ E Y E -Y E - εpE -εpE -( ẊE -ẊE ) = 0 (5.1) h -F E -F E + Ẇ E -Ẇ E = 0 (5.2)
where s E = ∆ s E for the first set of micro problems, and s E = (0, 0, W M E , 0, 0, 0) for the second set. To simplify the notations, let us drop the subscript E. Let us also introduce the following notations:

∆ = H -σ Y - εpE -Ẋ and ∆ = h -F + Ẇ
By injecting relations (5.1) and (5.2) into the weak expressions of the kinematic admissibility and of the state equations, we obtain the following dual formulation of the micro problem:

Problem 5.2. Find (σ, F ) ∈ F (0,T ) h,E,0 and Y ∈ Y (0,T ) h,E which verify: ∀(σ ⋆ , F ⋆ ) ∈ F (0,T ) h,E,0 , ∀Y ⋆ ∈ Y (0,T ) h,E , ΩE ×[0,T ] σ ⋆ Y ⋆ • K -1 0 0 Λ -1 σ Ẏ + H -σ Y dx dt + ∂ΩE ×[0,T ] h -F • F ⋆ ds dt = ∂ΩE ×[0,T ] ∆ • F ⋆ ds dt + ΩE ×[0,T ] ∆ • σ ⋆ Y ⋆ dx dt Remark 5.1. Problem 5.
2 is a classical parabolic evolution equation. In fact, the following developments could be applied to any evolution problem of this type.

Preprint 10

A. NOUY AND P. LADEVÈZE

The concept of radial approximation

The classical radial approximation technique consists in searching an approximate solution ( σ, F , Y) to the dual problem 5.2 as a sum of radial functions, which are products of scalar functions of time by functions of space:

( σ, F )(x, t) = m i=1 c i (t)(C i , C i )(x) Y(x, t) = m ′ i=1 d i (t)D i (x)
where ∀i, c i and

d i ∈ Z (0,T ) r , (C i , C i ) ∈ F h,E,0 and D i ∈ Y h,E .
The inelastic strain and other kinematic internal variables can be deduced from Relations (5.1) and (5.2). Of course, this approximation is not unique. Therefore, the following question arises: how can one obtain the "best" radial approximation?

Seeking the "best" approximation

Definition of the "best" approximation

The classical choice consists in injecting into the variational formulation 5.2 the following test functions:

(σ ⋆ , F ⋆ ) = m i=1 c ⋆ i (C i , C i ) + c i (C ⋆ i , C ⋆ i ) Y ⋆ = m ′ i=1 d ⋆ i D i + d i D ⋆ i
Thus, the approximation is defined by an orthogonality criterion.

The problem to be solved and the resolution technique

To simplify the presentation, let us consider only the case with no internal variable, in which one seeks a first-order approximation (i.e. m=1). Let u

= (σ, F ) ∈ F (0,T ) h,E,0 , v = (ε, W ) ∈ E (0,T )
h,E,0 and

< v, u > (0,T ) = ΩE ×[0,T ] ε : σ dx dt+ ∂ΩE ×[0,T ] W • F ds dt
The variational problem can be written formally as follows: find u ∈ F (0,T ) h,E,0 such that

< L(u), u ⋆ > (0,T ) =< f , u ⋆ > (0,T ) ∀u ⋆ ∈ F (0,T ) h,E,0
where L is the linear differential operator defined by L(u) = (K -1 σ+H -σ, hF ). f = (∆, ∆) is the right-hand side. The problem consists in finding (c, U) ∈

MULTISCALE COMPUTATIONAL STRATEGY WITH TIME AND SPACE HOMOGENIZATION: A RADIAL-TYPE APPROXIMATION TECHNIQUE FOR SOLVING MICRO PROBLEMS

11

Z (0,T ) r × F h,E,0 such that ∀(c ⋆ , U ⋆ ) ∈ Z (0,T ) r × F h,E,0 , < L(cU), cU ⋆ > (0,T ) =< f , cU ⋆ > (0,T ) (5.3) < L(cU), c ⋆ U > (0,T ) =< f , c ⋆ U > (0,T ) (5.4) Equation (5.
3) is a "space" problem whose operators depend on time functions. Equation (5.4) is a simple differential equation in time whose parameters depend on the space functions.

The classical technique for building these radial functions is iterative. It consists in calculating successively the space functions (with fixed time functions) and the time functions (with fixed space functions).

Limitations of the classical approach

Frequently, within the framework of the multiscale computational strategy, one observes that the classical algorithm for the construction of the approximation does not converge. Actually, this construction technique is unsuitable for the definition adopted for the approximation. Indeed, the problem defined by (5.3) and (5.4) can be interpreted as a generalized eigenvalue problem. Since operator L is not self-adjoint, this problem can have complex eigenvalues. Therefore, the proposed iterative algorithm, which can be viewed as a simple power method, generally does not converge. In order to circumvent this problem, the definition of the approximation could be extended to the complex domain. A more sophisticated resolution technique for the eigenvalue problem, inspired by Krylov-type techniques, could also be used. This type of technique constitutes an effective way to build a set of radial functions and could be used when an accurate solution of the evolution problem is necessary. However, the classical approach does not solve the true problem which is to find the best approximation. Moreover, it does not offer a simple criterion to evaluate the quality of the approximation. Sometimes, because of the lack of robustness of this approach, the LATIN method fails to converge; this tendency is particularly notable in the context of the multiscale computational strategy.

RESOLUTION OF THE MICRO PROBLEMS: A NEW RADIAL APPROXI-MATION TECHNIQUE

In order to circumvent the problems encountered with the classical radial approximation technique, one idea consists in defining the approximation u through the minimization of the residue:

| f -L(λU) | 2
This approach, based on the symmetrization of the variational problem, leads to a proper definition of the best approximation. It is appropriate for the spatial discretization framework. However, it is not easily transposable to the continuous framework. Indeed, in the case of partial derivative equations, it is not always easy to use the adjoint operator. Moreover, this formulation of the approximation requires a higher degree of regularity of the functions. This tends to make the a posteriori spatial discretization more expensive. In order to apply the radial approximation, we propose to maintain a completely mixed vision of this evolution problem. Then, a natural approach consists in formulating the micro problem as the minimization of the constitutive relation error associated with the search directions. With this new approach,

12

A. NOUY AND P. LADEVÈZE the approximation relates only to the search directions, which are the parameters of the strategy. Compared to the classical method, this approach certainly increases the computational cost, but it also improves the effectiveness and the robustness of the strategy.

The mixed formulation of the micro problem

The search directions can be interpreted as linear constitutive relations. Therefore, one can introduce the concept of constitutive relation error associated with the search directions. Thus, the micro problem 5.1 can be reformulated as the minimization of this error in the space S (0,T )

h,E,0 of the functions which are E-admissible to zero. The micro problem is expressed as follows: Problem 6.1. Find s ∈ S (0,T ) h,E,0 which minimizes the constitutive relation error e 2 CR,E (ss) associated with the search directions, defined by:

e 2 CR,E (s -s) = h -F -F + Ẇ -Ẇ 2 E,m + H -σ -σ Y -Y - εp -εp -( Ẋ -Ẋ) 2 E,M with ∆ 2 E,M = Ω E ×[0,T ] ∆ • M(t)∆ dx dt ∆ 2 E,m = ∂ΩE ×[0,T ] ∆ • m(t)∆ ds dt
where M and m are linear, symmetric, positive definite operators.

A natural definition for the constitutive relation error consists in choosing

M(t) = (1 - t T )H(t) -1 sym m(t) = (1 - t T )h(t) -1 sym

Radial decomposition of the variables

Here, the idea is to define a radial approximation s compatible with the definition of S (0,T )

h,E,0 , i.e.:

( σ, F )(x, t) = m i=1 c i (t)(C i , C i )(x) Y(x, t) = m ′ i=1 d i (t)D i (x) ( ε, W )(x, t) = m ′′ i=1 e i (t)(E i , E i )(x)
where ∀i, (c i , d i , e i ) ∈ Z (0,T ) r

3 and ((C i , C i ), D i , (E i , E i )) ∈ F h,E,0 × Y h,E × E h,E,0 .
Preprint

MULTISCALE COMPUTATIONAL STRATEGY WITH TIME AND SPACE HOMOGENIZATION: A RADIAL-TYPE APPROXIMATION TECHNIQUE FOR SOLVING MICRO PROBLEMS 13

The inelastic strain and the kinematic internal variables are deduced from the state equations:

ε p = m ′′ i=1 e i E i - m i=1 c i K -1 C i X = m ′ i=1 d i Λ -1 D i
This choice implies that the radial approximation s belongs to the space S (0,T )

h,E,0 . Since K and Λ are constant operators, ε p and X also have radial decompositions. Indeed, all the variables, static as well as kinematic, have radial decompositions. Thus, the new approach yields significant savings in terms of storage cost.

Search for the best approximation

Definition of the best approximation

The definition of the best radial functions is simply that they minimize the constitutive relation error e CR,E . Let us recall that with this new approximation technique, the search directions alone, which are the parameters of the algorithm, are approximated. The approximate solution verifies static and kinematic admissibility as well as the state equations, which was not the case of the classical approach. Then, the best first-order approximation (i.e. m = m ′ = m ′′ = 1) is defined by the following problem: Problem 6.2. Find an approximation s, defined by functions ((C, C), D, (E, E)) ∈ 0 and(c, d, e) ∈ Z (0,T ) r 3 , which minimizes the constitutive relation error (E,E)) constitute the solution of Problem 6.2, then they have the following properties:

F h,E,0 × Y h,E × E h,E,
e 2 CR,E ( s -s) = chC + ėE -∆ 2 E,m + K -1 0 0 Λ -1 ċC ḋD + H cC dD - ėE 0 -∆ 2 E,M 6
• the space functions are obtained in terms of (c, d, e) and of the data by expressing the stationarity of e 2 CR,E with respect to the space functions; • the functions (c, d, e) are such that the Rayleigh quotient R(c, d, e) defined by

R(c, d, e) = ∂ΩE ×[0,T ] {chC + ėE} • m ∆ ds dt+ ΩE ×[0,T ] K -1 0 0 Λ -1 ċC ḋD + H cC dD - ėE 0 • M∆ dx dt Preprint 14
A. NOUY AND P. LADEVÈZE is stationary; d,e) Proposition 6.1 constitutes a significant tool for the estimation of the error associated with a set of radial functions.

• e 2 CR,E ( s -s) = e 2 CR,E (s) -R(c, d, e) = ∆ 2 E,M + ∆ 2 E,m -R(c,

Technique of construction of the radial approximation

Table 1 presents the iterative scheme used to build the radial approximation. This scheme consists in carrying out the minimization on the time functions and on the space functions alternatively. The value of the Rayleigh quotient is an indicator of the convergence of the algorithm.

◮ Initialization (c (0) , d (0) , e (0) ) ◮ Loop on k = 0...k max • With the time functions (c (k) , d (k) , e (k) ) fixed, one seeks (C (k+1) , C (k+1) ), D (k+1) , (E (k+1) , E (k+1) )
which minimizes

e 2 CR,E ( s -s)
• Normalization of the space functions

• With the space functions (C (k+1) , C (k+1) ), D (k+1) , (E (k+1) , E (k+1) ) fixed, one seeks (c (k+1) , d (k+1) , e (k+1)
) which minimizes e 2 CR,E ( ss) • Ending criterion based on the convergence of R(c (k+1) , d (k+1) , e (k+1) ) Table 1: Algorithm for searching a set of radial functions The algorithm of Table 1 converges very rapidly. In practice, one limits the number of iterations to k max = 2 or 3. The resolution of the minimization problems with respect to the time and space functions does not present any particular difficulty. However, it can result in higher computational costs than with the classical approach.

-Minimization with respect to the space functions:

In practice, a standard displacement finite element approach is used. The minimization with respect to the space functions is partially dualized, which leads to seeking two finite element displacement fields simultaneously. Therefore, the size of the problem is twice that of a standard finite element calculation. In general, however, the resolution cost is not prohibitive because the sizes of the micro problems are small. Nevertheless, one can introduce an alternative resolution technique which consists in minimizing with respect to (C (k+1) , C (k+1) , D (k+1) ) and (E (k+1) , E (k+1) ) successively rather than over (C (k+1) , C (k+1) , D (k+1) , E (k+1) , E (k+1) ) simultaneously.

-Minimization with respect to the time functions:

The minimization with respect to the time functions leads to a global problem in time. In practice, this problem is assembled. Its size may be larger than that of the spatial problem, especially if the number of time functions being calculated simultaneously and the number of time steps are large. However, this problem leads to a symmetric block-tridiagonal linear system in which a block is associated with all the unknowns on MULTISCALE COMPUTATIONAL STRATEGY WITH TIME AND SPACE HOMOGENIZATION: A RADIAL-TYPE APPROXIMATION TECHNIQUE FOR SOLVING MICRO PROBLEMS 15 an interval of the partition T h . Such a system can be solved relatively inexpensively using a block-Cholesky technique.

Resolution of the micro problems

Relative error criterion

This new approximation technique comes with a natural error criterion which can be used to estimate the quality of the approximation s. This criterion is based on the constitutive relation error e 2 CR,E ( s-s). Let us define the following norm on S h,E :

| s | 2 E = H σ Y 2 E,M + εp -Ẋ 2 E,M + hF 2 E,m + Ẇ 2 E,m
One can define a relative error criterion ξ CR,E based on the norm of the previous iterate s n as follows:

ξ 2 CR,E ( s -s; s n ) = e 2 CR,E ( s -s) | s n | 2 E ( 6.1) 6.4.2

. Reutilization of the space functions: initialization step

The construction of the space functions is by far the most expensive step. Therefore, it can be interesting to store and reuse these functions. The space functions constructed at each iteration can be reused systematically during the subsequent iterations. Moreover, in the case of composite structures, the basis of the space functions is shared by a large number of identical substructures. Therefore, the initialization step of a micro problem on a particular substructure can reuse the space functions defined in other substructures during the previous iterations (or even the current iteration). Then, the resolution of a micro problem begins with an initialization step which consists in seeking a set of time functions associated with an already existing basis of space functions. This initialization consists in minimizing the constitutive relation error over all the time functions while keeping the space functions fixed.

The operators and the right-hand sides of the micro problems vary little from one iteration to the next, and even from one substructure to another. Therefore, this initialization step is often sufficient to obtain a good approximation of the solution of the micro problem. As we will see in the examples, only a few space functions are necessary to solve a problem with a large number of substructures. The initialization step leads to an initial approximate solution s. To find out whether new sets of functions are necessary, one can use the error criterion ξ CR,E defined by (6.1).

Remark 6.1. The solutions of micro problems in identical substructures can be very similar in terms of forces, stresses, strains and, possibly, internal variables. The solutions in terms of displacements can also be similar, except for the rigid body modes. Therefore, it is essential to add to the basis of functions of the initialization step the space functions corresponding to the rigid body modes of the boundary ∂Ω E . In order for the criterion ξ CR,E to be consistent among substructures, the norm of the interface's displacements in the denominator Preprint 16

A. NOUY AND P. LADEVÈZE must be replaced by a norm of {W }, the part of W with zero resultant and moment over ∂Ω E .

Construction of new sets of radial functions

One chooses a criterion ξ c CR for the verification of the search directions. If this criterion is not matched at the end of the initialization step, new sets of functions must be constructed. Usually, these functions are built one by one. When a new approximation s is obtained, if ξ CR,E > ξ c CR , a new set of functions is calculated using the algorithm of Table 1 with updated data. In practice, in order to reduce the calculation cost (and since a coarse approximation of the search directions is generally sufficient), one can limit the number of new sets of radial functions calculated for the resolution of a micro problem to 1 or 2.

EXAMPLE: COMPOSITE GANTRY

Description of the problem

Let us consider the 2D problem of a composite structure containing cracks (unilateral contact with Coulomb's friction characterized by parameter f = 0.2) (Figure 3). The structure is fixed on the bottom and subjected to forces F 1 and F 2 (Figure 4). The in-plane dimensions are 120 mm × 120 mm. We make the plane strain assumption. The structure consists of two types of cells denoted I and II. The materials are viscoelastic with Young's modulus E i , Poisson's ratio ν i and viscosity η i . Their constitutive relations are such that B i = 1 ηi K -1 i .

Type on macro grid T h ′ on the time level.

Description of the resolution procedure

Since the constitutive relation is linear, the natural search directions in the substructures are H + = H -= B. The characteristic length of the cells is L M = 4 mm. At the interfaces, we chose the same scalar search direction for all the interfaces: h

+ = h -= 1 2 ( 1 η1 + 1 T ) LM E1 .
The structure contains two families of cells. The space functions built by the radial approximation technique generated two common bases which were systematically reused for the initialization step of the micro problems. These bases were orthonormalized using the following scalar product:

((C ⋆ , C ⋆ ), (C, C)) F h,E,0 = ΩE C ⋆ : K -1 C dx ((E ⋆ , E ⋆ ), (E, E)) E h,E,0 = Ω E E ⋆ : KE dx Remark 7.1.
A more consistent procedure than a simple orthonormalization could have been used for sorting the space functions.

Remark 7.2. It is possible to solve the linear stage incrementally over the macro intervals (see [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF]). Therefore, the radial approximation consists in seeking a set of radial functions for each macro interval. In this case, the basis of the space functions is shared by all the macro intervals.

As we will see further, the criterion ξ c CR conditions the accuracy of the approximate solution obtained by the strategy. During the first iterations, s n is only a coarse approximation of the solution (error much greater than ten percent). Greater accuracy would be unnecessary for the verification of the search directions and would lead to the construction of too many new sets of radial functions. Therefore, one can use an adaptive technique which consists in increasing the accuracy ξ c CR as the iterations progress. Here, we took ξ c CR = 0.3 for the first 4 iterations. Then, we took ξ c CR = 0.05. If ξ CR,E > ξ c CR after the initialization step, we calculated a new set of radial functions using the algorithm of Table 1 with parameters k max = 3 and an ending criterion of 10 -2 for the convergence of the Rayleigh quotient. To evaluate the convergence of the algorithm, we calculated the distance between the s n and s n+1/2 . Then, we introduced the following convergence criterion:

err 2 n = E∈E | s n+1/2 -s n | 2 E E∈E | s n | 2 E (7.1)
Figure 5 shows the convergence curves of the iterative strategy with and without radial approximation of the micro problems. The algorithm with approximation had the same convergence rate as the reference algorithm during the first iterations, and converged towards an error of 2.5%, equivalent to the criterion ξ c CR chosen. If better accuracy were needed, one could, of course, modify criterion ξ c CR accordingly.

Construction of the radial functions

Figure 6 shows that very few new sets of radial functions were built at each iteration. Moreover, this number decreases as the iterations progress. Figure 7 shows the evolution of the error with respect to the number of radial functions constructed. Figure 8 shows the number of new sets of radial functions constructed per iteration and per substructure. During the first iterations, new functions were built primarily near the boundary and the edges, and around the high-gradient zones (near the cracks). After a few iterations, new sets of radial functions were constructed only near the cracks. These results account for the relevance of the new approximation technique and for the quality of the proposed error criterion.

Cost of the algorithm

Let us recall that at each iteration two micro problems must be solved for each substructure. The use of a classical incremental technique would have required the resolution of 468 × 60 × 2 = 56, 160 "space" problems per iteration. Here, the radial approximation technique required the construction of 72 sets of radial functions to reach an error of 2.5%. Since k max = 3, at most 216 space problems were solved. Thus, the number of space problems which had to be solved, although their size was indeed twice that of a classical finite element calculation on a cell, was divided by 3,000 for the same accuracy. Assembling the space and time problems, however, can be rather expensive, since this requires the computation of many integrals over time and space. To reduce this cost, it could be interesting to perform a radial decomposition of s and of the operators H and h. Indeed, if both the variables and the operators are decomposed radially, the integrals over the space-time domain can be split into products of integrals over space by integrals over time. A simplified construction technique for these radial decompositions can be found in [START_REF] Ladevèze | A method for computing time and space integrals within the framework of the LATIN method[END_REF].

CONCLUSION

A new radial approximation technique for the resolution of the micro problems of the multiscale computational strategy with time and space homogenization introduced in [START_REF] Ladevèze | A multiscale computational method with time and space homogenization[END_REF][START_REF] Ladevèze | Multiscale modelling and computational strategies for composites[END_REF][START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF] has been presented. This new technique does not have the drawbacks of the classical radial approximation technique [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]. Here, the approximation is based on the minimization of the constitutive relation error associated with the search directions of the computational strategy. This new approach comes with a simple and relevant error criterion to evaluate the quality of the approximation. It also leads to the construction of a relevant reduced basis of space functions which can be systematically reused. It is particularly suited to the case of composite structures containing many similar cells. This new technique drastically reduces the calculation and storage costs of the computational strategy without affecting its robustness and effectiveness. Another possible improvement of the computational strategy consists in approximating the macro problem by introducing a third scale [START_REF] Loiseau | On a multiscale computational strategy for the analysis of composite structures: discretization and performances[END_REF][START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF]. Together, these improvements enable the analyst to deal with structures containing very large numbers of cells.

Figure 1 :

 1 Figure 1: Decomposition into substructures and interfaces

Figure 2 :

 2 Figure 2: Scheme for one iteration

Preprint

8 A

 8 . NOUY AND P. LADEVÈZE 4.3.1. The micro problem defined in substructure Ω E and its boundary ∂Ω E , and over (0, T )

ME

  the macro basis functions of W (0,T ),M h,E . 4.3.2. The macro problem defined over Ω × (0, T )

. 3 . 2 .Proposition 6 . 1 .

 3261 Property of the approximationIf (c, d, e) and ((C, C), D,

2 Figure 3 : 3 Figure 4 :

 2334 Figure 3: Description of the problem Each cell is a substructure of the partitioned problem. The type-II cell is homogeneous, with a Type-1 material whose characteristics are E 1 = 50 P a, ν 1 = 0.3 and η 1 = 10s. It was meshed with 138 TRI3 elements. The type-I cells consists of a matrix of Type-1 material with inclusions of Type-2 material whose characteristics are E 2 = 250 P a, ν 2 = 0.2 and η 2 = 1000s. It was meshed with 262 TRI3 elements. Each interface was meshed with 8 elements. The time interval being considered is (0, 10s). Partition T h contains 60 intervals (N = 60) while T h ′ contains 3 intervals (N ′ = 3). The micro spaces are of type Z (0,T ) 0 on the micro grid T h . The macro spaces contain affine functions on the space level (4

Preprint

3 .

 3 Effectiveness of the new method 7.3.1. Convergence of the algorithm

Figure 5 :

 5 Figure 5: Error err n vs iteration with and without radial approximation

Preprint

Figure 6 :Figure 7 :

 67 Figure 6: Number of new sets of radial functions per iteration

Preprint

Figure 8 :

 8 Figure 8: Number of new sets of radial functions per iteration and per substructure

Preprint