Maria Isabel

Rua Taborda

Catherine Elissalde

U-Chan Chung

Mario Maglione

Egon Fernandes

Armaghan Salehian

Onuma Santawitee

Hélène Debéda
email: helene.debeda@ims-bordeaux.fr

Maria Isabel Rua Taborda

Key features in the development of unimorph Stainless Steel cantilever with screen-printed PZT dedicated to energy harvesting applications

Keywords: PZT, thick film, piezoelectric, energy harvester, stainless steel substrate, interfaces

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Due to its high electromechanical coupling coefficient (50% [START_REF] Priya | A Review on Piezoelectric Energy Harvesting: Materials, Methods, and Circuits[END_REF]) Pb(ZrxTi1-x)O3 (PZT) remains preferentially used as actuators, transducers or sensors in a massive industry production Although the hazardous lead content remains a significant drawback, PZT is also an attractive material for micro-scale energy harvesting applications [START_REF] Lei | MEMS-based thick film PZT vibrational energy harvester[END_REF][3] [START_REF] Aktakka | Thinned-PZT on SOI process and design optimization for piezoelectric inertial energy harvesting[END_REF]. Piezoelectric vibration EH based on polymers (PVDF, polyamides or parylene C) and lead-free materials (AlN, ZnO) are also investigated in the literature but the first ones present serious lifetime and aging issues and the second ones suffer from moderate piezoelectric coupling factors [START_REF] Marin | Multi-mechanism vibration harvester combining inductive and piezoelectric mechanisms[END_REF][6] [START_REF] Elfrink | Vibration energy harvesting with aluminum nitride-based piezoelectric devices[END_REF][8][9] [START_REF] Lin | Piezoelectric Nanogenerators Based on Zinc Oxide Nanowire Arrays[END_REF]. Many designs of piezoelectric vibration EH use a cantilever as the mechanical converter. This simple cantilever structure consists of a supporting body (passive substrate), a piezoelectric layer as the active material and a proof mass for both tuning the resonance frequency and increasing strain in the active material. Thin piezoelectric films (<1 μm) are less suitable when large power density is required due to their weak electromechanical coupling. Vibration EH based on screen-printed piezoelectric thick films (1 to 100µm) lead to better energy densities than thin layers and overcome the assembly issues inherent to bulk ceramics [START_REF] Hindrichsen | Advantages of PZT thick film for MEMS sensors[END_REF]. Early studies devoted to vibration EH using piezoelectric thick or thin films are mostly based on silicon [START_REF] Lei | MEMS-based thick film PZT vibrational energy harvester[END_REF][4][12] [START_REF] Kulkarni | Vibration based electromagnetic micropower generator on silicon[END_REF]. Although high normalized power densities are obtained, their brittleness makes them unfeasible to exploit ubiquitous ambient vibrations. Thus, metal-based substrates have received significant interest, especially stainless steel substrates over silicon ones, due to better mechanical strength, a lower Young's Modulus and the possibility to be pre-stressed in order to avoid tensile cracks of ceramic piezoelectric materials [3][14-21].

Recently, PZT layers supported on metallic substrates and its interface circuitry were developed, in order to collect and regulate the energy for smart grid applications [START_REF] Fernandes | Design, fabrication, and testing of a low frequency MEMS piezoelectromagnetic energy harvester[END_REF] [START_REF] Yang | A Simple Wireless Sensor Node System for Electricity Monitoring Applications: Design, Integration, and Testing with Different Piezoelectric Energy Harvesters[END_REF]. Power output tests results were presented, showing an improved performance compared to other MEMS-based piezo-electromagnetic devices. By integrating a screen-printed PZT layer sandwiched between two gold electrodes and supported on a stainless steel substrate (SS), vibrations resulting from the induced magnetic force led to harvested power. Moreover the novel design of the printed MEMS device, consisting in a centrally-supported modified zig-zag geometry has allowed decreasing the characteristic high resonance frequencies of MEMS structures and maintaining high strains, with a dissipation of approximately 9 μW across a 1 MΩ load, and a low fundamental frequency of 60 Hz, matching the fundamental frequency of the American power grid.

However, energy harvesting performance depends on a large extent on the microstructure, the properties of the active piezoelectric material and the quality of the interfaces [16][18][24] [START_REF] Xu | Screen printed PZT/PZT thick film bimorph MEMS cantilever device for vibration energy harvesting[END_REF], parameters inherently dependent on the fabrication process of the device. In particular, correct densification step of the PZT layer is required to limit the detrimental effect of porosity (leakage current) and to ensure high piezoelectric properties. Another issue to consider is the mechanical clamping effect arising from the substrate [21][24].

The quality of the interfaces between the different layers can be profoundly affected by the residual thermal stresses resulting from the mismatch of thermal expansion coefficients between the materials constituting the structure, leading to potential cracks in the active layer and bending of the structure [26][27]. Here, investigation of different PZT based EH is reported, focusing on interfaces, microstructure, and highlighting key guidelines for screen-printed EH processing. First, the complete EH structure with zig-zag geometry is described, a simple cantilever multilayer structure corresponding to a branch of the zig-zag EH structure is considered, in order to abridge any stainless steel supported cantilever structure. Then screen-printed disks and PZT bulk ceramics are also investigated. Microstructural features, dielectric, piezoelectric properties are compared to highlight the impact of porosity, passive layers, electrodes, and substrate thicknesses on the electromechanical performances of the multilayer energy harvester structure.

Experimental

Electrodes and piezoelectric layers shaping and deposition

Multilayer Thick film structures

The thick film screen-printing technology is selected to fabricate the EH device, as it represents a low-cost and alternative process to classical silicon manufacturing [START_REF] Hindrichsen | Advantages of PZT thick film for MEMS sensors[END_REF]. It requires specific inks, made of powders and organic binder mixed together, that are then transferred onto the substrate through open areas of screen meshes, previously patterned by lithography. The printed films with a thickness in the range 1-100 μm are dried and fired to achieve sintering and ensure adherence of layers. The minimum width and distance between the printed layers is classically 100µm. The semi-automatic screen-printer DEK Horizon used for the layers deposition allows alignments at +/-20µm. Compositions and preparation of the commercial or home-made polymer or inorganic pastes used are detailed in the material section.

The zig-zag EH and simple cantilever structures, both supported on a metallic substrate, and free-standing printed disks are considered.

Cantilever structures

The zig-zag EH microsystem consists of a three-dimensional thick-film structure made of the thick PZT film as the active layer, sandwiched between two gold electrodes (Fig. 1 (e)). The layers are supported on a 254µm thick stainless steel (SS) substrate. The software CleWin5 is used to design the electrodes and the piezoelectric layer following the substrate geometry. Short-circuits are avoided using a gap of 200 µm between the electrodes layers (Fig. 1 (a)). The layers design (Fig. 1 (b,c,d)) is converted separately into a suitable photo-plotter format, then translated into a patterned thick-film screen using a UV photo sensible polymer film. For the gold electrodes and PZT layers, the screen meshes and the polymer film thicknesses used are 325 mesh -15µm and 200 mesh -50 µm, respectively. The simple cantilever consists in the middle branch of the EH made of the same Au/PZT/Au/SS layers stacking as for the zig-zag (Fig. 1, 2 (e)). The successive layer deposition steps are illustrated in Fig. 2. For targeted thicknesses thinner than 40 microns, one layer is directly printed whereas for higher values of thicknesses, snap-off distance is modified and/or two layers are printed

Printed disks

Screen-printed disks are made of the same Au electrodes, and PZT printed layers. Once printed, these electroded PZT layers must be ultimately released from the substrate. This is achieved thanks to a sacrificial layer process classically used in MEMS process. The role of this sacrificial layer is to support the active layers during the densification process, without affecting the sintering process and with the minimum interactions with the active layers. Then, either during or after the thermal treatment, the sacrificial layer must be removed to obtain a free-standing active layer [START_REF] Debéda | Printed Piezoelectric Materials for Vibration-based Damage Detection[END_REF]. Here, the circular shaping of electroded printed PZT layers is aimed to further conduct piezoelectric characterizations on a simple disk, according to IEEE standards [START_REF]IEEE Standards on Piezoelectric Crystals:Measurements of Pieozolelectric Ceramics-1961[END_REF]. A sacrificial layer deposition needs to be first deposited on a ceramic substrate prior to the Au, PZT, and Au layers stacking (Fig. 3). This sacrificial polyester based layer is removed during the sintering process leading to free-standing piezoelectric disks ready for poling process and electrical characterizations. The screen chosen for the gold electrodes is the same as the one used for the cantilever fabrication. A 150µm thick stencil is used to achieve thick printed PZT layers in one step. The screen for the sacrificial layer is a 70µm screen mesh with a 50µm thick photosensitive film.

Ceramics

The PZT powder is prepared using the same paste formula as for the screen-printing process followed by a drying step at 120°C during 8h. The green pellets (13mm diameter) are uniaxially pressed with a pressure of 18.5± 2 MPa and then sintered following a firing cycle detailed in section 2.3. Afterward, Au electrodes (100nm) are sputtered on the ceramic disk surfaces for poling and electrical characterization.

Materials

