domain: phys chim

text:
Cécile Huez

Séverine Renaudineau

Florence Volatron

Anna Proust
email: anna.proust@sorbonne-universite.fr

Dominique Vuillaume
email: dominique.vuillaume@iemn.fr

Experimental observation of the role of countercations on the electrical conductance of Preyssler-type polyoxometalate nanodevices

Keywords: polyoxometalate, counterion, electron transport, molecular electronics, conducIve-AFM

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

Introduc1on.

Ions inserted (inadvertently or by a controlled way) in solid-state molecular juncIons (MJs) play an important role to modulate or control the electron transport properIes. For example, the ions can modulate the conductance of the MJs as theoreIcally predicted, 1,2 and observed in several MJs, this conductance modulaIon depending on the chemical nature of the ions and the details of the molecule-ion conformaIon. [3][4][5][6][7] Not only the conductance of the MJs, but also the shape of the current-voltage (I-V) behavior can be in situ switched between a linear (symmetric) and a recIfying (asymmetric) behavior by inserIng/removing the ions or moving the precise posiIon (electrical field driven) of the ions in the MJs, [8][9][10][11] opening the way to molecular-scale programmable funcIonal devices and synapIc (neuro-inspired) behavior. 12,13 Polyoxometalates (POM) are candidates of choice owing to their remarkable reversible mulI-redox states 14,15 and photo-sImulable redox properIes, 16,17 they have ajracted a growing interest with several recent demonstraIon of proof-of-principles as siliconintegrated POM-based memory devices 18 and unconvenIonal (in-memory and neuromorphic) compuIng devices. [19][20][21] In soluIon, it is well known that the nature of the counterions has a strong impact on the redox potenIals for the POMs, hence on the electron transfer rates. 22 In the specific case of POM-based molecular juncIons, the interacIons with the unavoidable counterions are also controlling many properIes of the POM-based materials and devices. 23,24 For example, for a single POM on a metallic surface, the electron cloud density around the POM, as observed with a scanning tunneling microscope (STM), depends on the counterions that modify the hybridizaIon between the POM and the metal surface. 25 A recent theoreIcal study shows that the presence of the counterion (whatever its nature, here tetramethylamonium (TMA) vs. Cs + ) moves the LUMO (lowest unoccupied molecular orbital) of W18O54(SO3) 2-based MJs closer to the Fermi energy of the electrodes (∼0.08 eV above the Fermi energy in both the cases). 26 Albeit this similar posiIon of the LUMO, the computed current-voltage (I-V) curves showed slight differences of the shape (more asymmetric I-V with Cs + ). The authors suggest that this feature might be due to small differences in the amplitude of the calculated transmission coefficient of electrons in the MJ, clearly calling to more joint experimental and theoreIcal work. They also note that the presence of the counterions does not create addiIonal new conducIon channels in the MJs, but rather that they modify the potenIal landscape "viewed" by the POMs that transmit the electrons through the juncIon.

Here, we synthesized a series of Preyssler polyoxometalates (POMs), [Na ⊂P5W30O110] 14-, 27 P5W30 for short, stabilized with four different counterions, H + , K + , NH4 + and tetrabutylammonium (N(C4H9)4 + or TBA + ). We formed, by electrostaIc deposiIon, monolayers of these POMs on ultraflat Au surface prefuncIonalized with a posiIvely charged self-assembled monolayer (SAM) of amine-terminated alkylthiol chains following the protocol already developed and reported in our previous work for other POMs. 17,28 The electron transport (ET) properIes of these POMs were measured by conducIve-AFM (atomic force microscopy) and the current vs. voltage curves of these MJs were staIsIcally analyzed (hundreds of I-V traces) to determine the energy posiIon of the molecular orbitals involved in the ET and the electronic coupling between the molecules and the electrodes. We demonstrate that the energy posiIon of the LUMO of the P5W30 with respect of the Fermi energy of the electrodes increases from ∼ 0.4 eV to ∼ 0.7 eV and that that electrode coupling energy evolves from ∼ 0.05 to 1 meV depending on the nature of the counterion in the order K + , NH4 + , H + and TBA + .

Synthesis and Structural Characteriza1ons

The synthesis of the P5W30 POMs with various counterions is schemaIzed in the figure 1a. The potassium salt of [Na⊂P5W30O110] 14-can be obtained under hydrothermal condiIons. 29 The sample we used corresponds to the molecular

f o r m u l a K 1 3 . 3 N a 0 . 7 [ N a ⊂ P 5 W 3 0 O 1 1 0 ] . 2 7 H 2 O .
T h e a m m o n i u m s a l t (NH4)14[Na⊂P5W30O110].25H2O has been prepared according to reference 30. The acid salt H14[Na⊂P5W30O110].44H2O was prepared by ion exchange of a soluIon of the potassium salt on a Dowex 50W-X8 proton exchange resin, following a published procedure. 31 Acido-basic ItraIon of heteropolyacids with caIonhydroxide is a common route to introduce other caIons and it has also been applied to this Preyssler-type phosphotungsIc acid. 32 In our hands, the acidobasic ItraIon of H14[Na⊂P5W30O110] with tetrabutylammonium hydroxide was not complete, despite careful colorimetric monitoring of the equivalence. We ended with the mixed salt TBA10H4[Na⊂P5W30O110].24H2O.

We end-up with 4 POMs: K + 13.3Na + 0.7[Na⊂P5W30O110] 14-(KP5W30 for short); H + 14[Na⊂P5W30O110] 14-(HP5W30 for short); TBA + 10 H + 4 [Na⊂P5W30O110] 14-(TBAP5W30 for short) and (NH4) + 14[Na⊂P5W30O110] 14-(NH4P5W30 for short). See more details, IR and 31 P NMR characterizaIons, TGA and elemental analysis in the SupporIng InformaIon (secIon 1).

Figure 1. (a) Schema'c of the synthesis routes (ca'ons and POM not at scale). (b)

Schema'c descrip'on of the two-step fabrica'on of the self-assembled monolayers.

These POMs (dissolved in H2O, or acetonitrile/isopropanol for HP5W30 and TBAP5W30) were electrostaIcally deposited on ultra-flat template-sIpped TS Au substrates funcIonalized with a 6-aminohexane-1-thiol hydrochloride (HS-(CH2)6-NH3 + /Cl -) self-assembled monolayer (SAM), or C6 SAM for short, following the protocol already developed and reported in our previous work for other POMs (Fig. 1b and also see secIon 2 in the SupporIng InformaIon). 17,28 The thicknesses (systemaIcally measured by ellipsometry, see secIon 3 in the SupporIng InformaIon) are summarized in Table 1.

Table 1.

Thickness measured by ellipsometry at each step of the fabrica'on of the molecular junc'on.

We note that the P5W30 POM has an oval shape with a long axis of 1.8 nm and a short axis of 1.3 nm. Thus the 1.4 nm measured thickness of the TBAP5W30 sample indicates a denser monolayer with possibly a small percentage of the TBAP5W30 molecule standing upright on the surface (verIcal POM orientaIon), while the thinner values for the 3 other samples can be ascribed to a less compact monolayer whatever the POM orientaIon (disordered monolayers with various POM orientaIons). The thicker TBAP5W30 SAM may also be caused by the bulky TBA ions (diameter ≈ 0.5-1 nm depending on the conformaIon of the butyl chains) mainly intercalated between the POMs and the electrode, while small caIons can be distributed around the POM.

Electron Transport Proper1es

Figure 2 shows the 2D-histograms (or "heat map") of the I-V curves (between 400 and 500) acquired by C-AFM at different locaIons on the POM monolayer for the four TS Au-C6/POM//Pt MJs ("-" denotes a chemical bond, "/" an electrostaIc contact and "//" a mechanical contact). The main feature is that these I-V distribuIons have almost the same level of currents at large voltages, i.e. |V|>1V, while the currents in the voltage range [-0.6 V ; 0.6 V] gradually increase with changing the counterion in the order from K + , NH4 + , H + to TBA + . This is clearly seen in Fig. 3a, plowng the mean current -voltage (Ī-V) for the 4 samples and their counterion dependence at several voltages (Fig. 3b). The mean current Ī is around 10 -8 A at +/-1.5V, while it varies over 2 decades at +/-0.5 V. This observaIon is also supported by plowng and comparing the current histograms The Ī-V shapes for all samples show a bump at large voltages, i.e. |V|>0.7 -1V, while the low-voltage Ī-Vs are well fijed between -0.6 and 0.6 V with a single-energy level (SEL) model given by the following analyIcal expression: [START_REF] Daja | Electronic transport in mesoscopic systems[END_REF][START_REF] Cuevas | Molecular Electronics: An introduc'on to theory and experiment[END_REF] (1)

with ε0 the energy of the molecular orbital (MO), here LUMO involved in the transport (with respect to the Fermi energy of the electrodes), Γ1 and Γ2 the electronic coupling energy between the MO and the electron clouds in the two electrodes, e the elementary electron charge, h the Planck constant and N the number of molecules contribuIng to the ET in the molecular juncIon (see details in secIon 6 the SupporIng InformaIon). Among several limitaIons (secIon 6 in the SupporIng InformaIon), the SEL model is a low temperature approximaIon and the temperature broadening of the Fermi funcIon in the electrodes is not taken into account. However, it was shown that it can be reasonably used at room temperature for voltages below the sharpened increase of the current (not observed here in the -0.6/0.6 V window) characterizing the transiIon between the off-resonance and resonant transport condiIons at which the broadening of the Fermi funcIon modify the I-V shape. [START_REF] Grüter | Resonant tunnelling through a C 60molecular juncIon in a liquid environment[END_REF][START_REF] Bourgoin | Molecular Electronics[END_REF][START_REF] Brunner | Gaining Microscopic Insight into Molecular JuncIons by Transport Experiments[END_REF] Fig 3c

shows typical examples of the SEL fits on the mean Ī-V and Fig. 3d shows the evoluIon of the fijed parameters versus the nature of the counterions, clearly revealing an increase of the molecular energy level ε0 from 0.42 eV (for KP5W30) to 0.69 eV (for TBAP5W30). The electrode coupling energies, Γ1 and Γ2, follow a same trend, increasing from ∼0.08 to 1-2 meV in the same counterion order (Fig. 3d and Table 2). The same trend is confirmed from a staIsIcal analysis by fiwng with the SEL model all the individuals I-Vs of the datasets shown in figure 2. The staIsIcal distribuIons of ε0, Γ1 and Γ2 for the 4

MJs are shown in Figs. S12 and S13 (SupporIng InformaIon) and the mean values of these distribuIons are summarized in Fig. 3d and Table 2.