A stainless steel substrate (SS301) with a Young Modulus of 193GPa [START_REF] Fernandes | Design, fabrication, and testing of a low frequency MEMS piezoelectromagnetic energy harvester[END_REF] is selected according to its compatibility with the high firing temperature of the screen-printed layers (900°C). This substrate is laser cut into the selected geometry. The electrodes and the polyester-based sacrificial layers are printed using commercial pastes from Electroscience Laboratory (ESL), ESL8836 and ESL244T, respectively. The gold paste contains 64wt % of Au and different oxides such as SiO2, Al2O3, CdO, PbO, to ensure the cohesion along with the compacity of the layer and its adherence to the substrate [START_REF] Lucat | Screen-printed varistors: New strategy for high non-linear coefficient[END_REF]. The PZT paste is prepared with a commercial organic binder from ESL (ESLV400) based on ethylcellulose whose volume percentage is adjusted to obtain a suitable viscosity for printing. The optimal viscosity of the ink has been found with a powder/ organic vehicle ratio (i.e. 85/14wt.%)

The viscosity of the different inks is controlled with a Brookfield rheometer in order to adapt the viscosity for screen-printing (less than 50 Pas.s for a shear rate of 100 s -1 ). The Pb(Zrx,Ti1-x)O3 commercial powder (hard navy PZ26 provided by Ferroperm) shows a perovskite-type structure with a predominantly tetragonal phase at room temperature. PZT powder is mixed with three oxides (1%wt CuO, 0.8%wt Li2CO3, and 1.2%wt Bi2O3 noted LBCu) expected to melt near 600°C and then (thus?) improve PZT densification [START_REF] Debéda | One-step firing for electroded PZT thick films applied to MEMS[END_REF].

After a milling step in order to eliminate agglomerates, the grain sizes of CuO, Li2CO3 and Bi2O3 powders are 5µm ±0.03, 19µm ±0.11 and 6µm ±0.03 respectively. The PZT powder is formed of agglomerates (10 µm) made of micrometer or sub-micrometer particles. The as described mixture PZT +3%wt LBCu paste is used to process ceramics, printed disks, and cantilevers.

Thermal treatment and poling step

After printing, all layers (electrodes, PZT, and sacrificial layer) are dried for 20min at 120°C. A slow drying rate 1°C/min until 120°C is preferred for thick layers (>100µm) to avoid cracks. The stacked printed layers are first isostatically pressed 5 min at 40 MPa and 65°C in order to improve the PZT densification [START_REF] Debéda | One-step firing for electroded PZT thick films applied to MEMS[END_REF].

Printed piezoelectric disks and cantilevers, as well as the ceramics, are then sintered 2 hours at 900°C using the same firing for the sake of properties comparison (Fig. 4). The heating rate of 1°C/min up to 400°C allows a gradual elimination of the organic binder and the elimination of the sacrificial layer (if present).

The ramp is then increased to 40°C/min up to 900°C for sintering. The LBCu sintering aid melts at 600°C.

The final PZT thicknesses are approximately 55 µm for the cantilevers, 100µm for the disk and 2-3 mm for ceramics.

Once densified, the samples are polarized. Ag wires of diameter 100µm are glued on the Au electrodes with Ag epoxy paste. Additional ultrasonic wire bonding of 20µm Au wires is led for the zig-zag cantilevers as detailed in previous work [START_REF] Fernandes | Design, fabrication, and testing of a low frequency MEMS piezoelectromagnetic energy harvester[END_REF]. Maximum voltages chosen to avoid dielectric breakdown are respectively of 0.25 kV, 0.35kV and 1kV for cantilevers, disks, and ceramics. This poling step is done under nitrogen and at 280°C (a temperature slightly below the Curie temperature), during 10min and voltage is maintained over temperature cooling.

Characterizations

Thickness and roughness profiles are measured with a contact profilometer KLA Tencor-Alpha Step.

Particle size distributions of the raw powders are evaluated using laser diffraction spectroscopy from wet dispersions (water) after ultrasound (1 minute) using Malvern Mastersize 2000. Phase structure of the PZT powders and ceramics are identified by X-ray powder diffraction (XRD) using a PANalytical X'Pert MPD diffractometer with Cu K radiation of wavelength  = 1.54056 Å. The XRD patterns are recorded at room temperature in the angular range of 2θ = 8-80°. The microstructures of the samples are observed using a scanning electron microscope (SEM, ZEISS EVO50). Cross polishing (Ilion + II 697) at low energy (5keV) is used to prepare proper cross-section surfaces of the EH multilayer for microstructural observations. Fracture surfaces are observed by SEM for ceramics. The SEM images are analyzed using ImageJ to extract an apparent porosity defined as the fractional pore area on the cross-sections. Chemical analysis is performed using Energy Dispersive X-ray Spectroscopy (EDX). Thermal dilatometry analyses (DIL Netzsch 402) under Argon and with a heating rate of 2°C/min until 600°C (maximum temperature of analyzer) are performed on both PZT ceramic (3 mm thick and 6 mm diameter) and four piled up discs (1 mm thick and 9 mm diameter) of stainless steel (SS301). The dielectric characterizations (permittivity ε', dielectric loss tan δ) are performed in the 100 Hz -1 MHz frequency range using an impedance analyzer (HP 4294A) with a home-made temperature controlling system. Impedance measurements are conducted using an Agilent E5061B. The electromechanical and piezoelectric coefficients are then calculated through the resonance and anti-resonance method using the standard IEEE model [START_REF]IEEE Standards on Piezoelectric Crystals:Measurements of Pieozolelectric Ceramics-1961[END_REF].

Results and Discussion

3.1 Microstructure analysis

Densification of PZT

The microstructures of PZT corresponding to the ceramic, the printed disk, and the simple cantilever are compared (Fig. 5). As mentioned in the experimental section, the PZT is mixed with similar content and nature of organic binder and sintering aid (3%wt LBCu). Whatever the design, the firing profile remains the same for all three structures. Highly densified ceramic is obtained at 900°C (relative density 97 %)

confirming efficient densification at such low temperature thanks to the use of the LBCu sintering aid. SEM images on a fracture surface of the ceramic confirm a low level of apparent porosity and reveal a mean grain size of about 3μm with a narrow grain size distribution (Fig. 5 (a)). The relative density is lower for printed structures (Fig. 5 (b,c)). An estimated apparent porosity of 14% and 20% is evaluated, for printed disks and cantilevers respectively. In this case, the grain size distribution is broader, thus a precise estimation of the average grain size cannot be reliable, considering the melted aspect of the grains or agglomerates. In the co-fired printed structures, the sintering shrinkage mismatch between the electrodes and the active PZT layer has to be considered. Differences in shrinkage behavior can generate stresses and defects and lead to constrained sintering compared to free sintering [32][33]. Retardation of PZT sintering kinetics can thus be expected. Moreover, for cantilevers type structures, mechanical clamping from the substrate is an additional parameter that contributes to constrained sintering. The specific preparation for SEM observation of cantilevers can explain the melted aspect of PZT grains (Fig. 5 (b)). In comparison, free surfaces in ceramics allow both more efficient debinding and densification providing a faster kinetic pathway for mass transport and leading to enhanced PZT grain growth.

Study of EH structure and interfaces in cantilever structures

When considering EH cantilevers structures, thermal stresses take place during firing cycle due to the difference of thermal expansion coefficients (CTE) between the different materials constituting the structure. Fig. 6. shows the complete energy harvester after firing with two different PZT layer thicknesses (55 and 100 µm). Bending of the structure is observed for the 100μm final thickness. Indeed, thermal stresses during the firing lead to significant deformation of the structure. This is confirmed by the high difference between the CTE of SS and PZT ceramic highlighted through dilatometry analysis (Fig. 7). In the high temperature range (500-600 °C), an average of CTE values of 16.7 and 4.7 ppm/°C are obtained, respectively.

Hsuen et al. have proposed a mathematical solution of the thermal stress in multilayered films on substrates.

For this purpose, the system's total strain is decomposed into a uniform bending component [START_REF] Hsueh | Thermal stress analyses of multilayered films on substrates and cantilever beams for micro sensors and actuators[END_REF]. The crosssection of our multilayered system is shown in Fig. 8, where the film layers of individual thicknesses (𝑡 1𝐴𝑢 , 𝑡 𝑃𝑍𝑇 , 𝑡 2𝐴𝑢 ) are bonded sequentially to the stainless steel (SS) substrate (thickness ts). The coordinate system is defined in such a way that the interface between the bottom Au layer (noted 1Au) and the SS substrate is located at z = 0, with 𝑧 = ℎ 𝑛 𝑙𝑎𝑦𝑒𝑟 . Given the conditions: (1) The multilayer structure is subjected to a temperature change (∆𝑇), with residual thermal stresses (σ) (eq. 1 and eq. 2), (2) the substrate's thickness is much more higher than the one of the active layer (𝑡 𝑠 ≫ 𝑡 𝑃𝑍𝑇 ), (3) the layers arrangement is asymmetrical and, (4) only the PZT and SS CTE mismatch are considered to contribute to final bending (Au layer effect is overlooked). The curvature 1 𝑟 ⁄ of the structure can be then evaluated using eq. 3.

𝜎 For a PZT thickness of 100µm, an absolute value of 2.69 mm for the theoretical maximum curvature is obtained with a structure concave facing up whereas the experimental curvature is 3.5 mm (Fig. 6 (b)).

Experimental curvature is evaluated at the highest curvature point containing the PZT active layer of the structure. To validate the trend predicted by the analysis, several cantilevers with different PZT thicknesses have been fabricated (35, 55, 60, 70, 80, 100 µm). The maximum deviation between the experimental and calculated values reaches 29% for a PZT thickness of 100 µm (Fig. 9), this high percentage of error can be explained considering that at higher 𝑡 𝑃𝑍𝑇 (80-100 µm) the condition of the analytical solution 𝑡 𝑠 ≫ 𝑡 𝑃𝑍𝑇 is no longer valid. Decreasing PZT thickness leads indeed to both lower curvature values and a better agreement between experimental and calculated data. Thermal stresses can thus be significantly reduced by adjusting the PZT thickness.

Cantilevers made with PZT thickness lower than 60 µm do not show significant bending (around 1 mm only) (Fig. 6 (c)). The selected PZT thickness is thus 55µm.

SEM observations of the cross-section of a non-deformed cantilever (tPZT55µm) allow to probe the good quality of the interfaces in terms of cohesion and composition (Fig. 10 (a)). Good adhesion of the printed layers is observed between the top gold electrode and the PZT layer without delamination (Fig. 10 (b)). The thickness of the top and bottom Au electrodes are 10 and 8µm ±2 µm, respectively. The top Au electrode roughness (300 nm) allows smoothing the PZT layer roughness.

The porosity within the active layer is homogeneously distributed without any significant evidence of densification gradients. Tensile stress due to the substrate effect does not seem to hinder densification close to the interface. However, the formation of an interfacial layer of about 20 µm between the bottom gold electrode and the SS substrate can be evidenced in Fig. 10 (c). This interfacial layer results from interdiffusion between the substrate and the bottom electrode during the sintering process. Bulk diffusion or diffusion across defect paths (dislocations in the metal, the porosity of gold due to thermal decomposition of organics) can explain the cross-diffusion between the two layers.

EDX analysis of the different layers (Fig. 11) reveals that the interface layer between bottom Au electrode and SS substrate contains 45% of the elements constituting the gold ink (Si, Pb, Cu) and 55% of elements of the SS substrate (Fe, Cr). The high level of oxygen confirms the partial oxidation of the substrate in the vicinity of the interface. Furthermore, the PZT layer analysis reveals the presence of Pb (68 wt%), Zr (21.5wt%), Ti (8.4wt%) and Cu and Bi from LBCu (1.6wt%)

The interface layer still ensures cohesion between the electrode and the substrate without the presence of Kirkendall porosity that could be more detrimental. However, such diffusion process at the electrode/substrate interface creates a passive layer whose mechanical properties would impact more globally the properties of the device and has to be accurately controlled.

Dielectric and electromechanical properties