The same results are obtained by analyzing the I-V curves with the transiIon voltage spectroscopy (TVS) method. [START_REF] Beebe | TransiIon from Direct Tunneling to Field Emission in Metal-Molecule-Metal JuncIons[END_REF][START_REF] Beebe | Measuring RelaIve Barrier Heights in Molecular Electronic JuncIons with TransiIon Voltage Spectroscopy[END_REF][START_REF] Huisman | InterpretaIon of TransiIon Voltage Spectroscopy[END_REF][START_REF] Mirjani | Advantages and limitaIons of transiIon voltage spectroscopy: A theoreIcal analysis[END_REF][START_REF] Bâldea | Ambipolar transiIon voltage spectroscopy: AnalyIcal results and experimental agreement[END_REF] The energy level and the electrode coupling energy determined by the fit of the SEL model and by the TVS method are in very good agreement (secIon 7, Figs. S14 and S15 in the SupporIng InformaIon).

Albeit the staIsIcal analysis shows a large distribuIon of values (e.g. 2 decades for the parameters Γ1 and Γ2 -Fig. S13 in the SupporIng InformaIon, and Table 2), the same trends and order of magnitude of these parameters are observed from the analysis of the mean Ī-V, the staIsIcal analysis of the full dataset, and using both the SEL and TVS methods, making the experimental observaIons conclusive.

Table 2. SEL parameters ε0, Γ1 and Γ2 fined on the mean Ī-V (Fig. 3a) and mean values ⟨ε0⟩, ⟨Γ1⟩ and ⟨Γ2⟩ deduced from the sta's'cal analysis (see text). For the energy level ⟨ε0⟩, the error is the standard devia'on of the normal distribu'on (Fig. S12). For the ⟨Γ1⟩ and ⟨Γ2⟩ parameters, the values in parenthesis are the min and max values at FWHM of the log-normal distribu'ons (Fig. S13).

Discussion

The low-bias conductance of the P5W30-based MJs increases by a factor ∼ 100 between the K + and the TBA + caIons (Figs. 3b-c). However, the weak increase of the energy posiIon of the LUMO (factor ∼ 1.6, Fig. 3d, Table 2) should have induced a decrease of the current due to a higher electron energy barrier (ε0-εF) at the molecule/electrode interface. A factor ∼ 4.5 is esImated from Eq. ( 1)

assuming the same molecule/electrode coupling energy for simplicity (Γ1=Γ2=0. meV) in the simulated I-V curves (Fig. S16 in the SupporIng InformaIon). The increase of the measured current is clearly related to the larger increase (factor ∼ 20, Fig. 3d and Table 2) of the molecule/electrode coupling energy. We examine several hypotheses to explain these results. The trend observed in soluIon for the size-depending behavior of the alkali salts of POMs is barely transposable: indeed, a posiIve shi} of the reducIon potenIal (which would correspond here to a stabilizaIon of the LUMO and a decrease of ε0) and a greater electron transfer rate to POMs associated with larger alkali metal caIons had been ascribed to an increase of the ion pairing energy related to a decrease of the caIon solvaIon and the formaIon of more inImate caIon/anion pairs. 22 Although we cannot completely exclude the presence of some remaining solvent molecules in our MJs, solvaIon is not really relevant here. In a recent theoreIcal report, the efficiency of caIons to compensate for the charging of a hexavanadate POM has been addressed by simultaneous addiIon of extra caIons and electrons: in the gas phase, Li + and H + were found the most efficient, followed by NH4 + and K + . [START_REF] Sorokina | What is the maximum charge uptake of Lindqvist-type polyoxovanadates in organic-inorganic heterostructures?[END_REF] Although the origin of this effect was not further discussed it could similarly impact the electron transport across the MJs. Last but not least, the tetrabutylammonium caIon is very different from the other caIons used in this study, and because of its volume and lower charge density the so}er caIon of the series. Among many factors, the energy level alignment in the MJs (LUMO of the P5W30 POM with respect to the electrode Fermi energy) may depend on the interface dipole. The larger TBA counterion could probably induce the largest local dipole (larger distance between the posiIve and negaIve charges), and if we assume that the counterions are mostly inserted between the POM monolayer and the electrode (due to steric hindrance, vide supra the thickness measurements), then a global dipole could exist at the interface with an orientaIon leading to an upwards larger LUMO offset from the Fermi level due to the dipole-induced vacuum level shi} (Fig. S17 in the SupporIng InformaIon).

On the contrary, smaller caIons create weaker dipoles randomly arranged around the POM, [START_REF] Fernández | Polyoxometalates with Internal CaviIes: Redox AcIvity, Basicity, and CaIon EncapsulaIon in [Xn+P5W30O110](15n)-Preyssler Complexes, with X = Na+, Ca2+, Y3+, La3+, Ce3+, and Th4+[END_REF] resulIng in a weaker (or even negligible) average dipole at the interface, thus a weaker LUMO shi}. Note that we cannot exclude, at this stage without further theoreIcal exploraIons, that the atomisIc details of the POM/electrode contact configuraIon play also a role (as also known in many other MJs). In the case of the W18O54(SO3) 2-based MJs, 26 a change by a factor ca.

4 in the amplitude of the current and a variaIon of the LUMO between ca. 0.4 and 1.2 eV was theoreIcally predicted whether the POM is oriented with its long axis perpendicular or horizontal between the electrodes and how (atomic configuraIon) the POM is connected to the Au electrode. From the thickness measurements (vide supra), we also suggest possible different orientaIon of the POMs, thus further detailed calculaIons (out of the scope of this work) are called to a reach a conclusion.

The evoluIon of the molecule/electrode coupling energy is most intriguing. Nevertheless, it has been recently demonstrated that the counterion can mediate the intermolecular electron transfer between two adjacent and nearby POMs, [START_REF] Werner | Hybrid Molecular Magnets with Lanthanide-and CountercaIon-Mediated Interfacial Electron Transfer between Phthalocyanine and Polyoxovanadate[END_REF] leading to a possible POM-caIon-POM electron conducIon channel. In the situaIon schemaIcally depicted above for the case of the TBA + counterion (Fig. S17 in the SupporIng InformaIon), we extrapolate that the same kind of mechanism could enhance the POM-caIon-electrode electronic coupling, compared to a more random organizaIon around the POMs for the other counterions (Fig. S17). Moreover, the fact that the POM/electrode interacIon also depends on the nature of the counterions has been observed by STM images of a single POM deposited on a metal surface. In this experiment, the size of the electron cloud density measured by STM largely extends the geometric size of the POM due to the hybridizaIon of the POM and the counterion with the metal surface, with a stronger hybridizaIon with a TBA + caIon rather than the POM encapsulated in a cyclodextrin cage (caIon free). 25 This hybridizaIon is also the key factor fixing the molecule/electrode energy coupling in the SEL model (Eq. 1). Thus, we hypothesis that the TBA + counterion could lead to a stronger POM/electrode hybridizaIon than the other ones.

Finally, we note that our energy values with the TBA + counterion (ε0 = 0.66-0.69 eV, Γ1 and Γ2 ≈ 1-2 meV) are of the same order as values determined by STM experiments on a single similar Preyssler POM ([DyP5W30O110] 12-) in the case of a strong molecule/electrode coupling (established by moving the STM Ip in close contact to the molecule in that case): ε0 ≈ 0.7 eV, Γ1 and Γ2 ≈ 1-10 meV. [START_REF] Sherif | Current recIficaIon in a single molecule diode: the role of electrode coupling[END_REF] At higher voltages (|V| > 0.7 V), the bump of the current cannot be fijed with the SEL model. We hypothesis that this "addiIonal" current can be due to a second level entering the energy window at higher voltages (Fig. S11 in the SupporIng InformaIon). Such a behavior was theoreIcally simulated in the case of a two-channel electron transmission in the MJs. [START_REF] Migliore | Nonlinear Charge Transport in Redox Molecular JuncIons: A Marcus PerspecIve[END_REF] This second level is tentaIvely ajributed to the LUMO+1 of the POMs. Our results suggest that the LUMO+1 is less sensiIve to the counterions, since the currents at |V| > 0.7 V are less dependent on them (Fig. 3b and Fig. S10 in the Supplementary InformaIon).

We also note that the POM samples with the NH4 + and H + countercaIons display a slightly higher current at +1V than at -1V (Fig. S10 and Table S1 in the SupporIng InformaIon), while no voltage dependence asymmetry is observed at lower voltages (Fig. S9). The mean asymmetry raIo (also known as recIficaIon raIo) R=Ī(+1V)/Ī(-1V) is weak (≈1.7 and ≈3.2 for HP5W30 and NH4P5W30, respecIvely) with a difference in log-Ī (log-mean current) of 0.24 for HP5W30 and 0.5 for NH4P5W30, smaller than log-σ (log standard deviaIon, in the range 0.84-1.1). It has been pointed out that weak R values must be considered with cauIon without a solid staIsIcal analysis. [START_REF] Nijhuis | Molecular RecIficaIon in Metal-SAM-Metal Oxide-Metal JuncIons[END_REF] A two-sample t-test on the current datasets at +1 and -1V was used to assess if these datasets staIsIcally and significantly differ from each other (see secIon 9 in the SupporIng InformaIon).