Dielectric measurements at different temperatures and frequencies are first carried out on PZT ceramics in order to determine the Curie temperature (Tc) (Fig. 12). The Curie temperature, determined at the maximum of capacitance, occurs at 330°C as expected for the PZ26 [START_REF]Pz26 ( Navy I ) Hard relaxor type PZT[END_REF]. Considering measurements performed 24 hours after poling at 0.5kV/mm, a relative permittivity of 1200 and dielectric losses lower than 0.14% were measured at room temperature for 1 kHz. The permittivity value is slightly lower than the one reported for commercial PZT ceramic (' 1300sintering at 1200°C). Nevertheless, the differences between the two ceramics in terms of grain size and poling conditions (2-8kV/mm for the commercial ceramic) can explain this difference in permittivity values. Dispersion of the permittivity at low frequencies near Tc is observed but there is no signature of conduction mechanisms activated at high temperature (T>400°C). However, the increase of tan from 380°C reflects the existence of residual space charges.

Capacitance measurements are also carried out for the printed disk. Relative permittivity close to 600 and dielectric losses  2.5% at room temperature and 1 kHz are obtained. Such lower permittivity values compared to PZT ceramics can be partly explained by the lower relative density in the printed disk but cannot be rigorously compared because the applied poling electric field is not the same. However, the level of porosity is also a parameter that determines the maximum value of the field applied during poling.

Impedance measurements are performed for printed disk and simple cantilever (Fig. 13 (a,b)). The obtained resonance frequencies (maximum of B curve) are fr416 kHz and fr356kHz for the printed disk and the cantilever, respectively. By using the resonant frequency fr and anti-resonant frequency fa, effective electromechanical coupling factor keff can be calculated to characterize an arbitrary resonator at any resonance frequency (fundamental or harmonic) thanks to the relation:

𝑘 𝑒𝑓𝑓 = √ (𝑓 𝑎 2 -𝑓 𝑟 2 ) 𝑓 𝑎 2 (eq.4)
Values of keff of 26% and 11% approximatively are obtained respectively for printed disk and cantilevers, respectively. The lowest keff value obtained for the cantilever clearly shows the effect of the passive metallic substrate. Cantilever piezoelectric coefficients have a contribution from d31 depending on substrate properties, which cannot be equivalent with the standard IEEE model [START_REF]IEEE Standards on Piezoelectric Crystals:Measurements of Pieozolelectric Ceramics-1961[END_REF]. The simple conductance measurement is nevertheless relevant as it is directly related to the piezoelectric response of the structure.

Table 2 reports the material characteristics and summarizes the dielectric and piezoelectric properties of the different PZT based ceramics and printed structures which are compared to a commercial ceramic from Ferroperm. Two examples of PZT layers deposited on metallic substrate investigated in literature are also considered. As expected, the bulk ceramics present the highest dielectric permittivity due to higher density.

Compared to the commercial one (sintering at 1200-1300°C), the ceramic sintered in the present work at 900°C presents slightly lower permittivity values, but that remains in the same range of the commercial ones despite a significant lower sintering temperature. Lower dielectric properties are also found for printed PZT in comparison with ceramics due to lower densification, as discussed in section 3.1.1. Dielectric losses are also increased in printed layers (by a factor 10).

Electromechanical coupling and piezoelectric coefficients k31, kp and d31, are deduced using the standard IEEE model [START_REF]IEEE Standards on Piezoelectric Crystals:Measurements of Pieozolelectric Ceramics-1961[END_REF]. Disk samples (ceramic or printed one), show comparable values to the commercial ones (table 2), comforting the efficiency of LBCu to reach high densification at low temperature [START_REF] Lakhmi | Study of screen-printed PZT cantilevers both self-actuated and self-read-out[END_REF]. For ceramics, even if thin electrodes contribute to higher electromechanical factor, the disk geometry (thicker than for the printed disk) probably introduces errors in the IEEE model. Thinner printed disks are more suitable for the IEEE model and also exhibit high piezoelectric and electromechanical coefficients, despite their lower level of densification and permittivity.

For cantilever structures, lower electromechanical coupling and piezoelectric coefficients are obtained, as confirmed by finite element simulation analysis for the zig-zag cantilever [START_REF] Fernandes | Design, fabrication, and testing of a low frequency MEMS piezoelectromagnetic energy harvester[END_REF]. Overall, lower properties are expected for cantilever multilayer structure due to the presence of the passive metallic substrate supporting the active PZT layers. The substrate induces mechanical clamping stresses that limit the active material densification. Lou-Moeller et al. [START_REF] Lou-Moeller | Screenprinted piezoceramic thick films for miniaturised devices[END_REF] also confirmed the limitations of the properties of printed PZT (PZ26) thick film when bounded on silicon and alumina substrates. Showing that sintering aids, inherent porosity and mechanical clamping from the substrate, are the main three contributions leading to lower properties values. Furthermore, referring to literature where very dense thin films are also deposited by aerosol on a cantilever SS substrate [START_REF] Lin | Piezoelectric micro energy harvesters based on stainless-steel substrates[END_REF], electromechanical coupling factor keff of around 10 are also obtained (table 2). In our case, even if our printed PZT layers present more porosity than the aerosol PZT layers, its higher thickness help to improve the electromechanical coefficient of the multilayer cantilever.

More dense and thicker PZT layers clearly improve the performances of the EH as shown in table 2 by the results reported by Colin et al on thinned bulk PZT glued on both sides of a SS cantilever [START_REF] Colin | Highly efficient low-frequency energy harvester using bulk piezoelectric ceramics[END_REF]. With such a bimorph configuration and even with a glue interface between the SS substrate and the PZT, an increase of more than a factor of two is obtained thanks to the optimized properties of the thinned bulk ceramic.

Finally, even if the electromechanical coupling coefficient reflects the energy conversion ability and is a relevant criterion to choose the active material, the figure of merit (FOM) (eq.5) should also be considered to compare different materials for energy harvesting applications [7][36]. Reduced permittivity is a way to increase the FOM, despite lower piezoelectric properties. Consequently, the printed structures exhibiting lower permittivity than bulk ceramics could be a good alternative for MEMS piezoelectric energy harvester. and the bottom gold electrode after sintering. This interface is considered as a passive layer that ensures cohesion between the electrode and the substrate, but this could limit the dynamics of the structure under performance. Bending arising from CTE differences between the different materials constituting the cantilever structures have been analytical and experimentally shown to be dependent on the thickness of the active piezoelectric material. Both modeling and experimental approaches have confirmed the possibility to reduce thermal stresses by decreasing the PZT layer thickness. As a result, the two main lines of optimization rely on the improvement of the densification and in better control of the interfaces quality (diffusion process, thermal stress). The key parameters to improve dielectric and electromechanical performances in complex ceramic EH structure are the grain size and the composition of the PZT layer, the nature and thickness of electrodes and substrate. The elaboration of bimorphs can also be considered as an additional option to cope with interfaces issues.

Figure captions

  Microstructural characteristics and functional properties of PZT based ceramics, printed disk, and cantilevers have been investigated and compared. Efficient densification of bulk ceramics and printed disks have been obtained at temperatures as low as 900°C with the use of LBCu sintering aid. The optimum permittivity values are obtained for bulk ceramics. Lower densification for cantilever structure is due to the mechanical clamping from the SS substrate that generates constrained sintering. As a result, reduced overall properties have been observed. Nevertheless, the piezoelectric effect has been demonstrated and both effective electromechanical coupling factor and piezoelectric coefficient are comparable to those of literature, confirming the excellent potential of such devices for EH MEMS applications. Compared to bulk ceramics, the main limiting features are the interfaces and bending issues. The microstructural analysis of the interfaces in the cantilever has highlighted the formation of an interphase layer between the substrate

Fig. 1 .Fig. 2 .Fig. 7 .Fig. 8 .Fig. 9 .

 12789 Fig. 1. Design of EH multilayer cantilevers (a-d) : (a) zig-zag cantilever with laser micromachined SS substrate with a zoom on the printed Au and PZT layers, (b) bottom Au electrode, (c) PZT layer, and (d) top Au electrode. All dimensions are in mm. Fig. 2. Fabrication screen-printing steps of EH multi-layer (modified zig-zag and simple cantilever): (a) bottom Au electrode printing, (b) active PZT layer printing (c) top Au electrode printing, (d) cross section of cantilevers. Top view with Au/PZT/Au printed layers of (e) zig-zag EH, and (f) simple cantilever structure equivalent to the middle branch of zig-zag cantilever. Dimension in mm. Fig. 3. Schematic view of the screen-printed layers for the disks fabrication. Fig. 4. Sintering cycle used for the densification of the different PZT based materials: ceramics, printed disks, cantilevers and zig-zag structure (sintering: 2h at 900°C). Fig. 5. SEM images of PZT + 3%wt LBCu for (a) fracture section of ceramic, (b) printed cantilever cross polished section, and (c) fracture section of free-standing printed disk. Fig. 6. (a) Modified zig-zag cantilever EH structure after firing, (b) with a PZT layer of 100μm thickness, (c) with a PZT layer of 55 μm thickness. Fig. 7. Thermal expansion coefficient curves of SS substrate (a), and PZT ceramic (b). Fig. 8. Cross section of multilayer structure with coordinate system used in the analyze [32]. Fig. 9. Comparison between calculated absolute values using eq.3 and the measured experimental values for the curvature of multilayer structure analysis Fig. 10. SEM analysis of (a) Au/PZT/Au/SS, (b) top Au electrode/PZT interface and, (c) interface layer between bottom Au electrode and SS substrate. Fig. 11. Energy Dispersive X-ray (EDX) analysis of SS/Au/PZT multilayer structure: chemical cartography of Au/SS interface layer after sintering. Fig. 12. Permittivity (a) and tan δ (b) for PZT ceramic at different frequencies and as a function of temperature (23<T(°C)<450). Fig. 13. Impedance measurements: Conductance (G) and Susceptance (B) of (a) printed disk (b) and printed cantilever.

Source du renvoi introuvable

  Where 𝜎 𝑠 and 𝜎 𝑃𝑍𝑇 are the thermal stresses in the substrate and PZT layer, respectively. 𝐸 ′ = 𝐸 (1 -𝑣 ⁄ )is the biaxial modulus with 𝐸 the Young's Modulus and 𝑣 is the Poisson's ratio of SS substrate or PZT layer individually. 𝛼 𝑃𝑍𝑇 𝑎𝑛𝑑 𝛼 𝑠 are the thermal expansion coefficients of the substrate and PZT layer, respectively. 1/r is the physical representation of the component curvature upcoming from the mismatch between the substrate and the active layer. Using the properties shown in Erreur !

	1 𝑟	=	6𝐸 𝑃𝑍𝑇 ′ 𝑡 𝑃𝑍𝑇 (𝛼 𝑃𝑍𝑇 -𝛼 𝑠 )∆𝑇 𝐸 𝑠 ′ 𝑡 𝑠 2	= ∑	1 𝑟 𝑖	(with 𝑟 𝑖 = PZT or S) (eq. 3)
			𝑠 =	𝐸 𝑠 ′ (3𝑧 + 2𝑡 𝑠 ) 3𝑟	(for -𝑡 𝑠 ≤ 𝑧 ≤ 0 ) (eq.1)
			𝜎 𝑃𝑍𝑇 =	𝐸 𝑠 ′ 𝑡 𝑠 2 6𝑡 𝑃𝑍𝑇 𝑟	(for 𝑡 𝑠 ≫ 𝑡 𝑃𝑍𝑇 ) (eq.2)

., the estimated curvature of the multilayer structure is calculated for different PZT thicknesses, considering a ∆𝑇 = 100°𝐶.

Table captions Table 1 .

 captions1 Elastic properties, average thermal expansion coefficient (500°C-600°C) and thicknesses of PZT and SS. Table2. Summary of properties for the different PZT based structures.

Acknowledgments

The authors would like to thank Dominique Denux (ICMCB -UMR 5026 CNRS-Université Bordeaux) for thermal analysis, Bernard Plano (IMS-UMR5218 CNRS-Université Bordeaux) and Pascale Garreta (Placamat-UMS 3626 CNRS-Université Bordeaux) for the SEM analysis and Xavier Hochart (EXELLIA-TEMEX) for the pressure step. Part of this work was supported by the grants from the Waterloo Institute of Sustainable Energy.


Pauline Seguinot

Anne Julien Ortiz

C Camarasa

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Impact of sequential inoculation with M. pulcherrima on wine properties and influence of the main characteristics of the must

Materials and Methods

Over the past 50 years, the use of selected S. cerevisiae starters to control alcoholic fermentation has broadened, to ensure the completion of fermentation and to avoid the production of undesirable off-compounds. Non-Saccharomyces yeasts, naturally predominant in grape must, are rapidly outcompeted. Although these species were long viewed as spoilage yeasts, their potential for improving the sensory quality of wines is now acknowledged. They possess some specific metabolism pathways and excrete hydrolytic enzymes involved in the release of varietal aromas. Since their use is recent, their behavior during fermentation, in particular the interactions with S. cerevisiae in mixed culture, remain little known.

The contribution of M. pulcherrima in sequential inoculation with S. cerevisiae to wine fermentation was investigated in this study. With this aim, we analyzed the impact of the main components of grape, i.e. concentrations of sugars, nitrogen and lipids, on the production of aroma compounds (higher alcohols, acetate and ethyl esters and thiols) during M. pulcherrima / S. cerevisiae sequential fermentation.

The comparison between sequential inoculation and pure cultures revealed the benefit of using M. pulcherrima in fermentation. Thanks to the Box-Behnken model, the impact and interaction of the parameters were assessed, showing that nitrogen was the most influencing nutrient on the sensory and organoleptic profile of wines. These results are of great interest for the management of M. pulcherrima in sequential inoculation in alcoholic fermentation. They pave the way for further investigations on the metabolic and transcriptomic origins of the phenotypic specificities observed, and the interactions taking part between M. pulcherrima and S. cerevisiae. • Quadratic effect of nitrogen on the release of 3MH and 4MMP

Experimental design Fermentation set-up

• Positive effect of nitrogen on the production of 3MHA 3MHA production using similar mechanisms than the other acetate esters

Response surfaces :

••••

  with central point in triplicate • Determination of the effects of parameters on aroma production 250 mL of synthetic must (Bely et al. 1990) + thiols precursors Samples at 90% of fermentation progress  Aroma analysis by GC-MS Inoculation : • Sequential fermentation : MP at 1.10 7 cells/mL + SC at 5.10 6 cells/mL after 48 hours of fermentation • SC (5.10 6 cells/mL) in pure culture as a control. Yeast strains : Metschnikowia pulcherrima Flavia® (MP) et Saccharomyces cerevisiae Lalvin QA23® (one parameter is fixed, the other two are represented on x and y axis • the studied aroma is represented on the z axis Positive effect of nitrogen • Interaction between lipids and nitrogen :  positive effect of lipids with low nitrogen concentrations  negative effect of lipids with high nitrogen concentrations Regulation of ATF1 by both nitrogen and lipids Positive effect of nitrogen : already observed with S. cerevisiae (Torrea et al. 2011, Garde-Cerdan and Ancin-Azpilicueta, 2008) • Negative effect of lipids Repression of lipids biosynthesis and thus ethyl esters production by exogenous lipids (Chirala et al. 1992) ? Quadratic effect of nitrogen, as for S. cerevisiae (Rollero et al. 2015, Mouret et al. 2014)  Except for phenylethyl alcohol : no significant effect of nitrogen


Jin Xie

Jinlan Xu

Zhenyu Dong

Gang Xu

Chongyang Deng

Bernard Mourrain

Jessica Yongjie

Zhang

Interpolatory Catmull-Clark volumetric subdivision over unstructured hexahedral meshes for modeling and simulation applications

Keywords: interpolatory volumetric subdivision, progressive-iterative approximation, volumetric material modeling, isogeometric analysis

Volumetric modeling is an important topic for material modeling and isogeometric simulation. In this paper, two kinds of interpolatory Catmull-Clark volumetric subdivision approaches over unstructured hexahedral meshes are proposed based on the limit point formula of Catmull-Clark subdivision volume. The basic idea of the first method is to construct a new control lattice, whose limit volume by the CatmullClark subdivision scheme interpolates vertices of the original hexahedral mesh. The new control lattice is derived by the local push-back operation from one CatmullClark subdivision step with modified geometric rules. This interpolating method is simple and efficient, and several shape parameters are involved in adjusting the shape of the limit volume. The second method is based on progressive-iterative approximation using limit point formula. At each iteration step, we progressively modify vertices of an original hexahedral mesh to generate a new control lattice whose limit volume interpolates all vertices in the original hexahedral mesh. The convergence proof of the iterative process is also given. The interpolatory subdivision volume has C 2 -smoothness at the regular region except around extraordinary vertices and edges. Furthermore, the proposed interpolatory volumetric subdivision methods can be used not only for geometry interpolation, but also for material attribute interpolation in the field of volumetric material modeling. The application of the proposed volumetric subdivision approaches on isogeometric analysis is also given with several examples.

Introduction

Volumetric representation describes the interior of an object in addition to its boundaries, which plays an important role in modeling and simulation applications. For example, volumetric modeling is required in advanced manufacturing technologies with heterogeneous materials, such as additive manufacturing of functionally graded material [START_REF] Massarwi | A B-spline based framework for volumetric object modeling[END_REF]]. On the other hand, in three-dimensional isogeometric analysis [START_REF] Lai | Integrating CAD with Abaqus: a practical isogeometric analysis software platform for industrial applications[END_REF]], volumetric modeling plays a key role as a geometric foundation for numerical simulation [START_REF] Zhang | Patient-specific vascular NURBS modeling for isogeometric analysis of blood flow[END_REF], [START_REF] Xu | Isogeometric computation reuse method for complex objects with topology-consistent volumetric parameterization[END_REF]]. Trivariate NURBS solids [Xu et al. (2013a), Xu et al. (2013b), [START_REF] Xu | Exact conversion from Bézier tetrahedra to Bézier hexahedra[END_REF], [START_REF] Xu | Efficient construction of multi-block volumetric spline parameterization by discrete mask method[END_REF], [START_REF] Xu | High-quality construction of analysis-suitable trivariate NURBS solids by reparameterization methods[END_REF], [START_REF] Pan | Volumetric spline parameterization for isogeometric analysis[END_REF], [START_REF] Pan | Low-rank parameterization of planar domains for isogeometric analysis[END_REF], [START_REF] Massarwi | Volumetric untrimming: precise decomposition of trimmed trivariates into tensor products[END_REF]] and trivariate T-spline solids [START_REF] Wang | Converting an unstructured quadrilateral/hexahedral mesh to a rational T-spline[END_REF], [START_REF] Zhang | Solid T-spline construction from boundary representations for genus-zero geometry[END_REF], [START_REF] Wang | Trivariate solid T-spline construction from boundary triangulations with arbitrary genus topology[END_REF], [START_REF] Zhang | Conformal solid T-spline construction from boundary T-spline representations[END_REF], [START_REF] Liu | Volumetric T-spline construction using Boolean operations[END_REF], Wei et al. (2017a)] have been used as modeling and numerical tools in isogeometric modeling and analysis. However, because of their tensorproduct structure, trivariate NURBS and T-spline solids have some limitations on the construction of analysis-suitable volumetric parameterizations from arbitrary complex geometries [START_REF] Wei | Truncated hierarchical tricubic C 0 spline construction on unstructured hexahedral meshes for isogeometric analysis applications[END_REF], [START_REF] Wei | Blended B-spline construction on unstructured quadrilateral and hexahedral meshes with optimal convergence rates in isogeometric analysis[END_REF], [START_REF] Chen | Constructing volumetric parameterization based on directed graph simplification of 1 polycube structure from complex shapes[END_REF]]. Hence, volumetric modeling for geometry with arbitrary topology is a key challenge to be addressed for material modeling and isogeometric analysis.

Subdivision is a simple and efficient technique for surface modeling in computer graphics, which can generate smooth surfaces from an input coarse control mesh with arbitrary topology [START_REF] Catmull | Recursively generated B-spline surfaces on arbitrary topological meshes[END_REF]]. There are some related works on volumetric generalization of subdivision surfaces. However, most of them focus on approximating subdivision [START_REF] Bajaj | A subdivision scheme for hexahedral meshes[END_REF]]. There is very few work on interpolatory volumetric subdivision. In recent work, the authors proposed the limit point formula of Catmull-Clark volumetric subdivision [START_REF] Xu | Exact conversion from Bézier tetrahedra to Bézier hexahedra[END_REF]], which provides a key theoretical foundation for interpolatory Catmull-Clark volumetric subdivision. In this paper, two kinds of interpolatory Catmull-Clark volumetric subdivision approaches over unstructured hexahedral meshes are proposed based on the limit point formula of Catmull-Clark subdivision volume. Main contributions of this paper can be summarized as follows:

• A simple and efficient interpolatory Catmull-Clark volumetric subdivision approach by the local pushback operation is proposed, and several shape parameters are involved in adjusting the limit volumes.

• A new interpolatory Catmull-Clark volumetric subdivision approach based on the progressive-iterative approximation is proposed. The convergence of the iterative process is also proven.

• The applications of the proposed interpolatory subdivision methods for volumetric modeling with heterogeneous materials and isogeometric analysis of heat conduction and linear elasticity problems are also presented.

The rest of the paper is structured as follows. Section 2 describes some related work on interpolatory subdivision surfaces and subdivision volumes. Section 3 reviews the Catmull-Clark volumetric subdivision rules and the corresponding limit point formula. The interpolatory Catmull-Clark volumetric subdivision based on the push-back operation is presented in Section 4. Section 5 describes the interpolatory Catmull-Clark volumetric subdivision based on progressive iterative approximation. Some comparison examples and applications in volumetric material modeling and isogeometric analysis are given in Section 6. Finally, we conclude this paper and outline future work in Section 7.

Related work

Interpolatory subdivision surface. Interpolatory subdivision surface can mainly be split into two categories. In the first one, the subdivision schemes consistently contain the old vertices in the refined mesh after each subdivision operation. Obviously, the limit surface generated in this way will contain all the original vertices. Butterfly scheme [START_REF] Dyn | A butterfly subdivision scheme for surface interpolation with tension control[END_REF]] and Kobbelt's scheme [START_REF] Kobbelt | Interpolatory subdivision on open quadrilateral nets with arbitrary topology[END_REF]] both fall into this category. Li et al. proposed a method for directly deducing new interpolation subdivision masks from the corresponding approximation subdivision masks [START_REF] Kobbelt | Interpolatory subdivision on open quadrilateral nets with arbitrary topology[END_REF]]. A unified interpolatory subdivision scheme is proposed for quadrilateral meshes based on local refinement operations in a way similar to that for approximating schemes generalizing splines of an arbitrarily high-order continuity in [START_REF] Deng | A unified interpolatory subdivision scheme for quadrilateral meshes[END_REF]]. Another idea to generate a subdivision surface interpolating vertices of an input control mesh is to employ approximating subdivision schemes by adjusting vertices of the original control mesh. [START_REF] Halstead | Efficient, fair interpolation using Catmull-Clark surfaces[END_REF]] proposed a method by building a global linear system with some fairness constraints to make the limit surface interpolate vertices of the given control mesh. The push-back method was firstly proposed in [START_REF] Maillot | A unified subdivision scheme for polygonal modeling[END_REF]] to deal with the shrinkage issue of approximating subdivision schemes. [START_REF] Deng | A simple method for interpolating meshes of arbitrary topology by Catmull-Clark surfaces[END_REF]] applied this idea to the Catmull-Clark surface subdivision scheme to derive an interpolation scheme based on the explicit limit point formula. Progressive-iterative approximation (PIA for short) and weighted PIA have been used for constructing interpolatory subdivision surfaces by adjusting the control mesh iteratively [START_REF] Chen | Progressive interpolation based on Catmull-Clark subdivision surfaces[END_REF], [START_REF] Cheng | Loop subdivision surface based progressive interpolation[END_REF], [START_REF] Deng | Weighted progressive interpolation of Loop subdivision surfaces[END_REF]]. However, since the PIA method requires the explicit limit point formula to realize the interpolation, currently there is no related work on interpolatory volumetric subdivision based on PIA. Subdivision surfaces have also been used broadly in design and isogeometric analysis [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF], [START_REF] Pan | A unified method for hybrid subdivision surface design using geometric partial differential equations[END_REF], [START_REF] Pan | Isogeometric analysis based on extended Catmull-Clark subdivision[END_REF]]. Truncation mechanism has been applied to Catmull-Clark subdivision to support local refinement [START_REF] Wei | Truncated hierarchical Catmull-Clark subdivision with local refinement[END_REF], [START_REF] Wei | Extended truncated hierarchical Catmull-Clark subdivision[END_REF]], and a hybrid non-uniform recursive subdivison scheme was developed to support isogeometric analysis with improved convergence rate [START_REF] Li | Hybrid non-uniform recursive subdivision with improved convergence rates[END_REF]]. Subdivision volume. Compared with subdivision surfaces, much less work has been done on volumetric subdivision. The first subdivision solid scheme was proposed by [START_REF] Maccracken | Free-form deformations with lattices of arbitrary topology[END_REF]], which refines a volume mesh into a region of three-dimensional space. Bajaj et al. proposed a volumetric extension of the Catmull-Clark subdivision surface and analyzed the smoothness property [START_REF] Bajaj | A subdivision scheme for hexahedral meshes[END_REF]]. A new volumetric subdivision scheme based on Box splines was developed in [START_REF] Chang | A new solid subdivision scheme based on Box splines[END_REF]]. Chang et al. also presented an interpolating scheme for recursive subdivision of meshes organized around octet-truss structures [START_REF] Chang | An interpolatory subdivision for volumetric models over simplicial complexes[END_REF]]. A natural generalization of the butterfly subdivision from 2D surface to 3D volume was proposed in [START_REF] Mcdonnell | Interpolatory, solid subdivision of unstructured hexahedral meshes[END_REF]] as a novel interpolating volumetric scheme. In the field of scientific visualization, an interactive visualization framework was propsoed by using radial basis functions with subdivision volumes [START_REF] Mcdonnell | Subdivision volume splatting[END_REF]]. Lin et al. proposed a constrained volume iterative fitting algorithm to fill a triangular mesh with an all-hex volume mesh [START_REF] Lin | Quality guaranteed all-hex mesh generation by a constrained volume iterative fitting algorithm[END_REF]]. In that paper, they fit the boundary of the initial all-hex mesh to the given triangular mesh by iteratively adjusting the boundary mesh vertices, and the movements of the boundary mesh vertices are diffused to the inner all-hex mesh vertices in each iteration without the application of limit point formula. As an emerging application, isogeometric analysis based on approximation CatmullClark subdivision volumes were proposed in [START_REF] Burkhart | Iso-geometric finite element analysis based on Catmull-Clark subdivision solids[END_REF]], which can achieve a unified representation of geometric modeling and physical simulation. In this paper, two kinds of interpolatory Catmull-Clark volumetric subdivision approaches over unstructured hexahedral meshes are proposed based on the explicit limit point formula, and their application in volumetric material modeling and isogeometric analysis are also illustrated.