The t-test results show that the null hypothesis (same mean value for the two datasets) can be rejected and this weak current asymmetry is staIsIcally significant. The origin of this effect may be due to different contact geometries of the POM to the electrodes (at the atomic level) in the presence of the different conterions, with an asymmetric hybridizaIon between the POM and the two electrodes for the LUMO+1 in the case of H + and NH4 + counterions (no asymmetry is observed at lower bias, -0.6/0.6V, involving only the LUMO in the electron transport mechanism (Fig. 3c and Fig. S9 in the SupporIng InformaIon). A similar effect was theoreIcally predicted for W 18 O 54 (SO 3 ) 2- POM with TMA (symmetric I-V) or Cs + (asymmetric I-V) counterions. 26 Further theoreIcal work (out of the scope of this work) would be necessary to clarify this point.

Conclusion

In conclusion, we have demonstrated that the electron transport properIes, at the nanoscale, of P5W30-based molecular juncIons depend on the nature of the counterions. From a staIsIcal study (hundreds of I-V traces) using a simple analyIcal model derived of the Landauer-Imry-Büwker formalism for charge transport in nanoscale devices, we have found that the energy posiIon of the LUMO of the P5W30 with respect of the Fermi energy of the electrodes increases from ∼ 0.4 eV to ∼ 0.7 eV and that that electrode coupling energy evolves from ∼ 0.05 to 1 meV depending on the nature of the counterion in the order from K + , NH4 + , H + to TBA + . We have suggested that these variaIons could be due to a counterion-dependent POM/electrode interface dipole and a molecule/electrode hybridizaIon mediated by the counterions. These non-trivial POM-counterionmetal electrode interacIons clearly require further theoreIcal studies to elucidate the influence of the counterions on the electron transport properIes of POM-based molecular devices.

Methods

Sample fabrica5on.

Bonom metal electrode fabrica'on. Template stripped gold ( TS Au) substrates were prepared according to the method previously reported. [START_REF] Hegner | Ultralarge atomically flat template-stripped Au surfaces for scanning probe microscopy[END_REF][START_REF] Weiss | Influence of defects on the electrical characterisIcs of mercury-drop juncIons: self-assembled monolayers of n-alkanethiolates on rough and smooth silver[END_REF][START_REF] Weiss | Si/SiO 2-Templated FormaIon of Ultraflat Metal Surfaces on Glass, Polymer, and Solder Supports: Their Use as Substrates for Self-Assembled Monolayers[END_REF] (secIon 2 in the SupporIng InformaIon).

Self-assembled monolayers. The SAMs on TS Au were fabricated following a protocol developed and opImized in our previous works for the electrostaIc immobilizaIon of POMs on amine-terminated SAMs. 17,28 (secIon 2 in the SupporIng InformaIon).

Spectroscopic ellipsometry.

The thickness of the SAMs was measured by spectroscopic ellipsometry (UVISEL ellipsometer (HORIBA), secIon 3 in the SupporIng InformaIon) CAFM in ambient condi5ons. We measured the electron transport properIes at the nanoscale by CAFM (ICON, Bruker) at room temperature using a Ip probe in plaInum/iridium (see secIon 4 in the SupporIng InformaIon).

Associated content

The SupporIng InformaIon is available free of charge at xxxxxx.

SUPPORTING INFORMATION 1. Synthesis Details

All the chemicals were used as supplied. K12,5Na1,5[Na⊂P5W30O110], 1 (NH4)14[Na⊂P5W30O110], 2 and H14[Na⊂P5W30O110], 3 were prepared according to published procedures and their purity checked by IR and 31 P NMR spectroscopy.

The synthesis of TBA10H4[Na⊂P5W30O110] was adapted from the literature. 4 IR spectra (KBr pellets) have been recorded in transmission from 250 to 4000 cm -1 on a Jasco FT/IR 4100 spectrometer (Resolu\on 4 cm -1 ). 31

Synthesis and characteriza/on of K12,5Na1,5[Na⊂P5W30O110]

To 33g of Na2WO4.2H2O dissolved in 30mL of water were added 26.5mL of H3PO4 (85%). This mixture was placed in a solvothermal synthesis bomb (Parr Model 4748), heated to 120°C overnight. Once the yellow solu\on has been cooled to room temperature, 15mL of H2O was added followed by 10g of KCl. The precipitate was membrane filtered and washed with 2M potassium acetate and methanol. Recrystalliza\on in hot water was necessary to purify the product (with a pause \me < 2h to avoid the forma\on of crystals of [a-P2W18O62] 6-). About 8g of crystals are then collected (yield about 30%). A second recrystalliza\on can be performed if the product s\ll contains impuri\es. For each step, the purity was assessed by 31 P NMR: the small signal at -12,55 ppm corresponds to [a-P2W18O62] 6-. 5 Thermogravimetric analysis gave 27 water molecules of crystalliza\on. IR (KBr, cm -1 ): n 3538 (m), 3432 (m), 1616 (m), 1164 (m), 1082 (m), 1017 (w), 982 (sh), 936 (s), 911 (s), 782 (s), 670 (sh), 571 (w), 538 (w), 464 (w), 349 (m), 311 (w), 286 (w). 31 P NMR (D2O): δ = -9.61 ppm. Elemental analysis: found (calculated for K13,3Na0,7[Na⊂P5W30O110].27H2O MW = 8475.07 g mol -1 ) K: 6.12 (6.14), Na: 0.42 (0.46), P: 1.76 (1.83), W: 65.23 (65.08) %.

Synthesis of (NH4)14[Na⊂P5W30O110]

The synthesis was performed following the published procedure and thermogravimetric analysis carried out on an air-dried sample gave 25 water molecules of crystalliza\on. IR (KBr, cm -1 ): n 3551 (m), 3443 (m), 3127 (m), 1615 (m), 1400 (s), 1165 (m), 1080 (m), 1017 (m), 982 (w), 934 (s), 910 (s), 777 (s), 571 (w), 539 (w), 470 (w), 359 (m), 311 (w). 31 P NMR (D2O): δ = -9.61 ppm.

Figure S3. IR spectrum of (NH4)14[Na⊂P5W30O110] (KBr pellet).

Synthesis of H14[Na⊂P5W30O110]

First, we prepared the Dowex column 50W-X8 in its acidic form. In a beaker, the very fine par\cles of the resin were separated from the mass by successive sedimenta\ons and seqlings in dis\lled water. The resin was then washed with dis\lled water un\l the orange colora\on resul\ng from the degrada\on of the polymer disappeared. Finally, the resin was placed in a column (diameter 1.5 cm, length 40 cm), washed three \mes with its dead volume of 1M HCl (∼3x20 mL)

and rinsed with water un\l the pH obtained for the eluent reached that of dis\lled water. 6 In a second step, 8.2g of K13,3Na0,7[Na⊂P5W30O110].27H2O were dissolved in 160 mL of H2O. The solu\on was transferred to the column and passed through the resin at a rate of one drop per second. Water (∼40mL) was then passed through the column, as long as the pH of the eluent remained acidic to ensure that all the product was collected. The eluent was then evaporated at a temperature of 60-65°C to recover the product as a yellow and s\cky solid, which was crushed twice with ethanol, then dried with diethylether (yield 7.8g).

According to thermogravimetric analysis, 44

Self-Assembled Monolayer

Ultraflat template-stripped gold surfaces ( TS Au), with rms roughness of ∼0.4 nm were prepared according to the method already reported. [7][8][9] In brief, a 300-500 nm thick Au film was evaporated on a very flat silicon wafer covered by its na\ve SiO2 (rms roughness of ∼0.4 nm), which was previously carefully cleaned by piranha solu\on (30 min in 7:3 H2SO4/H2O2 (v/v); Cau/on: Piranha solu\on is a strong oxidizer and reacts exothermically with organics), rinsed with deionized (DI) water, and dried under a stream of nitrogen. Clean 10x10 mm pieces of glass slide (ultrasonicated in acetone for 5 min, ultrasonicated in 2-propanol for 5 min, and UV irradiated in ozone for 10 min) were glued on the evaporated Au film (UVpolymerizable glue, NOA61 from Epotecny), then mechanically peeled off providing the TS Au film aqached on the glass side (Au film is cut with a razor blade around the glass piece).

The self-assembled monolayers (SAMs) of 6-aminohexane-1-thiol hydrochloride (HS-(CH2)6-NH3 + / Cl -) were prepared following a protocol op\mized and described in a previous work for the electrosta\c immobiliza\on of POMs on amineterminated SAMs. 10,11 The freshly prepared TS Au substrates were dipped in a solu\on of 6-aminohexane-1-thiol hydrochloride (Sigma-Aldrich) at a concentra\on of 10 -3 M in ethanol overnight in the dark. The samples were rinsed in ethanol for 5 min and then ultrasonically cleaned 5 min in deionized (DI)

water. These SAMs were treated by a PBS (phosphate-buffered saline, pH=7.4) solu\on for 2 hours, followed by ultra-sonica\on in DI water for 5 minutes. The substrates were finally washed with ethanol and dried under nitrogen flow. It was found that the PBS treatment removes the forma\on of aggregates on the aminoalkylthiol SAMs as well as avoids clustering of POMs during the electrosta\c deposi\on, likely because this treatment op\mizes the ra\o of NH3 + /NH2 on the surface. 10 The electrosta\c deposi\on of the POMs was done by immersion of these SAMs in a solu\on of POM at a concentra\on of 10 -3 M in H2O, except for HP5W30 and TBAP5W30 dissolved in ACN/isopropanol, for one to few hours. We checked by ellipsometry that the thickness of the POM layer was independent of the immersion \me when the immersion \me was longer than 1h.

Ellipsometry Measurements

We recorded spectroscopic ellipsometry data (on ca. 1 cm 2 samples) in the visible range using a UVISEL (Horiba Jobin Yvon) spectroscopic ellipsometer equipped with DeltaPsi 2 data analysis so]ware. The system acquired a spectrum ranging from 2 to 4.5 eV (corresponding to 300-750 nm) with intervals of 0.1 eV (or 15 nm). The data were taken at an angle of incidence of 70°, and the compensator was set at 45°. We fit the data by a regression analysis to a film-on-substrate model as described by their thickness and their complex refrac\ve indexes. First, a background for the substrate before monolayer deposi\on was recorded. We acquired three reference spectra at three different places of the surface spaced of few mm. Secondly, a]er the monolayer deposi\on, we acquired once again three spectra at three different places of the surface and we used a 2-layer model (substrate/SAM) to fit the measured data and to determine the SAM thickness.

We employed the previously measured op\cal proper\es of the substrate (background), and we fixed the refrac\ve index of the monolayer at 1.50. 12 We note that a change from 1.50 to 1.55 would result in less than a 1 Å error for a thickness less than 30 Å. The three spectra measured on the sample were fiqed separately using each of the three reference spectra, giving nine values for the SAM thickness. We calculated the mean value from this nine thickness values and the thickness incer\tude corresponding to the standard devia\on. Overall, we es\mated the accuracy of the SAM thickness measurements at ± 2 Å. 13

C-AFM Measurements

We measured the electron transport proper\es at the nanoscale by CAFM (ICON, Bruker) at room temperature using a \p probe in pla\num/iridium. We used a "blind" mode to measure the current-voltage (I-V) curves and the current histograms: a square grid of 10×10 was defined with a pitch of 50 to 100 nm. At each point, the I-V curve is acquired leading to the measurements of 100 traces per grid. This process was repeated several \mes at different places (randomly chosen) on the sample, and up to several thousand of I-V traces were used to construct the current-voltage histograms (shown in Fig. 3, main text).

The \p load force was set at ≈6-9 nN for all the I-V measurements, a lower value leading to too many contact instabili\es during the I-V measurements. Albeit

Details of the I-V measurements

From the datasets shown in Fig. 2, we extracted the distribu\on of the current measured at +/-0.5V (Fig. S9) and +/-1V (Fig. S10) to support the results shown in Fig. 3b (evolu\on of the current with the nature of the counterions). S1. S1.

Table S1. Values of log-mean (log-Ī), mean current Ī, and a log standard devia'on (log-𝜎) of the current dispersion (log-normal distribu'on) shown in Fig. S10.

Single Energy-Level Model