Catmull-Clark volumetric subdivision and limit point formula

Catmull-Clark volumetric subdivision is an volumetric extension of Catmull-Clark surface subdivision. Denote M 0 the initial control lattice and M i the hex-mesh obtained after i steps of Catmull-Clark volumetric subdivision. Let e i k , f i k , c i k be the edge point, the face point and the cell point of a vertex in M i . The Catmull-Clark volumetric subdivision rules for unstructured hexahedral mesh with arbitrary topology can be described as the following five steps:

Step.1 For each cell, insert a cell point (on the centroid) c i k , which is the average of all vertices of this cell.

Step.2 For each face, insert a face point f i k , which is derived from

f i+1 = c i 1 + c i 2 + 2A 4 (1)
where c i 1 and c i 2 are the cell points of two adjacent cells of this face, and A is the centroid of this face.

Step.3 For each edge, insert an edge point, which is derived by

e i+1 = C avg + 2F avg + M 4 (2)
where C avg is the average of cell points of all the incident cells of this edge, F avg is the average of the centroids of all the incident faces of this edge, and M is the midpoint of the edge.

Step.4 For each hexahedral cell, connect the cell points with all the face points, and connect all the face points with all incident edge points in this cell respectively. Then each hexahedral cell is divided into eight hexahedral cells, in which the cell points, the face points and the edge points are new vertices of M i+1 .

Step.5 Update each original vertex according to

v i+1 = C avg + 3F avg + 3E avg + v i 8 (3)
where C avg , F avg , E avg are the averages of the cell points, the face points and the edge points of all adjacent cells, faces and edges, and v i is the original vertex.

The above steps define the new geometry and the topology connectivity of new vertices. With these steps, a series of refined hex-meshes will be produced, which eventually converges to a limit volume called the Catmull-Clark subdivision volume.

Catmull-Clark subdivision surface has an explicit limit formula [START_REF] Halstead | Efficient, fair interpolation using Catmull-Clark surfaces[END_REF]], which plays an important role in theoretical analysis and application of subdivision surface. Recently, the authors proposed an explicit limit point formula of Catmull-Clark subdivision volumes [START_REF] Xu | Exact conversion from Bézier tetrahedra to Bézier hexahedra[END_REF]]. The main result is that the limit point v ∞ i of v 0 on M 0 can be computed in the following explicit form:

v ∞ i = 16(n -2)v 1 i + 4 n j=1 m j e 1 j + 4 3(n-2) j=1 f 1 j + 2(n-2) j=1 c 1 j 30(n -2) + 4 n j=1 m j (4)
where n is the valence of an interior vertex v i , m j is the valence of the j-th adjacent edge of v i , and

V 1 i = {v 1 i , e 1 1 , ..., e 1 n , f 1 1 , ..., f 1 3(n-2) , c 1 1 , ..., c 1 2(n-2) }
is the adjacency of v i . e i , f i , c i are the edge points, face points and cell points respectively after performing Catmull-Clark solid subdivision once.

In the following sections, based on this explicit limit point formula, two kinds of interpolatory Catmull-Clark volumetric subdivision approaches will be introduced.

Interpolatory Catmull-Clark volumetric subdivision with push-back operation

In this paper, we will introduce the interpolatory Catmull-Clark volumetric subdivision with the pushback operation. By using the limit point formula of Catmull-Clark subdivision volumes, from a given hexahedral control lattice, we can construct a new control lattice whose limit subdivision volume interpolates vertices of the original control lattice. Moreover, due to its local property, the interpolatory subdivision results can be obtained efficiently. By changing the shape parameters involved in this method, different subdivision volumes can be constructed to satisfy the interpolation condition.

Suppose that M 1 is the hexahedral control lattice obtained by performing one Catmull-Clark subdivision step on the initial hexahedral control lattice M 0 , V 1 i on M 1 is the corresponding vertex of V 0 i on M 0 . Then we define an increment as

∆ 0 i = V 0 i -V 1 i . (5)
For each edge E ij with end points V 0 i and V 0 j in M 0 , we can define

E = 1 2 (V 0 i + V 0 j ) + λ ij (∆ 0 i + ∆ 0 j ), 0 < λ ij < 1 (6)
as the edge point of this edge. ∆ 0 i and ∆ 0 j are the increments of two ending points respectively and λ ij is a shape parameter, which can be used as a degree of freedom to construct different limit volumes.

For each face F in M 0 , we can define

F = 4 j=1 V 0 j 4 + 2µ F 4 j=1 ∆ 0 j , 0 < µ F < 1 (7)
as the face point of this face, in which ∆ 0 j is the increment for vertices on this face, where µ F is also a shape parameter. Similarly, for each cell C in M 0 , we can define

C = 8 j=1 V 0 j 8 + 2γ C 8 j=1 ∆ 0 j , 0 < γ C < 1 (8)
as the cell point, in which ∆ 0 j is the increment of eight vertices of this cell. Here γ C is a shape parameter influencing the final result.

In order to satisfy the interpolation condition, the original vertices of the control lattice should be set as the limit points, that is,

V ∞ i = V 0 i .
Then we have a new hex-mesh with E , F , C as above and the following vertices:

V = [30(n -2) + 4 n j=1 m j ]V 0 i -4 n j=1 m j E j -4 3(n-2 j=1 F j - 2(n-2) j=1 C j 16(n -2) . ( 9
)
After the above steps, now we can obtain a new hexahedral control lattice whose limit subdivision volume interpolates all the vertices V 0 of M 0 . From the new control lattice, we can obtain a limit volume interpolating the original control points by applying the Catmull-Clark volumetric subdivision rules.

Figure 1 shows an interpolatory volume example with the push-back method. Figure 1(a) is the original control lattice. Figure 1(b) is the control lattice constructed by the push-back operation. Figure 1(c) shows the volumetric subdivision results from the modified control lattice. Figure 1(d) is the interior view of Figure 1(c). We can find that the final subdivision volume interpolates the original control lattice exactly. Figure 2 presents interpolatory volumetric subdivision examples with different shape parameters in the push-back method.

Progressive interpolation using Catmull-Clark volume subdivision

Progressive-iterative approximation is an iterative method with clear geometric meaning, which was first proposed by Qi et al. [START_REF] Qi | The method of numeric polish in curve fitting[END_REF]] and de [START_REF] De Boor | How does Agees smoothing method work[END_REF]]. A detailed review about progressive-iterative approximation was given in [START_REF] Lin | Survey on geometric iterative methods and their applications[END_REF]]. In this paper, based on the limit point formula, we derive a volumetric interpolating scheme based on Catmull-Clark volumetric subdivision by using progressive-iterative approximation.

Interpolation process

Denote by M k the volumetric control lattice from the original volume mesh after k steps of modification. For each vertex v k in M k , we compute its limit position of Catmull-Clark solid subdivision, namely according to the limit point formula. Then the distance d k between the counterpart of v k in M k , namely v 0 , and v k ∞ can be computed as

v k ∞ , (a) (b) (c) (d) (e) (f) (g) (h)
d k = v 0 -v k ∞ . (10)
Then we modify all v k ∈ M k by adding d k and get v k+1 by

v k+1 = v k + d k . ( 11
)
After several iterative steps, the final limit subdivision volume can achieve an interpolation of initial hexahedral control lattice. Figure 3 shows an interpolatory volume example with the PIA method. Figure 3(a)-3(d) show the iterative process of the control lattice, and Figure 3(e)-3(h) illustrate the corresponding iterative process of the limit volumes, in which the blue points are vertices to be interpolated.

Convergence analysis

In this section, we will prove the convergence of the proposed PIA approach for Catmull-Clark subdivision volumes. Since the convergence of PIA on Catmull-Clark subdivision surface has been proved in [START_REF] Chen | Progressive interpolation based on Catmull-Clark subdivision surfaces[END_REF]], we can firstly perform PIA on the boundary surface vertices of the volume mesh and fix them in the following steps. Then, by repeating the process (10)-( 11), we can get a sequence of volume meshes M k (k = 1, 2 . . .). We want to prove that as k tends to infinity, M k converges to the limit subdivision volume, which interpolates all the vertices in the original mesh M 0 .

The key to prove the convergence of M k is to show that when k tends to infinity, d k converges to zeros. For each interior vertex, we have By applying ( 12) to (10), we have

v k ∞ = 16(n -2)v k + 4 j m j e k j + 4 j f k j + j c k j 30(n -2) + 4 j m j . ( 12
) (a) (b) (c) (d) (e) (f) (g) (h)
d k+1 = v 0 - 16(n -2)v k+1 + 4 j m j e k+1 j + 4 j f k+1 j + j c k+1 j 30(n -2) + 4 j m j = v 0 - 16(n -2)v k + 4 j m j e k j + 4 j f k j + j c k j 30(n -2) + 4 j m j - 16(n -2)d k + 4 j m j d k ej + 4 j d k f k j + j d k c k j 30(n -2) + 4 j m j = d k - 16(n -2)d k + 4 j m j d k ej + 4 j d k f k j + j d k c k j 30(n -2) + 4 j m j (13)
in which d k = 0 when it corresponds to a vertex on the surface. Denote l as the number of vertices on the boundary surface, r as the number of all vertices in the given mesh. Equation( 13) can be written in matrix form as

          d k+1 1 d k+1 2 • • • d k+1 l d k+1 l+1 • • • d k+1 r           = (I -M )           d k 1 d k 2 • • • d k l d k l+1 • • • d k r           = (I -M ) k+1           d 0 1 d 0 2 • • • d 0 l d 0 l+1 • • • d 0 r           (14)
where d 0 i = 0 (i = 1, 2, . . . , l), I is the identity matrix of order r, M is a matrix with the following form:

M = 0 0 M 1 M 2 , ( 15
)
M 1 =            4mij 30(n l+1 -2)+4 j mj • • • 4 30(n l+1 -2)+4 j mj • • • 1 30(n l+i -2)+4 j mj • • • 4 30(n l+j -2)+4 j mj • • • • • • . . . . . . . . . . . . 1 30(n l+1 -2)+4 j mj • • • • • • • • • • • •            , (16)
M 2 =                     16(n l+1 -2) 30(n l+1 -2)+4 j mj • • • 4mij 30(n l+1 -2)+4 j mj • • • 4 30(n l+1 -2)+4 j mj • • • 1 30(n l+1 -2)+4 j mj • • • . . . . . . 4mij 30(n l+i -2)+4 j mj 16(n l+i -2) 30(n l+i -2)+4 j mj . . . 4 30(n l+j -2)+4 j mj . . . 1 30(n l+k -2)+4 j mj . . . • • • • • • 16(n l+r -2) 30(n l+r -2)+4 j mj                     (17)
where n i is the valence of the i-th vertex, and m ij is the valence of the i -j edge.

Because d k i = 0 (i = 1, 2, . . . , l) and M ij = 0 (i, j = 1, 2, . . . , l), proving that d k converges to zero. The problem reduces to proving that the absolute value of each eigenvalue of the matrix (I -M ) is less than 1, where M = M 2 . Furthermore, from [START_REF] Chen | Progressive interpolation based on Catmull-Clark subdivision surfaces[END_REF], [START_REF] Cheng | Loop subdivision surface based progressive interpolation[END_REF], [START_REF] Deng | Weighted progressive interpolation of Loop subdivision surfaces[END_REF]], we have the following similar proposition.

Proposition 1. Suppose that λ i (i = 1, 2, . . . , r -l) are the eigenvalues of M = M 2 , if λ i > 0 (i = 1, 2, . . . , r -l), then the absolute value of each eigenvalue of the matrix (I -M ) is less than 1.

In the above proposition, we suppose that the eigenvalues of M are all positive. Hereafter, we have λ i > 0 (i = 1, 2, . . . , r -l). According to [START_REF] Chen | Progressive interpolation based on Catmull-Clark subdivision surfaces[END_REF]], we know that the eigenvalues of the product of two positive definite matrices are positive. Hence, we can first decompose the matrix M into the product of A and B, in which

A =       1 30(n l+1 -2)+4 j mj 0 • • • 0 0 1 30(n l+2 -2)+4 j mj • • • 0 . . . . . . 0 • • • 1 30(nr-2)+4 j mj       , ( 18
) B =                  16(n l+1 -2) • • • 4m ij • • • 4 • • • 1 • • • . . . . . . 4m ij 16 (n l+i -2) . . . 4 . . . . . . 1 . . . • • • • • • 16 (n r -2)                  . ( 19
)
Obviously, A is diagonal and all elements on the diagonal line are positive, thus it is positive definite. Now the only thing that we need to prove is that the matrix B is positive definite.

Proposition 2. Matrix B in ( 19) is positive definite.

Proof. Matrix B is positive definite if and only if the corresponding quadratic form is positive for any non-zero vector X = (x l+1 , x l+2 , . . . , x r ):

f (x l+1 , x l+2 . . . , x r ) = X T BX > 0. ( 20
)
We can see that if vertices v i and v j are adjacent, then b ij = b ji = 4m ij . If they are not adjacent but lie in the same face, then b ij = b ji = 4. If they do not lie on the same face but lie in the same cell, then b ij = b ji = 1. Thus we can get

f (x l+1 , x l+2 . . . , x r ) = r i=l+1 16 (n i -2) x 2 i + eij 8m ij x i x j + fij 8x i x j + cij 2x i x j (21)
where e ij , f ij and c ij range through all edges, all diagonals of each face and all diagonals of each cell inside the volume (we do not include those on the surface).

In order to prove that ( 21) is positive, we construct the following equation. For all edges of the mesh, when x i corresponds to a surface vertex, set it to zero. Then we obtain

eij m ij (x i + x j ) 2 = r l+1 j m j x 2 i + eij 2m ij x i x j (22)
where the sum is taken for all edges e ij of the mesh (including surface edges). For all faces, similarly we have fijpq

(x i + x j + x p + x q ) 2 = r l+1 3 (n i -2) x 2 i + eij 2m ij x i x j + gij 2x i x j ( 23
)
where g ij means a pair of diagonal vertices on a given face. For all cells, we have cijpquxyz

(x i + x j + x p + x q + x u + x x + x y + x z ) 2 = r l+1 2 (n i -2) x 2 i + eij 2m ij x i x j + gij 4x i x j + hij 2x i x j
(24) where h ij means a pair of diagonal vertices on a given cell.

By applying ( 22)-( 24) to ( 21), we have

f (x l+1 , x l+2 . . . , x r ) = cijpquxyz (x i + x j + x p + x q + x u + x x + x y + x z ) 2 + 2 fijpq (x i + x j + x p + x q ) 2 + eij m ij (x i + x j ) 2 + r i=l+1   8 (n i -2) - j m j   x 2 i . ( 25
) (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)
Figure 4: Comparison examples. From left to right: The input hexahedral mesh; the butterfly interpolatory subdivision method; the push-back method; and the PIA method.

For each interior edge, m is equal to the number of adjacent cells. When computing j m j for a vertex, each adjacent cell will be counted three times because it must have three edges adjacent to the vertex. For each vertex, the number of adjacent cells is 2(n i -2). Thus we have

j m j = 6 (n i -2) . (26)
Substituting ( 26) into (25), we can obtain

f (x l+1 , x l+2 . . . , x r ) = cijpquuryz (x i + x j + x p + x q + x u + x x + x y + x z ) 2 + 2 fijpq (x i + x j + x p + x q ) 2 + eij m ij (x i + x j ) 2 + r i=l+1 2(n i -2)x 2 i > 0. (27)
Equation ( 27) shows that f (x l+1 , x l+2 . . . , x r ) is always positive for any non-zero vector X. The above proof shows that B is positive definite, and the convergence of PIA on Catmull-Clark subdivision solid is proved.

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k)

Examples and applications

In this section, we will compare the proposed interpolatory subdivision schemes with the butterfly volumetric subdivision approach [START_REF] Mcdonnell | Interpolatory, solid subdivision of unstructured hexahedral meshes[END_REF]], then we will present the applications of the proposed interpolatory subdivision methods for volumetric modeling with heterogeneous materials and isogeometric analysis of heat conduction and linear elasticity problems.

Comparison