The single-energy level (SEL) model (Eq. 1 main text), considers that: i) a single molecular orbital (MO) dominates the charge transport, ii) the voltage mainly drops at the molecule/electrode interface and iii) that the MO broadening is described by a Lorentzian or Breit-Wigner distribu\on. 16,17 The simple energy scheme (Fig. S11) is described by ε0 the energy of the MO involved in the transport (with respect to the Fermi energy of the electrodes), Γ1 and Γ2 the electronic coupling energy between the MO and the electron clouds in the two electrodes, e the elementary electron charge, h the Planck constant and N the number of molecules contribu\ng to the ET in the molecular junc\on (assuming independent molecules conduc\ng in parallel, i.e. no intermolecular interac\on [18][19][20] the \p/SAM interface in some cases when the Young modulus of the SAM is reasonably known, 14,15,[21][22][23] this is not the case here for the POM/alkyl SAM system for which the Young modulus has not been determined. Consequently, we use N=1 throughout this work. This means that the Γ1 and Γ2 values are "effec\ve" coupling energies of the SAM with the electrodes and they are used only for a rela\ve comparison of the POM SAMs measured with the same C-AFM condi\ons in the present work and they cannot be used for a direct comparison with other reported data (as for example from single molecule experiments). We also note that the exact value of N has no significant influence on the fiqed parameter ε0. The fits were done with the rou\ne included in ORIGIN so]ware, 24 using the method of least squares and the Levenberg Marquardt itera\on algorithm.

Transi1on Voltage Spectroscopy

We also used the transi\on voltage spectroscopy (TVS) [25][26][27][28][29] to analyze the I-V curves. Plo•ng |V 2 /I| vs. V (Fig. S14), 30 we determine the transi\on voltages VT+ and VT-for both voltage polari\es at which the bell-shaped curve is maximum.

This threshold voltage indicates the transi\on between off resonant (below VT) and resonant (above VT) transport regime in the molecular junc\ons and can therefore be used to es\mate the loca\on of the energy level. In Fig. S14, the thresholds VT+ and VT-are indicated by the ver\cal arrows (with values) and determined from the max of a 2nd order polynomial func\on fiqed around the max of the bell-shaped curves (to cope with noisy curves). The value of ε0-TVS is es\mated by: 27

(2)

We also determined an average value of the electrode coupling energy ΓTVS using this rela\onship: 31,32 (3)

with G(0) the zero-bias conductance (Fig. S15), G0 the conductance quantum

Two-sample t-test

To determine if the small current asymmetry at +1V and -1V for the HP5W30 and NH4P5W30 is sta\s\cally significant, we performed a two-sample t-test using the rou\ne provided by ORIGIN so]ware. 24 Since the t-test applies for normally distributed popula\ons and the currents are log-normal distributed (Fig. S10), we consider the decimal logarithm of the currents taken from the datasets shown in Fig. 2 (main text) at -1V and +1V as the two popula\ons 1 and 2 (P1, P2). We test the null hypothesis that the mean values are equal, mean(P1)-mean(P2)=0. At a significance level of 0.05, the data yielded a probability p-value that the null hypothesis is true of 4x10 -9 for HP5W30 and 6x10 -12 for NH4P5W30, which are less than our 0.05 significance level, and thus we can reject the null hypothesis. In other words, the weak current asymmetry can be considered as sta\s\cally significant.

  (here at +/-0.5V and +/-1V), Figs. S9 and S10 in the SupporIng InformaIon.

Figure 2 .

 2 Figure 2. 2D histograms (400-500 I-V traces) of the four samples in a semi-log |I|-V plot. The solid black lines are the calculated mean Ī-V curves.

Figure 3 .

 3 Figure 3. (a) Mean Ī-V curves for the four samples (from the datasets shown inFig. 3), (b) counterion dependence of the mean current Ī at +/-1.5 V and +/-0.5 V, (c) fits (solid lines) of the SEL model in the voltage window [-0.6 V ; 0.6 V] (fit R 2 = 0.995 for K + , 0.998 for NH4 + , 0.995 for H + , 0.999 for TBA + ) , (d) counterion evolu'on of the SEL model parameters ε0, Γ1 and Γ2.

  Synthesis details, IR and NMR characterizaIons; fabricaIon of self-assembled monolayers; ellipsometry measurements; protocol of C-AFM measurements; supplemental details on the I-V measurements (staIsIcal data); details on the single energy level (SEL) model and on the staIsIcal distribuIon of the model paramaters; transiIon voltage spectroscopy (details on the technique and results); supplementary figures for the discussion. Author Contribu5ons C.H., S.R., F.V. and A.P. synthesized and characterized the POMs and fabricated the SAMs. C.H. did all the C-AFM measurements, C.H. and D.V. analyzed the data. A.P. and D.V. conceived and supervised the project. The manuscript was wrijen by D.V. with the contribuIons of all the authors. All authors have given approval of the final version of the manuscript. Experimental observa1on of the role of counterca1ons on the electrical conductance of Preyssler-type polyoxometalate nanodevices. Cécile Huez, 1, ‡ Séverine Renaudineau, 2 Florence Volatron, 2 Anna Proust 2 and Dominique Vuillaume. 1 1) Ins'tute for Electronics Microelectronics and Nanotechnology (IEMN), CNRS, University of Lille, Av. Poincaré, Villeneuve d'Ascq, France 2) Ins'tut Parisien de Chimie Moléculaire (IPCM), CNRS, Sorbonne Université, 4 Place Jussieu, Paris, France ‡) Now at : Ins'tuto de Ciencia Molecular (ICMoL), University of Valencia, Spain.

  P NMR spectra (121 MHz) have been recorded on a Bruker Avance II 300MHz spectrometer and the chemical shi]s referenced to external 85 % H3PO4. Thermogravimetric analyses (TGA) have been performed on SDT-Q600 TA Instrument thermobalance under an air flow to assess the amount of water molecules. Elemental analyses have been carried out at Crealins and at the Ins\tut des Sciences Analy\ques, CNRS, Villeurbanne, France.

Figure S4. 31 P

 31 Figure S4.31 P NMR spectrum of (NH4)14[Na⊂P5W30O110] in D2O.

Figure S5 .

 S5 Figure S5. IR spectrum of H14[Na⊂P5W30O110] (KBr pellet).

Figure S6. 31 P

 31 Figure S6. 31 P NMR spectrum of H14[Na⊂P5W30O110] in D2O.

Figure S7 .

 S7 Figure S7. IR spectrum of TBA10H4[Na⊂P5W30O110] (KBr pellet).

Figure S8. 31 P

 31 Figure S8. 31 P NMR spectrum of TBA10H4[Na⊂P5W30O110] in CD3CN.

  larger than the usual load force (2-5 nN) used for CAFM on SAMs, this value is below the limit of about 60-70 nN at which the SAMs start to suffer from severe degrada\ons. For example, a detailed study (Ref.14) showed a limited straininduced deforma\on of the monolayer (≲ 0.3 nm) at this used load force. The same conclusion was confirmed by our own study comparing mechanical and electrical proper\es of alkylthiol SAMs on flat Au surfaces and \ny Au nanodots.15

Figure S9 .

 S9 Figure S9. Histograms of log(current) measured at 0.5V and -0.5V for the four samples. The solid lines are the fit by a log-normal distribu'on with the values of the log-mean (log-Ī) and a log standard devia'on (log-𝜎) marked on the panels

Figure S10 .

 S10 Figure S10. Histograms of log(current) measured at 1V and -1V for the four samples. The solid lines are the fit by a log-normal distribu'on with the values of the log-mean (log-Ī) and a log standard devia'on (log-𝜎) marked on the panels

Figure S11 .

 S11 Figure S11. Schema'c energy diagram of the molecular junc'on at 0V, at a moderate posi've voltage applied on the Au substrate V<0.6-0.7 V and at higher voltages V>0.7V.

Figure S12 .

 S12 Figure S12. Histograms of the energy values ε0 deduced from the fit of the SEL model on the complete I-V dataset shown in Fig. 2 (main text). The solid lines are the fits by a normal distribu'on, the mean values are shown in the panels.

Figure S13 .

 S13 Figure S13. Histograms of the coupling energy values Γ1 and Γ2 deduced from the fit of the SEL model on the complete I-V dataset shown in Fig. 2 (main text). The solid lines are the fits by a log-normal distribu'on, the mean values are shown in the panels.

(2e 2 /

 2 h=7.75x10-5 S, e the electron charge, h the Planck constant) and N the number of molecules in the junc\on. G is calculated from the slope of the I-V curve in its ohmic region (-50 mV/50 mV) and N=1 (see above). Note that ΓTVS is equivalent to the geometrical average of the SEL values (Γ1Γ2) 1/2 .31, 32

Figure S14 .

 S14 Figure S14. TVS plot from the mean Ī-V traces (Fig. 2 in main text) for the four molecular junc'ons.

Figure S15 .Figure S17 .

 S15S17 Figure S15. (a) First deriva've (with Savitzky-Golay smooth on 20 data points) of the mean Ī-V to determine the zero-bias conductance (horizontal dashed lines)

  ). Albeit this number can be es\mated using mechanical models of

			KP5W30	NH4P5W30	HP5W30	TBAP5W30
		log-Ī	-11.18	-10.24	-9.75	-9.35
	+0.5 V	Ī (A)	6.6x10 -12	5.7x10 -11	1.8x10 -10	4.5x10 -10
		log-σ	0.71	0.76	1.03	1.0
		log-Ī	-11.13	-10.24	-9.77	-9.4
	-0.5 V	|Ī| (A)	7.4x10 -12	5.7x10 -11	1.7x10 -10	4.0x10 -10
		log-σ	0.65	0.78	1.13	1.0
		log-Ī	-9.89	-9.00	-8.97	-8.89
	+1 V	Ī (A)	1.3x10 -10	1.0x10 -9	1.1x10 -9	1.3x10 -9
		log-σ	0.77	0.84	0.98	1.06
		log-Ī	-9.73	-9.50	-9.21	-8.95
	-1 V	|Ī| (A)	1.9x10 -10	3.2x10 -10	6.2x10 -10	1.1x10 -9
		log-σ	0.83	0.92	1.10	1.13

Acknowledgements.

We acknowledge support of the CNRS (France), project neuroPOM, under a grant of the 80PRIME program. We are thankful to Dr. Masahiro Sadakane for sharing his knowledge about the preparaIon of the tetrabutylammonium salt of the Preyssler anion.

Note

The authors declare no compeIng financial interest.

IPCM Jasco Infra-Rouge

Sev IPCM

determined by the SEL model fit and the TVS method.

Supplementary figures for discussion


Julien Diaz
email: julien.diaz@inria.fr

Patrick Joly
email: patrick.joly@inria.fr

ROBUST HIGH ORDER NON CONFORMING FINITE ELEMENT FORMULATION FOR TIME DOMAIN FLUID-STRUCTURE INTERACTION