The proposed interpolatory volumetric subdivision method can achieve a limit volumetric geometry with C 2 -continuity at regular part. McDonnell et al. has proposed a volumetric subdivision algorithm based on the well-known Lagrange polynomials, which is C 1 continuous at regular part. It can be seen as a natural generalization of the butterfly subdivision surface scheme to the volumetric case. Figure 4 presents three examples for the comparison of the butterfly method, the push-back method and the PIA approach. From these modeling examples, we can find that the PIA-based interpolatory volumetric subdivision results are much smoother than the butterfly volumetric subdivision results. From the cross-section of these models as shown in Figure 4(f)-4(h) and Figure 4(j)-4(l) , we can find that the PIA method can also achieve smoother interior mesh structure than the butterfly subdivision approach. Another comparison example is shown in Figure 5(b)-5(d).

Application in volumetric modeling with heterogeneous materials

Multi-material 3D printers are more and more popular for additive manufacturing, and the research on heterogeneous materials, especially the functionally graded material objects, is a hot topic in the field of intelligent manufacturing. The proposed interpolatory volumetric subdivision methods can be used not only for geometry interpolation, but also for material attribute interpolation in the field of volumetric material modeling with heterogeneous materials, which is a major benefit of our volumetric subdivision scheme.

Figure 6 shows a volumetric geometry and the material interpolation process by the interpolatory volumetric subdivision approach based on PIA, and the final volumetric material object can interpolate the material values at vertices of the input control lattice. Figure 5(e)-5(h), Figure 8(d)-8(i) and Figure 9(c)-9(e) illustrate several volumetric material modeling examples with heterogeneous material properties associated with geometry. In these examples, different material densities are assigned to different control points. By using our interpolatory volumetric approach, we can blend different densities smoothly and the resulting functionally graded material object can interpolate the material density at each control point exactly.

Application in isogeometric analysis

The proposed interpolatory volumetric subdivision methods have the following advantages for isogeometric analysis: (1) easy for imposing Dirichlet boundary conditions; (2) multi-resolution property; and (3) representation of complex geometry with arbitrary topology. Moreover, volumetric subdivision methods provide explicit regular shape functions for isogeometric application over complex geometry with arbitrary topology.In this paper, we also implement the isogeometric analysis framework with the proposed interpolatory volumetric subdivision for three-dimensional linear elasticity and heat conduction problems.

As an illustration problem, we consider the following three-dimensional heat conduction problem:

-∆T (x ) = f (x ) in Ω ⊂ R 3 T (x ) = 0 on ∂Ω ( 28
)
where ∆ is the Laplacian operator, Ω is the computational domain represented by the interpolatory Catmull-Clark subdivision volume, T (x ) is the unknown heat field, and f (x ) is the heat source function.

The IGA approximation of ( 28) is to find the numerical solution T h such that

Ω ∇T h • ∇vdxdy = Ω f • vdxdy, (29)
in which v is a test function. Each control point P i is associated with a basis function φ i . In isogeometric framework, the solution field T h has the following form

T h = n j=1 φ j T h j , (30)
in which T h j is the unknown control variable [START_REF] Pan | Isogeometric analysis based on extended Catmull-Clark subdivision[END_REF][START_REF] Pan | A unified method for hybrid subdivision surface design using geometric partial differential equations[END_REF]]. Then the numerical approximation problem (29) can be rewritten as

n j=1 T h j Ω ∇φ j • ∇φ i dxdy = Ω f • φ i dxdy, i = 1, • • • , n. (31)
Then we can obtain a linear system with unknown vector T, KT = b, in which K is the stiffness matrix and b is the load vector described as

K = Ω ∇φ j • ∇φ i dxdy and b = Ω f • φ i dxdy. ( 32
)
The stiffness matrix K and the load vector b over each element are evaluated by using three-dimensional Gauss-Legendre quadrature on the limit representation of the interpolatory Catmull-Clark subdivision volume. For the evaluation at the sampling points of the Gauss-Legendre quadrature, direct computation with cubic uniform B-spline basis functions is used for regular elements; for irregular elements, an approach based on diagonalization of the subdivision matrix is used to evaluate subdivision volumes at arbitrary parameter values [START_REF] Stam | Exact evaluation of Catmull-Clark subdivision surfaces at arbitrary parameter values[END_REF]; [START_REF] Burkhart | Iso-geometric finite element analysis based on Catmull-Clark subdivision solids[END_REF]]. For the convergence study of the IGA framework with interpolatory Catmull-Clark subdivision volumes, a heat conduction example of a sphere-like domain is presented in Figure 7. The initial level-0 hex-mesh in Figure 7(a) is subdivided three times to analyze the convergence, and the corresponding numerical error is computed with respect to a reference solution, which is solved with the level-4 hex-mesh as shown in Figure 7(d). After performing volumetric subdivision steps and IGA solving, we can find that the numerical error is concentrated in regions around irregular edges and vertices as shown in Figure 7(e)-7(h). Moreover, extraordinary edges and vertices also decrease the convergence rate as shown in Figure 7(i). The convergence rate is sub-optimal at 3.52 for the relative L 2 error.

Figure 5(i)-5(j), Figure 8(j)-8(k) and Figure 9(f)-9(h) illustrate several isogeometric simulation results for heat conduction and linear elasticity problems on complex geometry respectively.

Conclusion

With the help of limit point formula, two new methods are proposed to construct Catmull-Clark subdivision volumes that interpolate vertices of a given hex-mesh with arbitrary topology. The first method is based on the push-back operation, and the second method is based on progressive-iterative approximation with convergence proof. Some comparisons with examples are given, and applications in volumetric material modeling and isogeometric analysis are also presented to show the effectiveness of the proposed methods.

For the research on volumetric subdivision, there are still a few key issues to be addressed. Compared with subdivision surfaces, the extension of the subdivision scheme to higher dimensions is not trivial, such as the smoothness property and local-refinement theory, which will be our future research topics.

Figure 1 :

 1 Figure 1: Interpolatory Catmull-Clark volumetric subdivision example with the push-back method. (a) The original mesh; (b) new control lattices constructed by the push-back operation; (c) subdivision result after two levels of Catmull-Clark volumetric subdivision from the modified control lattice; (d) the interior view of the result.

Figure 2 :

 2 Figure 2: Examples of the push-back method with different parameters. (a) and (e): The original mesh; (b) and (f): interpolatory subdivision volume with λ ij = 0.05, µ F = 0.05 and γ C = 0.05; (c) and (g): interpolatory subdivision volume with λ ij = 0.25, µ F = 0.125 and γ C = 0.0625; (d) and (h): interpolatory subdivision volume with λ ij = 0.25, µ F = 0.25 and γ C = 0.25.

Figure 3 :

 3 Figure 3: Interpolatory volume example with the PIA method. (a): The original mesh; (b-d): new control lattices after one, three and ten PIA steps; (e-h): the corresponding iterative process of the limit volumes, in which the blue points are vertices of the original mesh.

Figure 5 :

 5 Figure 5: Comparison examples. (a-h): The input hexahedral mesh; the butterfly interpolatory subdivision method; the push-back method; and the PIA method. Blue, green and red vertices in (e) are assigned with material value 0, 50 and 100 respectively. (f-h): The corresponding material description over (b), (c) and (d) geometry, respectively. (i): Boundary condition for linear elasticity simulation on the interpolatory subdivision volume with the PIA method in (d); (j): IGA results for linear elasticity simulation with boundary condition in (i); (k): IGA results for heat conduction simulation.

Figure 6 :

 6 Figure 6: Volumetric geometry and material interpolation by the interpolatory volumetric subdivision approach based on PIA. (a): The original control lattice with different material values at vertices, where blue vertices have material value 0 and red vertices have material value 100; (b) -(d): volumetric subdivision results after one, three and ten PIA operations; (e): material distribution of subdivision volumes after two Catmull-Clark volumetric subdivision steps; (f) -(h): the corresponding material distribution after one, three and ten PIA operations with respect to (b) -(d).

Figure 7 :

 7 Figure 7: Convergence study on a sphere-like domain. (a): input control mesh; (b): interpolatory Catmull-Clark subdivision volume with the PIA method; (c):the interior view of (b); (d): reference solution to heat conduction problems; (e): error colormap of numerical solution on the level-0 subdivision volume; (f): error colormap of numerical solution on the level-1 subdivision volume; (g): error colormap of numerical solution on the level-2 subdivision volume; (h): error colormap of numerical solution on the level-3 subdivision volume; (i): convergence rate of relative L 2 error vs the number of degrees of freedom.

Figure 8 :Figure 9 :

 89 Figure 8: Bunny model. (a): The original mesh; (b): interpolatory volumetric subdivision result with the PIA method; (c): the interior view of (b); (d): the original control lattice with different material values, where blue, green and red correspond to 0, 50 and 100, respectively; (e) and (f): material distribution of the interpolatory subdivision solids by the PIA method; (g) the original control lattice with another set of material values at vertices; (h) material distribution of approximation Catmull-Clark subdivision volume with input (g); (i) material distribution of interpolatory Catmull-Clark subdivision volume with input (g); (j): boundary condition for linear elasticity simulation; (k): IGA results for linear elasticity simulation with boundary condition (j); (l): IGA results for heat conduction simulation.

Acknowledgment

The authors thank the anonymous reviewers for providing constructive suggestions and comments, that contributed to highly improve this work. The authors appreciate the joint support for this project by the National Natural Science Foundation of China (Grant Nos. 61761136010, 61772163), the Science Challenge Project (Grant No. TZ2016002), and the NSFC-Zhejiang Joint Fund for the Integration of Industrialization and Informatization (Grant No. U1909210).