Keywords: Fluid-structure interaction, high order finite elements, non-conforming meshes, local time-stepping

In this paper we present various numerical method for solving time-dependent fluid-structure interaction problem in two or three dimensions that we claim to be efficient, robust and highly accurate. These methods, based on mixed variational formulations are explicit and conservative and can be of arbitrary high order in space. Their accuracy will be illustrated via a comparison with analytical solutions in simple configuration.

Introduction

Fluid-structure interaction problems appear in numerous and various physical applications such as :

• elastic wave propagation for seismics [START_REF] Komatisch | Wave propagation near a fluid-solid interface : a spectral element approach[END_REF] (sea bottom) ; • ultra-sound waves in biological tissues, such as the bone 7 :

• vibrations of solid structures (a plane, a car, a musical instrument) immersed in surrounding fluids [START_REF] Hamdi | A displacement method for the analysis of vibrations of coupled fluid-structure system[END_REF][START_REF] Derveaux | Modélisation numérique de la guitare acoustique[END_REF] :

An abundant literature has been devoted on the development of efficient numerical methods for solving this type of problem in the frequency domain. Much less attention has been paid to time-dependent problems that introduce additional difficulties, in particular concerning the stability of the time-stepping procedure. Moreover, one is more and more led, in scientific computing, to consider large scale problems as domain decomposition problems for which the use of independent (and non-conforming) meshes in each domain is desired. For time-dependent problems the question is the use of local time-steps.

The objective of this paper is to present various numerical methods for the treatment of time dependent fluid-structure interaction problem. These methods are in the spirit of the meth-ods discussed in [START_REF] Joly | Variational Methods for Time-Dependent Wave Propagation Problems[END_REF] and in the continuation of several works at team ONDES at INRIA [START_REF] Tsogka | Modélisation mathématique et numérique de la propagation des ondes élastiques tridimenssionnelles dans des milieux fissurés[END_REF][START_REF] Fauqueux | Eléments finis spectraux et couches absorbantes parfaitement adaptées pour la propagation d'ondes élastiques en milieu transitoire[END_REF][START_REF] Rodríguez | Une Nouvelle Méthode de Raffinement de Maillage Spatio-Temporel pour l' Équation des Ondes[END_REF][START_REF] Derveaux | Modélisation numérique de la guitare acoustique[END_REF] . They share the following properties :

• their construction is based on (mixed) variational formulations of the complete system ; • they remain explicit except maybe at the interface ; • these are conservative methods : the stability is ensured through the conservation of an energy ; • they can be of arbitrary high order in space.

The outline of the article is as follows. In section 2 we motivate the practical interest of using in this context non-conforming grids and local time-step. In section 3 we present the two main variational formulation that will constitute the basis of our methods : the so-called primal-primal and dual-dual formulations. For each formulation, in section 4, a mixed finite element approximation is proposed for the space discretization. In section 5, the crucial question of the time discretization is considered with an emphasis on the stability analysis. Finally, in section 6, we present numerical experiments that aim at validating (section 6.1) our methods and give an example of possible applications (section 6.2).

Motivation of the work

Choice of the formulation of the problem

It is well-known that an acoustic medium can be regarded as the limit of an isotropic elastic medium when Lamé's coefficient µ tends to zero. As a consequence, when one has a computational code for elastodynamics equations with variable coefficients, it is natural to try to treat the acoustic medium as a particular region in which µ is zero (or very small) (see for instance [START_REF] Bossy | Evaluation ultrasonore de l'os cortical par transmission axiale : modélisation et expérimentation in vitro et in vivo[END_REF] ). There are several disadvantages to this :

• for computational codes based on H 1 finite elements for the displacement field, taking µ = 0 makes lose the H 1 coercivity result (only the L 2 norm of the displacement field is controlled) ; • for the discrete system (or equivalently µ very small) parasitic very small S-waves may be generated ; • a full stress tensor is computed in the fluid where a single pressure is a priori needed.

That is why we shall choose to treat here the fluid-structure interaction problem as a transmission problem between a system in velocity-pressure (v f , p) in the fluid and velocity-stress (v s , σ) in the solid as it is done for instance in [START_REF]Computational Methods for Fluid-Structure interaction[END_REF] .

Interest of non-conforming meshes and local time-stepping

With standard numerical methods, the choice of the discretization parameters (the space step h and the time step ∆t) is driven by the data of the experiment to be simulated. The two essential ones are the typical wavelength λ and the wave speed c. Let us explain below joly Robust High Order Non Conforming Finite Element Formulation for Time Domain Fluid-Structure Interaction 3 how this works in practice and why one is encouraged to use different space steps and time steps in the fluid and the solid media :

(1) for the space discretization, the usual criterion is to take h ≤ λ/N , where N is a given number of points per wavelength which depends on the considered scheme, the desired accuracy, the integration time... In a fluid structure interaction problem, at least three wavelengths have to be taken into account, namely λ f for the fluid, λ S for the shear waves and λ P for the pressure waves. Assume for simplicity that the number N required for the desired accuracy is the same for fluid and solid media and that one is in the typical situation where λ f < λ S < λ P . Taking h f = λ f /N = h s = λ S /N would ensure the same accuracy as h f = h s = λ f /N and would require less CPU time and memory ;

(2) for stability reasons, when one uses explicit schemes in time, as it is our objective, the time step is constrained by a so-called CFL (Courant-Friedrichs-Levy) condition :

c∆t h ≤ α,
where α depends on the numerical method. On the other hand, for a lot a numerical schemes, choosing c∆t/h too small increases needlessly the computational time and induces numerical dispersion phenomena, which could deeply affect the accuracy of the method. If one chooses the same time step in both solid and fluid media, it will satisfy the condition

∆t ≤ min α f h f c f , α s h s c P .
In a situation where c f c P (the typical ratio is one to three), we will have

∆t = α s h s c P α f h f c f .
Therefore using different time steps ∆t f and ∆t s in each grid will give us some flexibility to ensure the balance of the CFL numbers in both fluid and solid media.

About the coupling strategy between fluid and solid

Working with non-conforming meshes in space and time, the main difficulty is reduced to the treatment of the transmission conditions. In what we shall present here, we shall privilege the robustness of the system by imposing the stability of the coupling method. To do that, two basic ingredients are the following :

• base the space discretization on a variational formulation of the coupling system ;

• ensure the stability with the conservation of some appropriate discrete energy.

As far as space discretization is concerned, the well-known mortar element method is a possible approach to deal with the special non-conformity of the problem. However, such method involves the introduction of an additional interface unknown, namely a Lagrange multiplier.

And the discretization of this multiplier implies some constraints, linked in particular to the verification of some inf-sup condition. In this paper we shall privilege the methods that avoids the introduction of Lagrange multipliers.

Obtaining the variational coupled formulation

We denote by d the dimension of the problem (d = 2 or d = 3), x the space variable, t the time, by Ω f the fluid medium and by Ω s the solid medium. Both Ω f and Ω s are supposed to be open sets of IR d such that (see figure 1)

Ω f ∪ Ω s = IR d and ∂Ω f = ∂Ω s = Γ,

the interface between the two media

We lastly denote by n the unit normal vector to Γ, going from Ω s to Ω f .

Fluid

Solid 1. Illustration of the notation

n τ Ω f Ω s Γ Fig.

The variational formulations for the elastic medium

In what follows v s = v s (x, t) denotes the velocity field in Ω s , σ = σ(x, t) the stress tensor, ρ s (= ρ s (x)) > 0 the density and C = ((c ijkl )) (= C(x)) the elastic tensor that satisfy the usual symmetry and positivity relations, almost everywhere in Ω s :

c ijkl = c jikl = c ijlk = c klij and C ξ : ξ ≥ α ξ : ξ for some α > 0. ( 1
)
where σ : ξ denotes the usual scalar products over tensors : σ : ξ = σ ij ξ ij = tr(ξ t σ). The stress-velocity form of linear elastodynamics equations in Ω s is :

∂σ ∂t -C ε(v s ) = 0, (2a)
ρ s ∂v s ∂t -div σ = 0, (2b)
         (S)
where ε(v s ) = 1 2 (∇v s + ∇v t s ) and div is the usual divergence operator for tensors.

We consider below two mixed variational formulations for (S).

The primal formulation. It corresponds to looking for v s in H 1 (Ω s ) = (H 1 (Ω s )) d and σ merely in

L 2 (Ω s ) = {σ | σ ij ∈ L 2 (Ω s )}.
(we do not impose a priori any symmetry conditions) :

Find (v s , σ) : t ∈ IR + → H 1 (Ω s ) × L 2 (Ω s ) such that for all (w s , ξ) ∈ H 1 (Ω s ) × L 2 (Ω s ), d dt Ωs σ : ξ - Ωs C ∇v s : ξ = 0, (3a)
d dt Ωs ρ s v s • w s + Ωs ∇w s : σ - Γ (σn) • w s = 0. (3b)            (SP) Remark 3.1.
It is worthwhile to notice that :

• the equivalence between (S) and (SP) follows from C ε(v s ) = C ∇v s , because of (1) ; • at this level, we do not pay any attention, in equation 2b, to the boundary term which does not make sense in the prescribed functional spaces. This will be treated when taking into account the transmission conditions (see section 3.3).

The dual formulation. It corresponds to looking for v s merely in L 2 (Ω s ) = (L 2 (Ω s )) d and σ in the space of symmetric H div tensor field :

H sym div (Ω s ) = {σ ∈ L 2 (Ω s ) | div σ ∈ L 2 (Ω s ) and σ ij = σ ji }.
We recall here that, thanks to (1) the map ε → C ε is invertible when restricted to symmetric tensors and we denote by A the corresponding inverse map (A also satisfies (1)).

Find (v s , σ) : t ∈ IR + → L 2 (Ω s )×H sym div (Ω s ) such that for all (w s , ξ) ∈ L 2 (Ω s )×H sym div (Ω s ), d dt Ωs A σ : ξ + Ωs div ξ • v s - Γ (ξ n) • v = 0, (4a)
d dt Ωs ρ s v s • w s - Ωs w s • div σ = 0. (4b)            (SD) Remark 3.2.
It is worthwhile to notice that :

• to obtain equation (4a) we have used the fact that div is the formal adjoint of ε( • ) when acting on symmetric tensors fields ; • here again, in equation (4a), the boundary term, which is not a priori well defined, will be treated in section 3.3.

The variational formulations for the fluid medium

In what follows v f = v f (x, t) denotes the velocity field in Ω f , p = p(x, t) the pressure field, ρ f (= ρ f (x)) > 0 the density and c f (= c f (x)) > 0 the sound speed.

The pressure-velocity form of the acoustic equations is :

∂p ∂t + c 2 f ρ f div v f = 0, ( 5a
)
ρ f ∂v f ∂t + ∇p = 0. (5b)          (F)
As for the solid, we can write two mixed variational formulations in the fluid :

The primal formulation.

Find (p, v f ): t in IR + → H 1 (Ω f ) × L 2 (Ω f ) such that for all (q, w f ) ∈ H 1 (Ω f ) × L 2 (Ω f ), d dt Ω f p q c 2 f ρ f - Ω f ∇q • v f - Γ q (v f • n) = 0, (6a)
d dt Ω f ρ f v f • w f + Ω f ∇p • w f = 0. (6b)            (FP)
The dual formulation.

Find (p, v f ): t in IR + → L 2 (Ω f ) × H div (Ω f ) such that for all (q, w f ) ∈ L 2 (Ω f ) × H div (Ω f ), d dt Ω f p q c 2 f ρ f + Ω f q div v f = 0, ( 7a
)
d dt Ω f ρ f v f • w f - Ω f p div w f - Γ p (w f • n) = 0. (7b)            (FD)
Robust High Order Non Conforming Finite Element Formulation for Time Domain Fluid-Structure Interaction 7

The coupling between the fluid and the solid

Let us first recall that the natural transmission conditions on the interface Γ are the continuity of normal velocities and of normal stresses :

v s .n = v f .n, ( 8a
) σn = -p n. (8b)      (TC)
Coupling the primal formulation in one domain with the dual formulation in the other one would involved the introduction of a Lagrange multiplier as for mortar elements. The two other possibilities we are considering here are multiplier free.

The primal-primal formulation We wish to couple (SP) to (FP). We use the continuity of normal stresses condition (8b) in the boundary integral of (3b) :

Γ (σ n) • w s = - Γ (p n) • w s . (9)
Note that the right hand side of ( 9) makes sense since p is searched in H 1 . In the same way we use the continuity condition of normal velocities (8a) in the boundary integral of (6a) :

Γ q (v f • n) = Γ q (v s • n). (10)
Once again the integral of the right hand side of (10) makes sense since q and v s both belong in H 1 . The primal-primal formulation is then :

Find (p, v f , σ, v s ) : t ∈ IR + → X such that for all (q, w f , ξ, w s ) in X, d dt Ω f p q c 2 f ρ f - Ω f v f • ∇q - Γ q (v s • n) = 0, (11a)
d dt Ω f ρ s v f • w f + Ω f ∇p • w f = 0, (11b)
d dt Ωs ρ s v s • w s + Ωs σ : ∇w s + Γ p n • w s = 0, (11c)
d dt Ωs σ : ξ - Ωs C ∇v s : ξ = 0, (11d)
                                   (PP) where X = H 1 (Ω f ) × L 2 (Ω f ) × L 2 (Ω s ) × H 1 (Ω s ).
The dual-dual formulation We wish here to couple (SD) with (FD). Using the transmission conditions is slightly more delicate than for the primal-primal formulation. Let us first explain how we treat the boundary term in (4a).

From the identity

u = (u • n)n + n ∧ (u ∧ n), the dot product (ξ n) • v s in the boundary integral of (4a) is equal to : ((ξ n) • n)(v s • n) + ((ξn) ∧ n) • (v s ∧ n).
Let us also notice that the normal stress continuity condition (8b) is equivalent to :

(σ n) ∧ n = 0, ( 12a
) p = -(σ n) • n. (12b)        (8b) Γ p (w f • n) = - Γ (σ n) • n (w f • n). (14)
Note that there is still a problem of functional analysis since, if we take (σ, ξ) in H sym div,0 (Ω s ) and (v s , w s ) in H div (Ω s ), the two right hand sides of ( 13) and ( 14) appear a priori as integrals of products of functions in H -1 2 (Γ). That is why we shall look for σ and v s in smaller (and not closed) subspaces, namely :

   H div (Ω f ) = {v ∈ H div (Ω f ) | v • n ∈ L 2 (Γ)}, H sym div,0 (Ω s ) = {σ ∈ H sym div,0 (Ω s ) | (σn) • n ∈ L 2 (Γ)},
equipped with their natural Hilbert space norm.

Finally our dual-dual formulation is : joly Robust High Order Non Conforming Finite Element Formulation for Time Domain Fluid-Structure Interaction 9

Find (p, v f , σ, v s ) : t in IR + → X such that for all (q, w f , ξ, w s ) in X, d dt Ω f p q c 2 f ρ f + Ω f q div v f = 0, (15a)
d dt Ω f ρ f v f • w f - Ω f p div w f + Γ (σ n) • n (w f • n) = 0, (15b)
d dt Ωs A σ : ξ + Ωs div ξ • v s - Γ (ξ n) • n(v f • n) = 0, (15c)
d dt Ωs ρ s v s • w s - Ωs w s • div σ = 0. (15d)                                    (DD) where X = L 2 (Ω f ) × H div (Ω f ) × H sym div,0 (Ω s ) × L 2 (Ω s )

Space discretization

We present in this section the finite elements used for the space discretization of the variational formulation. For the dual-dual formulation ,we have chosen , in the solid, finite elements constructed from Tsogka's 1 Q div k+1 -Q k finite elements. In the fluid, we use the classical Raviart-Thomas' finite elements. For the primal-primal formulation we have chosen Fauqueux's 2 Q k finites elements in both fluid and solid.

In the following we denote

• T h f (resp. T hs ) a mesh of Ω f (resp. Ω s ) composed of quadrangles (if d = 2) or hexaedra (if d = 3) K f i (resp. K s i ) ; • F f i (resp. F s i ) : K → K f i (resp. K s i
) is a one to one vectorial mapping transforming the boundary of K onto the boundary of

K f i (resp. K s i ) ; • K is the unit square (or cube) [0 ; 1] d .

Dual-dual formulation

In the solid. Tsogka's finite elements are specifically defined for regular meshes, therefore we suppose that all the elements of T hs are squares. As a consequence, curved interfaces shall be staircase discretized. For a sake of simplicity we restrict ourselves to present the finite element of order 1, Q div 1 -Q 0 in dimension 2 and we refer the reader to 1 for the extension to higher order element or to dimension 3.

The Q div 1 space is the space of stresses which are piecewise Q 1 and satisfy : • σ 11 is discontinuous through the horizontal edges of the mesh and continuous through the vertical ones ;

• σ 22 is discontinuous through the vertical edges of the mesh and continuous through the horizontal ones ;

• σ 12 = σ 21 is a continuous function.
The degrees of freedom of the stresses are all associated to the nodes of the elements (see figure 2). In order to be conforming in the space H sym div,0 we have to impose that σ 12 = 0 for all nodes at the fluid-solid interface.

Note that (σn)n belongs to L 2 (Γ) for all σ in Q div 1 .

The Q 0 space is the space of piecewise constant velocity fields (the degrees of freedom are located at the center of the element -see figure 2). In the fluid. The pressure is piecewise constant (it belongs to Q 0 ) and Raviart-Thomas's finite elements are used for the velocities (see figure 3).

p v f 2 v f 2 v f 1 v f 1
Semi-discretized equations. Applying the standard finite element procedure leads to the following ordinary differential system where the vectors P , V f , V s and Σ are composed of the degrees of freedom of p, v f , v s and σ : The matrices D f and D s represent discrete divergence operators and B and B * discrete trace operators. Since p and v s are piecewise constant, the matrices M P and M s are diagonal. We apply the mass-lumping technique describe in 3 in order to obtain the matrices M f diagonal and M Σ (5 × 5) block-diagonal.

M P dP dt + D f V f = 0, ( 16a
)
M f dV f dt -D * f P -BΣ, = 0 (16b) M Σ dΣ dt + D * s V s + B * V f = 0, (16c)
M s dV s dt -D s Σ = 0. (16d)                               

Primal-primal formulation

In view of using Fauqueux's finite element in the solid, we have to modify the formulation (PP) by introducing an auxiliary function ∂γ ∂t = ∇v s . The formulation, equivalent to (PP), that will be used for the discretization is the following :

Find ((p, v f , σ, v s ), γ) : t ∈ IR + → X × L 2 (Ω s ) such that for all (q, w f , ξ, w s ) in X, d dt Ω f p q c 2 f ρ f - Ω f v f • ∇q - Γ q (v s • n) = 0, (17a)
d dt Ω f ρ f v f • w f + Ω f ∇p • w f = 0, (17b)
d dt Ωs ρ s v s • w s + Ωs σ : ∇w s + Γ p n • w s = 0, (17c)
Ωs σ : ξ - Ωs C γ : ξ = 0, (17d)
d dt Ωs γ : ξ - Ωs ∇v s : ξ = 0, (17e)
                                             (PP') where X = H 1 (Ω f ) × L 2 (Ω f ) × L 2 (Ω s ) × H 1 (Ω s ).
The H 1 functions p and v s are approximated by continuous functions whose respective restrictions p i and v j s to K f i and K s j are such that

p i = p i • F f i and v j s = v j s • F s j are Q k in K.
Moreover, in view of mass lumping, the degrees of freedom associated to p i and v j s are chosen to coincide with the Q k Gauss-Lobatto quadrature points.

The fields v f and σ are approximated by discontinuous fields whose respective restrictions v i f and σ j to K f i and K s j are such that (this is the Piola transform, J l i is the Jacobian of F l i , l = f, s) :

v i f = |J f i |DF f i -1 v i f • F f i and σ j = |J s j |DF s j -1 σ j • F s j are Q k in K . ( 18
)
With obvious notation we end up with the following ordinary differential system ;

M P dP dt -G * f V f -BV s = 0, (19a)
M f dV f dt + G f P = 0, (19b)
M s dV s dt + G * s Σ + B * P = 0, (19c)
M Σ Σ -C Γ = 0, (19d)
M Σ dΓ dt -G s V s = 0. ( 19e
)                                     
where the matrices G f and G s are discrete gradient operators and C approximates the multiplication by tensor C. The mass lumping procedure is once again applied, so that the matrices M P , and M s are diagonal, M f and M Σ are (d×d) block-diagonal and C is (d 2 ×d 2 ) block-diagonal.

Remark 4.1. In the practical computation, the unknown Γ is eliminated and the two equations (19d) and (19e) are replaced by

M Σ dΣ dt -CM -1 Σ G s V s = 0. ( 20
)
Remark 4.2. Let us recall that an interesting property of the choice (18) is

K f i ∇p i • v j f = b K ∇ p i • v j f and K s i ∇v i s : σ j = b K ∇ v i s : σ j .
As a consequence, the matrices G f and G s are (quasi) independent on the meshes T h f and T hs (only their topology is involved) : all the information about the coefficients and the geometry of the meshes is contained in the diagonal matrices M P , M f , M s , M Σ and C.

Time discretization

We now concentrate ourselves on the semi-discretized scheme obtained from the primalprimal formulation, but the following can also be applied to the dual-dual formulation.

We intend to obtain second order time-centered schemes. The most natural way to do so is to compute P et V f at shifted instants (n and n + 1 2 for instance -see figure 4), so that (19a) and (19b) could be explicitly solved when the coupling terms cancel out. We proceed in the same way for Σ and V s and compute (19c) and (19d) explicitly (see figure 5). We will first present the time-discretization in the case of a conforming scheme (∆t f = ∆t s ) before considering a particular non-conforming case (namely ∆t f = 2∆t s ).

Using the same time step in both domains

Two choices are possible :

(1) compute V s and P at shifted instants (see figure 6) and obtain a completely explicit scheme ; (2) compute V s and P at the same instants (see figure 7), which leads to an explicit scheme for the internal nodes but not on the interface.

t P P P V f V f Σ Σ Σ V s V s Fig. 6. Time discretization, choice (1) t P
P P V f V f V s V s V s Σ Σ Fig. 7. Time discretization, choice (2)
Fully explicit scheme. We consider choice (1) and apply the following scheme :

M P P n+1 -P n ∆t -G * f V n+ 1 2 f -B V n+ 1 2 s = 0, (21a)
M f V n+ 1 2 f -V n-1 2 f ∆t + G f P n = 0, (21b)
M s V n+ 1 2 s -V n-1 2 s ∆t + G * s Σ n + B * P n = 0, (21c)
M Σ Σ n+1 -Σ n ∆t -CM -1 Σ G s V n+ 1 2 s = 0. (21d)                                    (I)
Since only the inversion of the mass matrices is involved, the scheme is fully explicit. The disadvantage is that the stability condition of the scheme seems to be affected by the transmission conditions (it does not coincide with the stability condition imposed in the two interior schemes -see below).

Semi-explicit scheme. We consider choice (2) and apply the following centered discretization procedure :

M P P n+ 1 2 -P n-1 2 ∆t -G * f V n f -B V n+ 1 2 s + V n-1 2 s 2 = 0, (22a)
M f V n+1 f -V n f ∆t + G f P n+ 1 2 = 0, (22b)
M s V n+ 1 2 s -V n-1 2 s ∆t + G * s Σ n + B * P n+ 1 2 + P n-1 2 2 = 0, (22c)
M Σ Σ n+1 -Σ n ∆t -CM -1 Σ G s V n+ 1 2 s = 0. (22d)                                    (II)
It is obvious from equations (22b) and (22d) that V f and Σ can be computed explicitly, while P and V s are obtained by solving a system of the form :

M P P n+ 1 2 - ∆t 2 BV n+ 1 2 s = F n+ 1 2 , (23a)
M s V n+ 1 2 s + ∆t 2 B * P n+ 1 2 = G n+ 1 2 , (23b)
        
where F n+ 1 2 and G n+ 1 2 are known from previous time steps. The invertibilty of (23a) and (23b) is guaranteed by the fact that

M P + ∆t 2 4 BM -1 s B *
is symmetric positive definite. Moreover, as the matrices M P and M s are both diagonal and as the matrix B only "sees" the interface, the scheme is implicit only for the degrees of freedom of P and V s that are located on the interface. The "boundary matrix" to be inverted is accordingly very sparse.

On the stability analysis of schemes I and II. The stability of both schemes can be analyzed by energy techniques. If n is an integer, we introduce the two discrete energies :

E n+ 1 2 f = 1 2 M P P n+ 1 2 , P n+ 1 2 + M f V n+1 f , V n f and E n+ 1 2 s = 1 2 M s V n+ 1 2 s , V n+ 1 2 s + CΓ n+1 , Γ n ,
where the sequence Γ n is defined by the following equation :

M Σ Γ n+1 -Γ n ∆t -G s V n+ 1 2 s = 0, (24)
joly Robust High Order Non Conforming Finite Element Formulation for Time Domain Fluid-Structure Interaction 15 so that

M Σ Σ n = CΓ n (25) E n+ 1 2 f and E n+ 1 2 s
are respectively discrete equivalent of the continuous energies

E f = Ω f p 2 c 2 f ρ f + Ω f ρ f |v f | 2 and E s = Ωs ρ s |v s | 2 + Ωs C -1 σ : σ = Ωs ρ s |v s | 2 + Ωs C γ : γ.
Each of this discrete energies would be conserved in the absence of coupling. With the coupling scheme II (respectively I) one can show the conservation of the global discrete energy :

E n+ 1 2 II = E n+ 1 2 f + E n+ 1 2 s
, ( respectively E

n+ 3 4 I = E n+1 f + E n+ 1 2 s -∆t(BV n+ 1 2 s , P n+1 ) ).
Therefore, the L 2 stability of the scheme follows from the positivity of each energy.

For scheme II, it is clear that the positivity of E .

For the positivity of E

n+ 1 2 f
, we rewrite the scalar product

M f V n+1 f , V n f as M f V n+1 f + V n f 2 , V n+1 f + V n f 2 - ∆t 2 4 M f V n+1 f -V n f ∆t , V n+1 f -V n f ∆t .
Then equation (22b) leads to

M f V n+1 f -V n f ∆t , V n+1 f -V n f ∆t = G * f M -1 f G f P n+ 1 2 , P n+ 1 2 and E n+ 1 2 f becomes E n+ 1 2 f = 1 2 M f V n+1 f + V n f 2 , V n+1 f + V n f 2 + M P - ∆t 2 4 G * f M -1 f G f P n+ 1 2 , P n+ 1 2 .
Using the same way it is clear that :

E n+ 1 2 s = 1 2 C Γ n+1 + Γ n 2 , Γ n+1 + Γ n 2 + M s - ∆t 2 4 G * s M -1 Σ CM -1 Σ G s V n+ 1 2 s , V n+ 1 2 s . Thus E n+ 1 2

II

will be positive if the matrices

M P - ∆t 2 4 G * f M -1 f G f and M s - ∆t 2 4 G * s M -1 Σ CM -1 Σ G s
are definite positive. It has been shown 2 that the first one was so if

∆t ≤ α r h f c f √ d (26)
and that the second one was so if

∆t ≤ α r h s c P √ d . ( 27
)
where α r depends on the order of the finite element method, as described in 2 , for instance α 5 0.101.

Finally the scheme II is stable (i.e. the energy E 1 2 is positive) if the following stability condition is satisfied :

∆t ≤ α r min h s c P √ d , h f c f √ d . ( 28
)
For scheme I it is clear that condition (28) is a necessary condition for the positivity of E I (consider separately V n+ 1 2 s = 0 and P n+1 = 0). However, one also has to ensure that the additional term ∆t(BV

n+ 1 2 s , P n+1 )
can not be too large, which provides an additional condition on the matrix ∆tB that must not be too large. This leads to an additional CFL-like condition of the type

∆t ≤ α c h s c P √ d . ( 29
)
with α c ≤ α r . In this sense the stability condition of scheme I is affected by the coupling procedure and that is why, for robustness purpose the scheme II will be preferred.

Remark 5.1. The numerical experiments we have achieved seem to show that α c is close to α r if one uses the primal-primal method whereas it could be much less than α r if one uses the dual-dual method.

Local time-stepping

We present here an adaptation of the multiplier free mesh refinement technique [START_REF] Joly | Variational Methods for Time-Dependent Wave Propagation Problems[END_REF][START_REF] Rodríguez | Une Nouvelle Méthode de Raffinement de Maillage Spatio-Temporel pour l' Équation des Ondes[END_REF] to the fluid-structure coupling problem. We describe here the case where ∆t f = 2∆t s = 2∆t (scheme 1-2), but this could be extended to any rational ratio p -q by adaptating the results presented in [START_REF] Rodríguez | Une Nouvelle Méthode de Raffinement de Maillage Spatio-Temporel pour l' Équation des Ondes[END_REF] . Here again several schemes could be consider. We shall present a scheme that possesses the same characteristics and the same robustness properties that the scheme II presented above :

• the scheme is explicit in the interior domain and implicit on the interface ;

• the scheme conserves some discrete energy and the stability condition is not affected by the coupling procedure.

Robust High Order Non Conforming Finite Element Formulation for Time Domain Fluid-Structure Interaction 17

We choose to compute the unknown P at time t 2n = 2n∆t, V f at time t 2n+1 , V s at time t n and Σ at time t n+ 1 2 (see figure 8). In the time interval [t 2n ; t 2n+2 ] we shall use the following system :

• in the fluid :

M P P 2n+2 -P 2n 2∆t -G * f V f 2n+1 -B[V s ] 2n+1 = 0 (30a) M f V f 2n+1 -V f 2n-1 2∆t + G f P 2n = 0, (30b)
           (III.f)
• in the solid (first time step) :

M s V 2n+1 s -V 2n s ∆t + G * s Σ 2n+ 1 2 + B * [P ] 2n+ 1 2 = 0, (30c)
M Σ Σ 2n+ 1 2 -Σ 2n-1 2 ∆t -CM -1 Σ G s V 2n s = 0 (30d)            (III.s1)
• in the solid (second time step) :

M s V 2n+2 s -V 2n+1 s ∆t + G * s Σ 2n+ 3 2 + B * [P ] 2n+ 3 2 = 0 (30e) M Σ Σ 2n+ 3 2 -Σ 2n+ 1 2 ∆t -CM -1 Σ G s V 2n+1 s = 0 (30f)            (III.s2) where [V s ] 2n+1 , (respectively [P ] 2n+ 1 2 and [P ] 2n+ 3
2 ) represents some approximations of V s at time 2n + 1 (respectively P at time 2n + [START_REF] Tsogka | Modélisation mathématique et numérique de la propagation des ondes élastiques tridimenssionnelles dans des milieux fissurés[END_REF] 2 and at time 2n + 3 2 ) that remain to be define. This terms are determined to get an energy conservation result. An energy can be deduced in the fluid :

t t P V f P V s V s V s Σ Σ
E 2n f = 1 2 M P P 2n , P 2n + M f V 2n+1 f , V 2n-1 f ,
and in the solid :

E 2n s = 1 2 M s V 2n s , V 2n s + CΓ 2n+ 1 2 , Γ 2n-1 2 ,
where the sequence Γ n+ 1 2 is defined by the following equation :

M Σ Γ n+ 1 2 -Γ n ∆t -G s V n s = 0. (31)
Denoting E 2n = E 2n f + E 2n s (that is obviously a positive energy under the condition (28)), it can easily be shown that :

E 2n+2 -E 2n 2∆t = B * [P ] 2n+ 3 2 V 2n+2 s + V 2n+1 s 4 + B * [P ] 2n+ 1 2 V 2n+1 s + V 2n s 4 -B [V s ] 2n+1 P 2n+2 + P 2n 2
As a consequence a conservative scheme is obtained by choosing :

[P ] 2n+ 1 2 = [P ] 2n+ 3 2 = P 2n+2 + P 2n 2 (32) and [
V s ] 2n+1 = V 2n+2 s + 2V 2n+1 s + V 2n s 4 . (33)
The system (30) then becomes :

• in the fluid :

M P P 2n+2 -P 2n 2∆t -G * f V f 2n+1 -B V 2n+2 s + 2V 2n+1 s + V s 2n 4 = 0 (34a) M f V f 2n+1 -V f 2n-1 2∆t + G f P 2n = 0 (34b)            (III.f)
• in the solid (first time step) :

M s V 2n+1 s -V s 2n ∆t + G * s Σ 2n+ 1 2 + B * P 2n+2 + P 2n 2 = 0 (34c) M Σ Σ 2n+ 1 2 -Σ 2n-1 2 ∆t -CM -1 Σ G s V s 2n = 0 (34d)            (III.s1)
A post-treated scheme. Clearly, the transmission conditions (32) are only first order accurate. A consequence of that is the appearance of parasitic evanescent waves confined at the interface. To overcome this difficulty, in 6 the following post-treatment is proposed.

• the standard mean of the unknowns over two time steps in the fluid :

P 2n+1 = P n+2 + P n 2 and V 2n f = V 2n+1 f + V 2n-1 f 2
• the following weighted mean of the unknowns over three time steps in the solid :

V 2n s = V 2n+1 s + 2V 2n s + V 2n-1 s 4 and Σ 2n+ 1 2 = Σ 2n+ 3 2 + 2Σ 2n+ 1 2 + Σ 2n-1 2 4
Then scheme III can also be reinterpreted as follows :

• in the fluid :

M P P 2n+1 -P 2n-1 2∆t -G * f V 2n f -B V 2n+1 s + V 2n-1 s 2 = F 2n , (37a)
M f V 2n f -V 2n-2 f 2∆t + G f P 2n-1 = 0, (37b)
           (III'.f)
• in the solid (first time step) :

M s V 2n s -V 2n+1 s ∆t + G * s Σ 2n-1 2 + B * P 2n+1 + 3P 2n-1 4 = 0, (37c)
M Σ Σ 2n-1 2 -Σ 2n-3 2 ∆t -CM -1 Σ G s V 2n-1 s = 0, (37d)
           (III'.s1)
• in the solid (second time step) :

M s V 2n+1 s -V 2n s ∆t + G * s Σ 2n+ 1 2 + B * 3P 2n+1 + P 2n-1 4 = 0, (37e)
M Σ Σ 2n+ 1 2 -Σ 2n-1 2 ∆t -CM -1 Σ G s V 2n s = 0. (37f)            (III'.s2)
Remark 5.2. The new scheme has the same stability property as scheme III, however we remark that the coupling terms are this time second order accurate, that results into a gain of accuracy.

Numerical Results

In this last section we present various numerical results that first aim at showing that the methods proposed in this paper provide accurate numerical solutions (via the comparison with exact solutions), in particular for interface phenomena (section 6.1). Our next results

Robust High Order Non Conforming Finite Element Formulation for Time Domain Fluid-Structure Interaction 21 aim at proving that our methods are able to handle efficiently complex problems such as the propagation of ultra-sounds waves in the bone (section 6.2).

Comparison with analytical solutions

In two dimensions. In order to validate both the primal-primal and the dual-dual code in two dimensions, we have compared the solution obtained with both methods to an analytical one obtained by the Cagniard-de Hoop technique [START_REF] De Hoop | The surface line source problem[END_REF][START_REF] Poncelet | Réponse transitoire à une source impulsionnelle en espace en milieu anisotrope. Calcul de fonctions de Green par la méthode de Cagniard-de Hoop[END_REF] in a simple configuration : two half-planes separated by a straight interface, one constitutes the fluid medium and the second one constitutes the solid medium. Since this is not the main purpose of this article, we shall not give here neither the details of the calculation of this analytical solution (which are rather technical), nor the expression of this solution (which is rather complicated and has no particular interest in the context of this paper). Let us simply mention that the Cagniardde Hoop method provides a closed form explicit expression of the Green's function of the problem associated to a point source in the fluid. Then, a general solution can be obtained by a simple convolution with the source term that can be computed with a Matlab code with an arbitrary accuracy. In the experiments we consider, the physical parameters for the fluid are c f = 1500 m s -1 and ρ f = 1000 kg m -3 and the physical parameters for the solid are c S = 1800 m s -1 , c P = 4000 m s -1 and ρ s = 1850 kg m -3 . For the numerical computation the size of the computational domain in each medium is 20 mm × 5 mm and Perfectly Matched Layers are used to simulate infinite half spaces. The source function is a point source in space, (located in the fluid at 2mm from the interface -see figure 9) whose time variation is given as the first derivative of a Gaussian whose dominating frequency is 1 MHz.

We compare two numerical calculations : • for the dual-dual methods we use Q div 1 -Q 0 finite elements, regular meshes in each domain but non-conforming : the space step is 0.1 mm in the fluid and 0.12 mm in the solid which corresponds to 15 points of calculation per wavelength in both fluid and solid (for S-wave). We use the same time step in the two grids chosen in order to satisfy CFL condition in both the fluid and the solid : ∆t = 45 ns ;

• for the primal-primal method we use Q 5 finite elements. Here the global mesh is regular and conforming, the space step being equal to 1 mm. Conversely we use different time step in the fluid (∆t f = 35.6 ns) and in the solid (∆t s = 17.8 ns) so that both CFL conditions are optimally satisfied.

We use two types of representation for the numerical results :

• snapshots at a given time (namely t = 4µs) of the modulus of the velocity field (see figures 10, 11 and 12). The results presents a very good agreement between the three solutions ; • to get a better idea on the accuracy, we also represent the variation in time of the two components of the velocity field at a receiver located inside the solid whose precise position is indicated in figure 9. The numerical solution obtained by the primal-primal method (the red curve) is compared to the analytical one (the blue curve) in figures 13 and 14 (we have represented the solution computed by the primal-primal method but the dual-dual one presents the same accuracy). The curves are perfectly superimposed, showing once again the accuracy of the method. In three dimensions. A similar experiment has been achieved in the case of two homogeneous half spaces in order to validate the primal-primal code with Q 5 finite elements in three dimensions (the three dimensions version of the dual-dual code is still not existing). The fluid medium is such that c f = 1 m s -1 and ρ f = 1 kg m -3 and the solid medium is such that c S = 2 m s -1 , c P = 4 m s -1 and ρ s = 1 kg m -3 . The size of each domain is 10 m × 5 m × 10 m. The source function in space is a point source, located in the fluid at 0.66 m from the interface -see figure 15 chosen in the fluid and in the solid. Once again, to show that the coupling does not affect the accuracy, the space step is chosen by using the same criterion as for a single fluid or a single solid : h = 0.66 m and the time-step is chosen to satisfy the CFL condition in both the fluid and the solid : ∆t = 9 ms.

Snapshots of the modulus of the velocity field obtained by the numerical method at time t = 3 s are represented in figure 16. The analytical solution obtained by the Cagniard-de Hoop method is represented in figure 17. We also have computed the first and the second components of the velocity in the fluid along the time thanks to a receiver located as shown in figure 9. The numerical solution (the red curve) is compared to the analytical one (the blue curve) in figures 18 and 19. One observes a similar precision as the one obtained in two dimensions. The slight deterioration of the results after t = 3 s have to be attributed to the fact that the parameters into the PML domain have not been chosen sufficiently well.

More realistic experiments

The propagation of ultra-sounds through cortical bones, which is of particular interest for the diagnoses of osteoporosis 7 , can be modelized via a fluid-structure interaction problem. The skin, the muscles and the marrow behaves merely like water while the bone behaves like an anisotropic solid medium. The following experiments are applications of the primal-primal method to such problems, first in two dimensions, then in three dimensions.

In two dimensions. We have first computed the propagation of a wave through a section of a bone. The figure 20 represents the configuration of the experiments :

• the inner disk is the marrow and is considered as filled with water (c f = 1500 m s -1 and ρ f = 1000 kg m -3 ), its radius is 5 mm ; • the circlet is the bone regarded as an isotropic solid with c S = 1700 m s -1 , c P = 3400 m s -1 and ρ s = 1850 kg m -3 (the anisotropy only concerns the third dimension). The external radius of the circlet is 10 mm ; • the exterior medium that represents the muscles and the skin, is also assimilated to water. • the computational domain is limited with the help of PML to simulate an infinite domain.

The mesh of this domain is represented in figure 21. The source function in space is a line Both 2D and 3D of numerical results (that have been compared -with a good agreement -to the solution obtained with finite difference code 7 ) give a good idea of the complexity of wave propagation inside bones.

11 σFig. 2 .

 112 Fig. 2. Finite element in the solid : the stresses (left and center) and the velocities(right)

Fig. 3 .

 3 Fig. 3. Finite element used in the fluid

jolyf 2 Fig. 4 .s 2 Fig. 5 .

 2425 Fig. 4. Time discretization in the fluid

Fig. 8 .

 8 Fig. 8. Time discretization, scheme III

Fig. 9 .

 9 Fig. 9. Configuration of the experiment

Fig. 10 .Fig. 14 .

 1014 Fig. 10. Dual-dual formulation Fig. 11. Primal-primal formulation Fig. 12. Analytical Solutions

Fig. 15 .

 15 Fig. 15. Configuration of the experiment

Fig. 16 .Fig. 19 .

 1619 Fig. 16. Primal-primal formulation Fig. 17. Analytical Solution

July 6 ,Fig. 20 ..

 620 Fig. 20. Configuration of the experiment Fig. 21. Mesh

Fig. 22

 22 Fig. 22. t=2 ms Fig. 23. t=4 ms

Fig. 26 .

 26 Fig. 26. The solution upon a 3D view

  Robust High Order Non Conforming Finite Element Formulation for Time Domain Fluid-Structure Interaction 11

	July 6, 2004 17:41 WSPC/130-JCA	joly

If we strongly impose (12a) in the functional space for tensors, namely :ξ ∈ H sym div,0 (Ω s ) = σ ∈ H sym div (Ω s ) | (σ n) ∧ n = 0 ,the normal velocity condition (8a) can then be used in the boundary integral (4a) :Γ (ξ n) • v s = Γ ((ξ n) • n)(v s • n) = Γ ((ξ n) • n)(v f • n).(13)In the same way, for the boundary integral of (7b), as soon as σ is searched in H sym div,0 (Ω s ), we can use (12b) to write :