Downstream of the vocal folds, it is often assumed that all the kinetic energy of the glottal jet is dissipated by turbulence. Although this assumption is reasonable for normal voicing, there are many configurations for which a constriction or even an occlusion of the vocal tract occurs and leads to a significant pressure recovery. Such a phenomenon occurs during the production of vowel-bilabial plosive sequences. In this paper, we will first present results from measurements on a human speaker during production of vowel-bilabial plosive-vowel sequences. These measurements concern intra-oral pressure, acoustical pressure, EGG and lips video recording of a female subject. The relative timings between the motion of the lips, the glottal signals and the intra-oral pressure are analyzed and related with Voice Onset/Offset Times. A theoretical model for air flow in the lips is validated by experiments on a mechanical replica of the vocal tract including a self-oscillating glottis and a controlled moving constriction to mimic the closure of the lips. This theoretical model is used to achieve numerical simulations of an /apa/ utterance. The differences between the intra-oral pressure observed in vivo and the simulated one lead to an investigation on the influence of the cheeks expansion during bilabial plosive production.

Introduction

The context of this work is physical modeling of speech production. This approach allows to study the mechanical and aeroacoustic mechanisms governing speech production. The goal is to understand the physical phenomena but also to achieve modeling in order to predict the consequences. For speech, numerical simulations based on these theoretical models are used to synthesize sound signals. The development of such a tool has applications in many fields (medical, communication, teaching ...). Speech sounds are produced by various disturbances of the flow of air through the vocal tract. The present study focuses on the bilabial plosives, especially on the air flow in the upper part of vocal tract, downstream to the vocal folds. Many studies on the fluid mechanical interactions in the vocals folds have been published. They propose theoretical air flow models which are useful to model voiced sound. In this paper, we investigate on the air flow in the constriction shaped by lips, during the production of a bilabial plosive. The physical modeling of the flow in the whole vocal tract provides prediction of the evolution of the intra-oral pressure, which is, in bilabial plosive utterance, responsible for the creation of an audible acoustic disturbance. Otherwise, we study the effect of the expansion of oral cavity volume under the intra-oral pressure effect, during bilabial plosive production. The physical origin of this phenomenon is related to the elasticity of the cheeks. This deformation of the cheeks causes a variation in the volume of the oral cavity and therefore in the intra-oral pressure and in oral flow. In the first part, we present in vivo measurements during /apa/ utterances. Then three different theoretical air flow models are exposed and compared to results from in vitro measurements. Experimental data allows to evaluate the relevance of this models for the air flow in the lips during the production of vowel-plosive-vowel. Numerical simulations of an /apa/ sequence using the air flow model validated before, are achieved and compared to in vivo measurements. Finaly, the influence of the cheeks expansion during bilabial plosive production is identified and quantified through in vivo and in vitro measurements.

In vivo measurements

Experimental setup

In vivo measurements for /apa/ sequences presented in this paper are part of a data base initially achieved by [START_REF] Maillou | Approche physique de la voix d'une personne sourde[END_REF]. These measurements were performed using EVA TM system and a video camera. The EVA device [START_REF] Giovanni | Analyse objective des dysphonies avec l'appareillage eva, état des lieux[END_REF] allows to perform aeroacoustic measurements during speech production. These measurements concern 4 synchronized signals :

• acoustic pressure, recorded by a microphone placed in front of the mouth of the subject,

• electroglottographic signals, measured by two electrodes placed on the larynx,

• Intra-Oral Pressure (IOP) measured by means of a tube of 5 mm diameter, inserted through the corner of the lips,

• lip parameters derived from the video recording of the subject face.

The pressure sensor of EVA have been calibrated in static and dynamic configurations.

Data analysis

The figure 1 shows a typical measurement observed during vowel-bilabial plosive-vowel production. IOP rises significantly, when the lips close, until a maximum value of 1000 Pa. At the lips opening, IOP drops quickly to atmospheric pressure. At this moment, a burst, which is audible as a bilabial plosive, appears in the acoustic pressure signal. As IOP falls, electroglottographic and acoustic pressure signals shows that the vocal folds start to oscillate again to produce the second vowel.

3 Theoretical models and experimental validation

Theoretical models

The objective of this section is to validate a theoretical model for the air flow in the upper part of the vocal tract, during a plosive production. The constriction shaped by the lips is modeled a rectangular duct of length L (in the flow direction), of aperture h l and of width W. We make the assumption of a one-dimensional flow, u = u(x, y, t) e x , where u is the air velocity field and e x is the unit vector parallel to the vocal tract axis. Thus, pressure are supposed to be uniform in a cross section normal to the flow direction. p(x) represents pressure relative to the atmospheric pressure. According to the Reynolds number, based on the aperture h l ,the air flow is expected to be laminar. Moreover, we make the assumption of a quasi-steady flow (low Strouhal The pressure p(x > 0) inside the lip constriction is predicted from the pressure P 1 = p(x = 0) in the vocal tract, upstream of the lips. Three different theoretical models are under investigation for the flow in the lips :

• the steady Bernoulli equation, (Bernoulli),

• boundary layer model for steady flow, using linear velocity profile (Van Zon),

• the steady Bernoulli equation corrected for viscous pressure losses, given by the lubrification theory of Reynolds (Reynolds).

Details of these models, can be founded in [START_REF] Deverge | Influence of collision on the flow througgh in-vitro rigid models of the vocal folds[END_REF].

Inviscid Bernoulli equation

The steady and inviscid Bernoulli equation links the pressure difference P 1 across the lips to the volume velocity Φ l or velocity u = Φ l /(LW), inside the lips.

p(x) = P 1 - ρ 2 Φ l Wh l 2 (1)
The vocal tract area is considered to be large compared with the labial area, thus, we assume that the volume velocity in the vocal tract (for x = 0) is negligible. Since the lips constriction is represented by a uniform duct, the volume velocity conservation implies that p(x = L) = p(x) = 0. In this case, the steady inviscid Bernoulli equation will always predict a intra-labial pressure equal to zero.

Boundary layer equations

The boundary layer theory assumes the existence of a frictionless core , far enough from the wall, with a uniform velocity u e (x). Van Zon proposes to uses a linear velocity profile u(x, y) = u e (x)y/δ(x), where δ(x) is the thickness of the boundary layers. The corresponding volume velocity becomes :

Φ l = Wu e (x)(h l -δ(x)), (2)
δ(x) satisfies the non-linear equation :

4δ(x) h l + 9 ln 1 - δ(x) h l + 5δ(x) h l -δ(x) = 6νxW h l Φ l , (3)
where ν denotes the kinematic viscosity. Equations ( 2) and ( 3) are applied at the exit of the lip constriction (x = L). Then the frictionless core velocity is determined by applying Bernoulli equation ( 1) across the lips, u e (x = L) = 2P 1 /ρ. The set of two equations (2, 3) can be solve to determinate δ(x = L) and Φ l . The value of volume velocity Φ l allows to calculated the viscous boundary layer δ(x) at any position x using equation 3. Finally u e (x) is obtained by using (2) and the pressure is calculated by using Bernoulli equation.

p(x) = P 1 - ρ 2 u 2 e (x) ( 4 )

Steady Bernoulli equation corrected for friction

The principle of this model is to add to inviscid Bernoulli equation, a term ΔP ν (x) accounting off viscous losses. The Bernoulli equation becomes :

p(x) = 1 2 ρ Φ l Wh 2 -P 1 + ΔP ν (x) ( 5 )
For a uniform duct the lubrification theory of Reynolds yields :

ΔP ν (x) = 12ρνΦ l Wh 3 l x (6)
Combination of equations ( 5) and (6) gives a quadratic equation for Φ l which is solved for x = L. The pressure profile p(x) inside the lips is obtained by using equation (5).

Experimental Setup

Validation measurements were performed on a 1:3 scale replica of the human phonatory system. Figure 2 schematises the experimental setup.

The replica is composed of various elements:

• a pressure tank, cuboid about 0.6 m 3 volume, supplied by a compressor,

• a replica of the vocal folds made of latex tube filled with water, able to self-oscillate,

• a replica of lips in metal, whose motion of the upper part is controlled either by a motor or manually,

• a metal uniform tube of length 12 cm with an internal diameter of 2.5 cm for modeling the trachea, (e) : metal replica of lips.

• a uniform tube of length 24.5 cm in metal, including a plexiglas tube of 17 cm length, with the same diameter as the trachea, to represent the vocal tract.

The dimension of the straight channel is L = 2 cm and W = 3 cm. Three differential pressure sensors are used to measure P 0 , P 1 = p(x = 0) and P 2 = p(x = L/2) which correspond respectively, to subglottal, intra-oral and intralabial pressure. Without oscillations, the initial aperture of the glottis h g is 0.7 mm. The aperture h l varies from 0 to 1.2 mm. The motion imposed by a motor is sinusoidal for a frequency of 0.5 Hz.

Theoretical models are evaluated in two different experimental conditions : in case of steady flow and in case of an unsteady flow (with vocal folds oscillations).

Results

Steady flow experiments

In figure 3 is presented an example of the comparison between experimental data and the prediction obtained using the theoretical models described in section 3.1.

Pressure P 2 measured inside the replica of lips (for x = L/2) stay finite even when the channel is closed. This result shows that a complete closure can not be achieved because of the roughness of the mechanical lips surface. Although boundary layer model seems to be more accurate, this theory fails to explain the experimental pressure P 2 for h l < 0.3 mm. This result is due to the fact that, below this value of aperture, the predicted boundary layer becomes larger than the aperture. Corrected Bernoulli equation using Reynolds theory reproduces an acceptable behavior of pressure in mechanical lips. However this theory overestimates the value of P 2 . The relative difference with experimental data, estimated for h l = 0 is of order of 20 %. When h l (t) = 0, the corrected Bernoulli equation predicts P 1 /P 2 = 1/2. The inviscid steady Bernoulli equation is not able to reproduce the behavior since this prediction is always P 2 = 0.

Unsteady flow experiments

In this part, the pressure of water within the latex tubes which represent the vocal folds, is increased to allow self- For small values of apertures (h l < 0.3 mm), predictions by Van Zon theory is less accurate than for steady flow experiments, especially during the closure of mechanical lips (0.4 < t < 0.7). The prediction from corrected Bernoulli equation, as in the steady flow case, shows good agreement with measured pressure P 2 but overestimates the value. The relative difference with experimental data is of order of 20 % again. We conclude from this result that the corrected Bernoulli equation for friction, using Reynolds theory seems to be a reasonable approximation for predict air flow in the lips during bilabial plosive production. This theory will be used in numerical simulations for vowel-plosive-vowel utterance.

Numerical simulations

In this part, we achieve numerical simulation of an /apa/ utterance. We use the physical modeling synthesizer developed by [START_REF] Ruty | Modèles d'interactions fluide/parois dans le conduit vocal[END_REF]. Mechanical behavior the vocal folds is modeled by a two-mass model. The aerodynamic part of the simulation is based on a two-dimensional description of the flow taking into account for the formation of a jet at the glottis and the dissipation by turbulence. The flow at the lips is modeled by the corrected Bernoulli equation, using Reynolds theory of viscous flows. The input parameters of the simulation are subglottal pressure and labial parameters (aperture and width). Subglottic pressure is chosen as a constant equal to the maximum value of IOP measured in vivo (1000 Pa).

Labial parameters, h l and W are directly taken from in vivo measurements. Since the video rate was 50Hz only, a polynomial interpolation was performed between each data point. Simulation results are presented in figure 5. The behavior of simulated IOP is qualitatively similar to the experimental one. However, differences can be observed.

dΦ g /dt IOP [Pa] h l [m] (a) (b) (c)
• Simulated onset and offset times of IOP are different for these measured on the subject but these differences might not be relevant since they are within the accuracy given by the video acquisition system (20 ms).

• During the closure of the lips, the simulated increase in IOP is significantly higher than the measured one. A possible explanation for this departure could lie in the inflation and deflation of the subject's cheeks during the vowel-plosive-vowel production. This hypothesis will be test in the next section.

Influence of the cheeks expansion

The objective of this section is to highlight the influence of the cheeks expansion on the intra-oral pressure during the production of bilabial plosives. The assumption made here is that the increase of oral cavity volume causes an additional depression which concurrence the increase of intra-oral pressure. Firstly this phenomenon will be illustrated by in vivo measurements and reproduced in the laboratory.

In vivo measurements

These measurements consist in recording simultaneously acoustic pressure and intra-oral pressure during /apa/ utterances using EVA station. A subject is asked to repeat 10 times the /apa/ sequences in two different ways :

• in a natural way: N,

• by placing hands on the face to prevent the cheeks from inflating : H. Figure 6: Influence of cheeks expansion for in vivo measurements during /apa/ utterances : evolution of IOP evolution at the closure of the lips for a natural condition (N) and with hands on the face to prevent the cheeks from inflating (H). The increase in IOP is modeled, between the square markers, as a straight line.

We want to quantify the differences in IOP time evolution under both conditions. The increase in IOP that appears at the lips closure, is measured using a straight line whose slope is determined by a linear regression. Figure 6 According to table 1, the increase rate in IOP is stronger for the condition H. The expansion of the cheeks has thus a measurable influence in the case of bilabial plosives. In the next part, this experience is reproduced in the laboratory.

In vitro measurements

These measurements are made with the replica of human phonatory system presented in section 3.2. We focus on the evolution of P 1 at the lips closure. The vocal tract is represented either by a rigid plexiglas tube of 16 cm of length or deformable tube made of latex of 10 cm of length and of thickness 0.2 mm. We compare the increase of P 1 with the rigid or the deformable tubes, for a similar lips movement, h l (t)

The increase rate coefficient a P1 is estimated from P 1 , by a linear regression. An example of results is shown in figure 7 and the corresponding increase rate coefficients are presented in table 2 Table 2: Pressure differences ΔP 1 and increase rate coefficients a P1 of P 1 , for the rigid tube and the deformable tube made of latex

In vivo and in vitro increase rate coefficients are not directly comparable because of the difference in duration of the lips closure. This duration is 200 ms for in-vitro measurement (figure 7) and 50 ms in human (figure 6). Variation of IOP slope due to the stiffness of the oral cavity is 51% for in vitro measurements and 77% for in vivo measurements. This difference can be explained by the fact that latex is stiffer than human cheeks.

Conclusion

Acoustical pressure, electroglottographic signals and intraoral pressure have been measured on a female speaker during /apa/ utterances. Experiments on a replica of phonatory system allow to validate the Bernoulli equation corrected for friction, using the Reynolds lubrification theory, for the air flow in the lips. Numerical simulations of aerodynamical phenomenons governing speech production are achieved for a vowel-bilabial plosive-vowel utterance, by using this theoretical model. The comparison between simulated IOP and IOP measured on human highlights the importance of the expansion of supra-glottal cavity volume during the production of a bilabial plosive. The role of the cheeks expansion has been identified and quantified through in vivo and in vitro measurements. This expansion of the walls of oral cavity causes a reduction of increase of the intra-oral pressure at lips closure. The obtained results are stimulating and lead to the conclusion for the need of the physical modeling of the cheeks expansion. For this purpose, the estimation of elasticity parameters of the latex tube which represents the cheeks is a necessary step to validate a elasticity model of oral cavity walls. In addition, the use of a high sampling rate camera for lip parameters extraction could allow to determinate if the timing difference between the simulated IOP and the measured one is due to a limited time precision of the measure or a lack in physical model used for simulation.

Figure 1 :

 1 Figure 1: In vivo measurements for /apa/ sequence, acoustic pressure P a , electroglottographic signal (EGG), intra-oral pressure (IOP) and labial aperture h l

Figure 2 :

 2 Figure 2: Replica of human phonatory system. (a) : pressure tank. (b) : replica of vocal folds. (c): plexiglas tube. (d): motor and mechanical lips displacement sensor.(e) : metal replica of lips.

Figure 3 :

 3 Figure 3: Experimental and theoretical results for steady flow. (a) : aperture h l of the mechanical lips. (b) : Experimental pressure P 2 and theoretical pressure P 2 computed from experimental pressure P 1 using three theoretical models : Bernoulli inviscid equation (Bernoulli), Bernoulli equation corrected for friction (Reynold), boundary layer equations (Van Zon)

  0

Figure 4 :

 4 Figure 4: Experimental and theoretical results for unsteady flow, during vocal folds replica oscillation. (a) : aperture h l of the mechanical lips. (b) : Experimental pressure P 2 and theoretical pressure P 2 computed from experimental pressure P 1 using three theoretical models : Bernoulli inviscid equation (Bernoulli), Bernoulli equation corrected for friction (Reynold), boundary layer equations (Van Zon)

Figure 5 :

 5 Figure 5: Numerical simulation of an /apa/ utterance. (a): h l polynomial interpolation of lips aperture measured in vivo, represented by the crosses. (b): dΦ g /dt time derivative of glottis volume velocity Φ g simulated. (c): Comparison between simulated IOP and IOP measured in vivo.

Figure 7 : 1 Tubes ΔP 1 [

 711 Figure 7: Influence of the expansion of latex tube during in vitro vowel-plosive-vowel utterances for the rigid tube and the deformable tube made of latex, (a): lips aperture h l , (b): evolution of pressure P 1

Table 1 :

 1 illustrates the determination of the increase rate coefficient a iop , presented in table1. A student test made on the 10 repetitions allows to discriminate, at 1% level, the IOP increase between the two conditions. IOP differences ΔIOP and increase rate coefficients a iop of IOP for /apa/ utterance during in vivo measurement for both conditions of production.

	Conditions ΔIOP [Pa] a iop [Pa.s -1 ]
	N	520	11100
	H	740	19600

  .

	1	x 10 -3
									Latex
									Rigid
	0.5
	0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8
	3000
	2500
	2000								Latex
									Rigid
	1500	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8

Acknowledgments

This work is supported by the Rhône-Alpes region by means of CIBLE PhD grant and by ANR Plasmody project.


Strong ground motion in the very near field as an example of the 208 Iiwate-Miyagi earthquake

Hideo Aochi

STRONG GROUND MOTION IN THE VERY NEAR FIELD AS AN EXAMPLE OF THE 2008 IWATE-MIYAGI EARTHQUAKE

AOCHI H. 1,2 (1) Disaster Prevention Research Institute, Kyoto University, Japan (2) Risk Division, Bureau de Recherches Géologiques et Minières, Orléans, France (h.aochi@brgm.fr)

The extreme ground motion in the near field is actually an engineering interest, because important damages are observed even if any standard construction code is applied. Such extreme ground motion is identified by a PGV (Peak Vround Velocity) more than 100 cm/s or a PGA (Peak Ground Acceleration) more than 800 gal (cm/s/s), for instance (Kawase, JAEE, 1998). The 2008 Mw6.9 Iwate-Miyagi, Japan, earthquake is remembered for the maximum recorded acceleration of 3800 gal at the surface sensors of the station IWTH25 (Kik-net, NIED) located above the ruptured thrusting fault plane. An acceleration of more than 1000 gal is observed even at the borehole sensor (280 m depth) at the same station. The PGV also reaches about 100 cm/s. We focus on the fact that an extreme ground motion is found even at depth, and this should be originated from the source process. In particular, the PGV is found at the very beginning of the ground motions, about 2 seconds after the arrival of S wave, as the first main pulse. We study the mechanics of this strong ground motion, supposing that there might have been related to some very local heterogeneity of the source process. We compute the wave propagation from the published kinematic source models, but due to the resolution limit of the inversion model, the strong first pulse is hardly reproduced in terms of its huge amplitude and short duration. Then we try to simulate the beginning of the rupture process (about the first 3 seconds) kinematically and dynamically. We still encounter two difficulties: The station is unfavorably located with respect the S-wave radiation from the supposed fault plane, and the dynamic model is difficult to give a strong heterogeneity in the beginning without affecting significantly the later rupture process.


Vincent Martin
email: vincent.martin@utc.fr

Agnès Drochon
email: agnes.drochon@utc.fr

Odette Fokapu
email: odette.fokapu@utc.fr

Jean-Frédéric Gerbeau
email: jean-frederic.gerbeau@inria.fr.

MagnetoHemoDynamics in the Aorta and Electrocardiograms

Keywords: Blood ow, electrophysiology, Magnetohydrodynamics (MHD), MRI, nite elements

This paper addresses a complex multi-physical phenomemon involving cardiac electrophysiology and hemodynamics. The purpose is to model and simulate a phenomenon that has been observed in MRI machines: in the presence of a strong magnetic eld, the T -wave of the electrocardiogram (ECG) gets bigger, which may perturb ECG-gated imaging. This is due a magnetohydrodynamic (MHD) eect occurring in the aorta. We reproduce this experimental observation through computer simulations on a realistic anatomy, and with a three-compartment model: inductionless magnetohydrodynamic equations in the aorta, bidomain equations in the heart and electrical diusion in the rest of the body. These compartments are strongly coupled and solved using nite elements. Several benchmark tests are proposed to assess the numerical solutions and the validity of some modeling assumptions. Then, ECGs are simulated for a wide range of magnetic eld intensities (from 0 to 20 Tesla).

Introduction

The interaction between the blood ow and a strong magnetic eld generates an induced electric eld in the body. In a Magnetic Resonance Imaging (MRI) machine, this MagnetoHemoDynamic (MHD) eect mainly occurs in the aortic arch, because the blood velocity there is large and roughly orthogonal to the magnetic eld in this location. The induced voltage may alter the electric potential measured in the ElectroCardioGrams (ECG). It has indeed been observed that it makes the T -wave larger [START_REF] Gaey | Alterations in the rat electrocardiogram induced by stationary magnetic elds[END_REF], which in turn may perturb the ECG-gated image acquisition in MRI. This eect gets worse with higher magnetic eld, which is the current trend of MRI (from typically 3 Tesla nowadays for humans, it may reach 10 Tesla in the future). Dierent ways of overcoming this gating problem have been proposed, for instance using VectoCardioGrams (VCG) [START_REF] Chia | Performance of QRS detection for cardiac magnetic resonance imaging with a novel vectorcardiographic triggering method[END_REF], or acoustic cardiac signal [START_REF] Frauenrath | Feasibility of cardiac gating free of interference with electro-magnetic elds at 1.5 Tesla, 3.0 Tesla and 7.0 Tesla using an MR-stethoscope[END_REF], ECG signal processing [START_REF] Tse | 12-lead ECG in a 1.5 Tesla MRI: Separation of real ECG and MHD voltages with adaptive ltering for gating and non-invasive cardiac output[END_REF]. However, a better understanding and prediction of the MHD eect would have many potential benets. First, it may help improve the quality of the ECG, which is not only used for gating but also for clinical purposes, in particular for interventional MRI or MRI stress tests. Second, the MHD artifact might be used to infer some physical properties of the system, using inverse problem techniques. This paper is devoted to the modeling and numerical simulation of this phenomenon, and can be viewed as a rst step in this direction.

The interaction of the aortic blood ow and a magnetic eld has already been investigated in the literature. An overview of the in vitro experiments, in vivo studies as well as MHD calculations in blood ows can be found in [START_REF] Abi-Abdallah | Eects of static magnetic eld exposure on blood ow[END_REF]. Gaey et al [START_REF] Gaey | Alterations in the rat electrocardiogram induced by stationary magnetic elds[END_REF] and Tenforde et al [START_REF] Tenforde | Cardiovascular alterations in Macaca monkeys exposed to stationary magnetic elds: experimental observations and theoretical analysis[END_REF] observed animal ECG alterations, manifested mainly as T -wave elevations. Later, the same group provided some simulations of the voltage and current distributions, and the resulting reduction in blood ow, by means of a two-dimensional nite element analysis in two sections of the ascending and descending aorta with a conductive wall [START_REF] Kinouchi | Theoretical analysis of magnetic eld interactions with aortic blood ow[END_REF]. They studied the case of a steady state ow with an external magnetic eld transverse to the ascending and descending aorta.

More recently, Sahakian et al published a series of papers where the magnetohydrodynamic eect during cardiac MRI was investigated both experimentally [START_REF] Nijm | Characterization of the magnetohydrodynamic eect as a signal from the surface electrocardiogram during cardiac magnetic resonance imaging[END_REF] and by computer simulations, using the Comsol Multiphysics software [START_REF] Nijm | Inhomogeneous human torso model of magnetohydrodynamic blood ow potentials generated in the mr environment[END_REF][START_REF] Nijm | Extraction of the magnetohydrodynamic blood ow potential from the surface electrocardiogram in Magnetic Resonance Imaging[END_REF]. Computer simulations of the MHD eect, with experimental validations, are also proposed in [START_REF] Kainz | Development and validation of a magnetohydrodynamic solver for blood ow analysis[END_REF]. Another group addressed the problem of ECG simulations in the presence of the magnetic eld [START_REF] Gupta | Simulation of elevated T-waves of an ECG inside a static magnetic eld (MRI)[END_REF] using the closed-form solution of Gold [START_REF] Gold | Magnetohydrodynamic pipe ow. part 1[END_REF], which corresponds to a steady, fully developed ow of a conducting uid in a straight cylinder subjected to a transverse magnetic eld.

In this approach, the aortic arch was divided into slices. In each slice, Gold's solution was used, with a dierent angle for the magnetic eld. The MHD signal was computed from these slices via a transfer function and then directly added to the ECG waveform that was taken from the ECGSIM software [START_REF] Huiskamp | The depolarization sequence of the human heart surface computed from measured body surface potentials[END_REF].

In the present study, state-of-art mathematical models are considered in each physical compartment of the problem on a realistic anatomy: bidomain equations for the cardiac electrophysiology; inductionless MHD equations for the aortic blood ow; electrostatic equations in the torso, which includes bones and lungs. In contrast to what is usually done in the literature, these three compartments are strongly coupled. This allows us to evaluate the possible interaction between the electric potential induced in the aorta and the extracellular potential in the myocardium. The present article signicantly enriches the preliminary results presented in [START_REF] Martin | Magnetohemodynamics eect on electrocardiograms[END_REF]. It contains many details about the numerical methods and several new results. Due to the complexity of the model and the multiphysical coupling, a particular eort has been placed on verifying the computer simulations: several new test cases are proposed to assess the MHD solutions and the coupling algorithm with the torso. In addition, a numerical justication of the inductionless approximation of the MHD equations in the blood ow regime is presented. Finally, new numerical simulations of ECG are performed for a wide range of magnetic elds (up to 20 Tesla) and new conclusions are drawn about the link between the ow rate, the magnetic eld and the ECG perturbation.

Models

Let Ω be a domain representing the whole body (Fig. 1) partitioned into three subdomains, the aorta Ω A , the heart Ω H and the complement Ω T , hereafter sometimes called torso for the sake of brevity: Ω = Ω A ∪ Ω H ∪ Ω T , where ω denotes the closure of an open domain ω, and ∂ω its boundary. The boundary of the aorta is assumed to be entirely included in the torso and divided into an inlet, four outlets and a wall:

∂Ω A = Γ inlet ∪ Γ wall ∪ Γ outlets ⊂ Ω T .
The inlet is a planar cut in the ascending aorta. The rst three outlets Γ outleti , i = 1, 2, 3, are respectively a section of the brachiocephalic, the left common carotid, and the left subclavian arteries. The outlet Γ outlet4 is a section of the abdominal aorta (Fig. 7). The heart boundary is included in the torso ∂Ω H ⊂ Ω T , thus the boundary of the torso is partitioned as follows: ∂Ω T = ∂Ω A ∪ ∂Ω H ∪ Γ ext , where Γ ext is the body surface. Throughout the article, n will denote a normal unit vector on given surfaces.

Let us briey review the main stages of the electrical activity of the myocardium. During the depolarization, an electrical activation wavefront propagates throughout the heart. This generates an electric potential dierences in Ω T . The recording of these potential dierences gives the ECG. The rst event seen in the ECG, called the P -wave, corresponds to the depolarization of the atria.

It is not modeled in the present study. Next, the ventricles are depolarized, which corresponds to the QRS-complex in the ECG. The R-wave, which is the largest positive deection in the lead I, is used for MRI gating. This depolarization triggers the mechanical contraction of the myocardium, which ejects the blood from the left ventricle into the aorta (systole). Next, the ventricles return to their rest state (repolarization), which gives the last deection of the ECG, called the T -wave. The peak ow in the aortic arch, and therefore the peak MHD eect, takes place after the S-wave (Fig. 1).

Figure 1. Left: the three subdomains (aorta, heart and rest of the body, including some bones and the lungs), schema of the magnetic eld B and velocity eld u in the aortic arch. Right: sketch of the lead I of the ECG and aortic ow rate.

MHD in the aorta

In the presence of an electric eld E and a magnetic eld B, Ohm's law in the blood is written j = σ bl (E + u × B), where j denotes the current density, σ bl the electrical conductivity of the blood and u the blood velocity. Ignoring the displacement current 0 ∂ t E, the Maxwell-Ampère equation gives another expression of the current density:

j = curl B µ , ( 1
)
where µ is the magnetic permeability, supposed to be that of the vacuum (µ = µ 0 = 4π10 -7 H m -1 ).

Inserting these two expressions in the Maxwell-Faraday equation

∂ t B + curl E = 0, (2)
and coupling them to the incompressible Navier-Stokes equations that describe blood ow in large arteries [START_REF] Formaggia | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF], we obtain the usual MHD equations:

           ρ ∂u ∂t + u • ∇u -η∆ u + ∇p = 1 µ 0 curl B × B, div u = 0, ∂ B ∂t + curl 1 µ 0 σ bl curl B = curl (u × B), (3)
where p denotes the pressure of the blood, η its viscosity and ρ its density. The right-hand side of (3) 1 corresponds to the Lorentz force j × B, with j given by [START_REF] Abi-Abdallah | Eects of static magnetic eld exposure on blood ow[END_REF].

System (3) is typically used for uids with high electrical conductivity, like liquid metals, but is not very suitable for blood ows. Indeed, the magnetic Reynolds number Rm, which is the dimensionless ratio of magnetic advection to magnetic diusion, is extremely small in that case (typically of order 10 -8 during systole, cf. Table A2). Thus, the magnetic eld is almost unaected by the blood ow. In this regime, it is enough to assume that B is equal to B, the external imposed magnetic eld, neglecting the MHD eect on B. In an MRI device, B is mainly homogeneous and constant in time. Thus, if B ≈ B, the Maxwell-Faraday equation [START_REF] Abi-Abdallah | Alterations in human ECG due to the magnetohydrodynamic eect: a method for accurate R peak detection in the presence of high MHD artifacts[END_REF] gives curl E = 0, and there exists an electric potential φ A such that E = -∇φ A . With this expression of E, and neglecting again the variations of B caused by the induction, Ohm's law becomes j = σ bl (-∇φ A + u × B).

Using this relation in the Lorentz force, and writing the conservation of the electrical charges div j = 0, we obtain a simplied form of MHD equations in Ω A , sometimes called inductionless MHD equations,

       ρ ∂u ∂t + u • ∇u -η∆ u + ∇p = -σ bl ∇φ A × B + σ bl (u × B) × B, div u = 0, div (σ bl ∇φ A ) = div (σ bl u × B) . ( 4
)
In the blood ow regime, system (4) is an excellent approximation of (3) for a wide range of B, as will be shown by numerical simulations in Section 4.1. Moreover, it is numerically better conditioned, less expensive to solve, and it provides the electric potential directly, which is the quantity of interest for the ECG. Fluid boundary conditions. A time dependent plug ow is applied on Γ inlet , a no-slip condition is applied on Γ wall , and a pressure is applied on Γ outleti, i=1,..,4 . The outlet pressures p i are obtained by coupling three-element Windkessel models to the Navier-Stokes equations:

p i = π i + R p,i Q i where Q i is the volumetric ow rate, Q i = Γ outlet i u • n dS

Electrophysiology in the heart

We denote by φ H and V m the extracellular potential and the transmembrane potential, and by σ i and σ e the intra-and extracellular conductivity tensors. In the heart Ω H , the electrical activity is governed by the bidomain equations [START_REF] Sachse | Computational Cardiology: Modeling of Anatomy, Electrophysiology, and Mechanics[END_REF][START_REF] Sundnes | Computing the electrical activity in the heart[END_REF]:

   A m C m ∂V m ∂t + I ion (V m , w) -div (σ i ∇V m ) -div (σ i ∇φ H ) = A m I app , -div ((σ i + σ e )∇φ H ) -div (σ i ∇V m ) = 0, (5)
where A m is a constant representing the rate of membrane area per volume unit and C m the membrane capacitance per area unit. The term I ion (V m , w) represents the ionic current through the cell membranes, and I app is a given source term. The dynamics of the ionic variable w and the ionic current I ion are governed by the phenomenological two-variable model proposed by Mitchell and Schaeer in [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF]. The ventricles are subdivided into four dierent regions where the coecient governing the Action Potential Duration takes dierent values. The atria are not considered here, which explains the absence of the P -wave in our ECGs. The initial activation starts from the septum. Further details can be found in [START_REF] Boulakia | Mathematical modeling of electrocardiograms: A numerical study[END_REF]. On the boundary of the heart, the intracellular medium is isolated:

σ i ∇V m • n + σ i ∇φ H • n = 0 on ∂Ω H . (6)

Electric diusion in the body

To obtain the ECG, the electric potential has to be computed on the surface of the body. For that purpose, the electrostatic equation satised by φ T is solved in

Ω T : div (σ T ∇φ T ) = 0, in Ω T , σ T ∇φ T • n = 0, on Γ ext , ( 7
)
where the conductivity σ T takes dierent values in the dierent regions of the torso (see Table A3).

Coupling the three models together

To close the system, the continuity of the potential and the normal current are imposed at the interfaces between the aorta and the torso ∂Ω A and between the heart and the torso

∂Ω H σ bl ∇φ A • n = σ T ∇φ T • n, φ A = φ T ,
and

σ e ∇φ H • n = σ T ∇φ T • n, φ H = φ T . ( 8
)
Equations ( 4), ( 5) and ( 7), supplemented with the uid boundary conditions and ( 6), are strongly coupled together by the transmission conditions [START_REF] Frauenrath | Feasibility of cardiac gating free of interference with electro-magnetic elds at 1.5 Tesla, 3.0 Tesla and 7.0 Tesla using an MR-stethoscope[END_REF]. This makes the global model particularly challenging from a scientic computing standpoint.

Remark 2.1 In the literature, it is sometimes assumed that the extracellular medium is isolated, i.e. σ e ∇φ H • n = 0 on ∂Ω H . This assumption considerably simplies the problem since it decouples the heart and the torso. But it does, however, have some impact on the ECG [START_REF] Boulakia | Mathematical modeling of electrocardiograms: A numerical study[END_REF]. In particular, it signicantly increases the amplitude of the ECG. Since our goal is to capture the modication of the amplitude of the signal, we wish to avoid a priori simplications that might aect the ECGs.

Numerical methods

All the equations are discretized with P1 nite elements. This section is devoted to a brief presentation of the time schemes and the coupling algorithm. Let N be a positive integer, T > 0 the nal time and τ = T /N the time step. Given the initial conditions, the discrete solutions u n h ,

p n h and φ n A,h are computed in the aorta, V n m,h and φ n H,h in the heart and φ n T,h in the torso, for n = 1, . . . , N . The MHD model of the aorta ( 4) is solved with a fractional step method: given u n h , p n h and

φ n A,h , compute (u n+1 h , p n+1 h , φ n+1 A,h ) in Ω A such that            ρ u n+1 h -u n h τ + u n h • ∇u n+1 h -η∆ u n+1 h -σ bl u n+1 h × B × B = -∇p n h -σ bl ∇φ n A,h × B, ∆ p n+1 h = ρ τ div u n+1 h , div σ bl ∇φ n+1 A,h = div σ bl u n+1 h × B . ( 9
)
The advection is stabilized with a streamline diusion term. In (9) 1 , the pressure and the potential are treated explicitly, thus making it possible to solve (9) 1 , and then (9) 2 and (9) 3 separately.

During this study, implicit monolithic schemes were also considered. But with the parameters of blood ows, and with the time step needed for the bidomain equations, algorithm [START_REF] Gaey | Alterations in the rat electrocardiogram induced by stationary magnetic elds[END_REF] proves to be stable and accurate. Since it is by far the most ecient, it is the algorithm used in the following.

The electrical heart model ( 5) is discretized in time by combining a second order Backward Dierentiation Formula (BDF2) implicit scheme with an explicit treatment of the ionic current, as proposed in [START_REF] Boulakia | Mathematical modeling of electrocardiograms: A numerical study[END_REF].

The interface conditions (8) strongly couple equations ( 4), ( 5) and ( 7) together. This coupling is performed using a Dirichlet-Neumann domain decomposition algorithm, the interface linear problem being solved by a matrix-free GMRES algorithm. The matrix-vector product of the interface problem is computed as follows: at each time step, the potentials coming from the aorta and the heart are sent to the torso where a Dirichlet problem is solved; then the torso sends back the variational form of the electrical current to the aorta and to the heart where a Neumann problem is solved; these iterations are repeated until convergence before the next time step. The weak coupling strategy essentially consists in doing only one iteration of the domain decomposition algorithm, starting with the homogeneous Neumann problems in the aorta and the heart.

Numerical results

Verication tests

In this section, we consider some tests to verify our MHD solver, the validity of the inductionless approximation of the MHD equations in blood ow regime and the coupling strategy. This conguration is similar to a Hartmann ow [START_REF] Gerbeau | Mathematical methods for the Magnetohydrodynamics of liquid metals[END_REF][START_REF] Moreau | Magnetohydrodynamics[END_REF], however it diers in two ways.

Firstly, the standard Hartmann solutions usually refer to the original MHD equations (3), whereas its inductionless approximation (4) is considered here. Secondly, the static magnetic eld B is not imposed along the usual direction Oy but along Oz, in order to avoid trivial uncoupling.

The velocity can be written as u = [u(y), 0, 0] , while the pressure is linear in x, p = -Gx, and the electric potential only depends on y, φ = φ(y). Thus system (4) is reduced to Coupling with a Poisson equation To verify the multi-domain coupling algorithm, we supplement

-ηu (y) + σB 2 0 u(y) -G = -σB 0 φ (y) φ (y) = -B 0 u (y).
Ω with the domain D = [0, L x ] × [L y , H y ], (H y > L y )
where no uid ows, but where the potential satises a Poisson equation (Fig. 3). This conguration roughly mimics the aorta and the torso. In D, the potential satises -div (σ∇φ) = 0. With the boundary conditions φ(0) = φ 0 and φ(H y ) = φ 0 + δφ (where φ 0 and δφ are two given constants) and using the transmission conditions

=      CB 0 y 3 3 -CB 0 L y y 2 2 + G -2Cη σB 0 y + φ 0 , y ≤ L y , G -2Cη σB 0 (y -H y ) + φ 0 + δφ, y ≥ L y .
Again, the analytical and the numerical solutions are in very good agreement (Fig. 3).

4.1.2.

Comparison with an analytical solution: non-conducting straight tube Gold [START_REF] Gold | Magnetohydrodynamic pipe ow. part 1[END_REF] provides the solution to the classical MHD (3) in the case of a steady ow of an incompressible, viscous, electrically conducting uid through a insulating circular pipe in the presence of a transverse uniform magnetic eld. The analytical solution (mean ow, velocity, electric potential along a diameter) is computed as a series of modied Bessel functions of the rst kind [START_REF] Abi-Abdallah | Eects of static magnetic eld exposure on blood ow[END_REF]. The analytical solutions is approximated by truncated series.

For comparison, we solve the inductionless MHD (4) numerically in a physical conguration close to the blood ows. The domain Ω is a straight cylinder of radius R = 1 cm and oriented along Oz. The viscosity is η = 0.05 Poise and the conductivity is σ = 5.10 -3 S cm -1 . A pressure gradient G = -dp dz = 24 g cm -2 s -2 is imposed between the inlet and the outlet. A static magnetic eld of 0, 10, 100 T is imposed along Ox. The wall is insulating and the potential does not depend on z. The boundary conditions are thus

p(z = L) = LG on Γ L , p(z = 0) = 0 on Γ 0 ,
and

u = 0 on Γ wall , ∇φ • n = 0 on ∂Ω. ( 11
)
Three dierent levels of mesh renement were considered, denoted by M i, i=1,2,3 containing respectively 88 000, 380 000 and 1 500 000 tetrahedra. Table 1 shows the ux, the maximum axial velocity, the transverse dierence of electric potential and the maximum electrical current modulus.

The impact of the magnetic eld on the ow is very small (about 2% of ow reduction between B = 0T and B = 10T ). Fig. 4 and5 show the very good agreement of the velocity, the potential and the electrical current with the analytical solution. These results conrm the validity of the inductionless approximation for MHD in blood ows with B ranging from 0 to 100T . Abi-Abdallah et al [START_REF] Abi-Abdallah | Eects of static magnetic eld exposure on blood ow[END_REF] also reached the same conclusion.

Transient simulations starting from 0 are performed on the mesh M 2 . In Fig. 6, the maximum velocity and the potential are represented as a function of time for B = 0, 10, and 100 T . The stationary analytical solutions are also reported. As expected, the transient values actually tend to the correct numerical stationary state. purposes, we compared the results from [START_REF] Kinouchi | Theoretical analysis of magnetic eld interactions with aortic blood ow[END_REF] with those given by the solution of the 3D MHD model (4) coupled to the diusion in the torso [START_REF] Formaggia | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF], and we obtained identical results, see [START_REF] Martin | MagnetoHemoDynamics in Aorta and Electrocardiograms[END_REF].

B(T ) Ha

Flux(cm 3 s -1 ) Vmax(cm s -1 ) ∆φ(mV ) ||j|| max (A m -2 ) exact MHD 0 0

MHD eects on the ECG

This section contains our main results. The three models MHD in the aorta (4), bidomain in the heart [START_REF] Chia | Performance of QRS detection for cardiac magnetic resonance imaging with a novel vectorcardiographic triggering method[END_REF] and diusion in the torso [START_REF] Formaggia | Cardiovascular Mathematics. Modeling and simulation of the circulatory system[END_REF] are coupled together and solved on a realistic anatomy.

The body mesh contains three regions with dierent conductivities: bones in the torso, lungs and remaining parts. A heart, limited to two ventricles, and a large portion of an aorta are inserted into the body. The aorta comes from [START_REF] Moireau | External tissue support and uid-structure simulation in blood ows[END_REF], all the other parts come from Zygote Media Group, Inc. The three subdomains aorta, heart, body were meshed using ghs3d [START_REF] George | Fully automatic mesh generator for 3d domains of any shape[END_REF], see Table 2 and Fig. 7. The mesh of the heart was chosen so as to be ne enough to capture the front propagation (renement intermediate between R2 and R3 of Table 5 in [START_REF] Boulakia | Mathematical modeling of electrocardiograms: A numerical study[END_REF]). Mesh convergence was checked on the ECG: less than 2% change in the peak value of the ECG was observed with more rened meshes, see Fig. 8 right. The time step is τ = 0.25ms for all computations in order to model the electrophysiology accurately. As observed in [START_REF] Boulakia | Mathematical modeling of electrocardiograms: A numerical study[END_REF], larger time steps produce undesirable spurious oscillations during the QRS complex. This time step was checked to be small enough for the MHD.

The computation was performed over 3 cardiac cycles. The data used for the computations are given in Tables A1A3. The conductivities are taken from [START_REF] Potse | A comparison of monodomain and bidomain reaction-diusion models for action potential propagation in the human heart[END_REF]. The blood in the ventricles is assumed to be still and is treated as a passive diusive domain.

The ECGs are obtained by plotting linear combinations of the values of φ T on 9 standard points of the body surfaces (Fig. 7 right). See [START_REF] Malmivuo | Principles and Applications of Bioelectric and Biomagnetic Fields[END_REF] p. 277 for a denition of these quantities. The ECGs obtained without MHD (B = 0T ) are considered a good reference to evaluate the impact of MHD.

Various improvements could be made in the ECG simulation. In spite of the aws observed in this simulated ECG (see e.g. the discussion in [START_REF] Boulakia | Mathematical modeling of electrocardiograms: A numerical study[END_REF]), we believe that it can be considered reasonable for the purpose of this paper, which is to study the perturbations of the ECG by the MHD eect.

Electrocardiograms with MHD eects

In Fig. 8 (left), we show the ow Q imposed at the inlet of the aorta. We also show the ECG in the same picture, in order to show how the velocity pulse is located in time with respect to the cardiac cycle: we set the pulse so that the ejection of blood in the aorta immediately follows the QRS complex (taken from [START_REF] Despopoulos | Color atlas of physiology[END_REF]). The T -wave occurs at approximately t ≈ 250ms.

The 12-lead ECGs represented in Fig. 9 were obtained by solving the full coupled problem with B = 0, 1, 2, 3, 5, 8 and 10T . With B = 0, the results are similar to those presented in [START_REF] Boulakia | Mathematical modeling of electrocardiograms: A numerical study[END_REF][START_REF] Chapelle | Numerical[END_REF]: the peak during the QRS complex is approximately 1.8mV for the lead I, and the peak during the T -wave reaches approximately 0.8mV . The precordial leads V 3 , V 4 and V 5 are almost unaected by the magnetic eld. In all other leads, the T -wave is altered by the MHD in the aorta, and this eect is amplied when B increases. The MHD eect is particularly strong for the limb leads I, II, III, the augmented limb leads aVR, aVL, aVF and for the precordial lead V 1 . In Table 3, we report for various B the extremal values of the electric potential in the aorta and on the body surface, the peak values of the leads I, II and III, and the maximum electrical current intensity in the aorta. In lead I (Fig. 9, top left), the T -wave peak increases from 0.8mV when B = 0T , to 1.6mV when B = 3T , and up to 3.85mV when B = 8T . Leads II and III are less aected. The maximum value increases quasi-linearly with B, see Fig. 8, right. These results are in relatively good agreement with those of [START_REF] Gupta | Simulation of elevated T-waves of an ECG inside a static magnetic eld (MRI)[END_REF], although the ECG given in their Fig. 1 shows a more pronounced MHD eect than ours for lead II. Dierent experimental values are reported in the literature for the T -wave elevation during MRI examinations: Nijm et al [START_REF] Nijm | Characterization of the magnetohydrodynamic eect as a signal from the surface electrocardiogram during cardiac magnetic resonance imaging[END_REF] measured a T -wave elevation of order 0.2mV for healthy subjects inside a 3T MRI scanner. Abi-Abdallah et al [START_REF] Abi-Abdallah | Alterations in human ECG due to the magnetohydrodynamic eect: a method for accurate R peak detection in the presence of high MHD artifacts[END_REF] found increases of between 0.2 and 0.7mV (about 300% increase) for subjects exposed to a 1.5 Tesla magnet.

0

Wave

B

Ha

φ min A φ max A ∆φ A φ min S φ max S ∆φ S Lead I Lead II Lead III ||j A || max (T ) (mV ) (mV ) (mV ) (A m -2
) R all -4.9

1.9 and of a quantity proportional to the magnetic eld B and to the aortic ux Q:

D i (B) ≈ D i (B = 0) + α i B Q (mV ), (12)
where D i is the i th lead and the coecients α i , i = 1, . . . , 12, are given in Table 4. The equation ( 12) is a good approximation for almost all leads, except some precordial leads, particularly V 3 and V 4 .

V 1 (D 7 ) V 2 (D 8 ) V 3 (D 9 ) V 4 (D 10 ) V 5 (D 11 ) V 6 (D 12
) α i -0.43 10 -3 -0.21 10 -3 -0.10 10 -3 -0.025 10 -3 0.10 10 -3 0.20 10 -3 Table 4. Proportionality coecients in the numerical relation [START_REF] Gold | Magnetohydrodynamic pipe ow. part 1[END_REF]. α i in 10 cm -1 = mm -1 .

It might be feared that the MHD eect in the aorta induces a perturbation of the electrical activity of the heart, and possibly lead to brillation. According to our computations, this does not happen, even for strong magnetic eld: the extracellular potential in the myocardium is only slightly modied by the aortic electrical activity, see Fig. 10. We also performed a computation using the weak coupling strategy, see Fig. 11 and Remark 2.1. This simplication results in a modication of the amplitude in the QRS complex and in the Twave, but with no clear pattern, as it can either increase or decrease their relative amplitudes, see the leads I and II in Fig. 11 for B = 10T . It is therefore questionable whether the weak coupling simplication is able to assess the MHD eects on the ECG satisfactorily. 4.2.2. MHD eect in the aorta and in the body In Fig. 12, the electric potential on the skin is represented at various time instants for B = 0, 1, 3, 8T . One can see the eects of the electrical activity of the heart particularly during the QRS complex and the T -wave. The MHD eect in the aorta clearly modies the potential around t = 200ms. The velocity inside the aorta is almost unaected by the MHD eect, even for large magnetic elds [START_REF] Martin | Magnetohemodynamics eect on electrocardiograms[END_REF]. This observation is consistent with the results of [START_REF] Abi-Abdallah | Eects of static magnetic eld exposure on blood ow[END_REF] and of Section 4.1.2. The electric potential on the boundary of the aorta is represented in Fig. 13 during the diastole (t ≈ 200ms). The heart activity slightly modies the potential of the aorta (Fig. 13, left). The maximum voltage due to MHD appears, as expected, across the aortic arch, where the velocity is roughly orthogonal to the magnetic eld. The predicted values are of the order of 10 mV to 100 mV , which is in agreement with the calculations of [START_REF] Abi-Abdallah | Eects of static magnetic eld exposure on blood ow[END_REF] and of [START_REF] Kinouchi | Theoretical analysis of magnetic eld interactions with aortic blood ow[END_REF]. The electrical current is represented in Fig. 14 and its maximum value is reported in Table 3 for various B. The current intensity is the largest at the top and at the bottom of the aortic arch, near the aortic wall. This is in agreement with the Gold simulations (Fig. 5). The maximum electrical current density in the aorta is approximately 0.6A m -2 at B = 1T , 1.8A m -2 at 3T , and 4.9A m -2 at 8T . This value is roughly twice as large as in the Gold test case (for Ha = 1 it is ≈ 5.1A m -2 in the aorta whereas it is 2.9A m -2 in the ideal tube, see Table 1). This dierence can be explained by a larger maximum velocity, a more complex geometry (non-innite tube, outlets, etc.), and a coupling with the surrounding tissue (conducting walls). These values can also be compared to the current densities generated by the cardiac activity that reach approximately 1A m -2 inside the heart [START_REF] Tenforde | Magnetically induced electric elds and currents in the circulatory system[END_REF]. However, the induced current signicantly decreases with the distance from a small region of the aortic arch. In [START_REF] Kinouchi | Theoretical analysis of magnetic eld interactions with aortic blood ow[END_REF][START_REF] Tenforde | Magnetically induced electric elds and currents in the circulatory system[END_REF], the authors predicted a maximum current of 1.15A m -2 at 5T , 2.20A m -2 at 10T and 3.20A m -2 at 15T . These values are two to three times smaller than ours. Two factors can explain these dierences. Firstly, as recalled in Section 4.1.3, the orientation of B in [START_REF] Kinouchi | Theoretical analysis of magnetic eld interactions with aortic blood ow[END_REF] is not along the main body axis. Secondly, the ow pattern obtained in our realistic aortic arch solving the nonlinear Navier-Stokes equations is of course very far from a steady state ow in an innite cylinder. The induced electrical current has therefore a much more complicated behavior than in the simplied models.

Conclusions

A model made up of three components has been proposed to simulate the eects of MHD in the aorta and its impact on ECGs: electrophysiology in the heart, MHD in the aorta and electrostatics in the rest of the body. The main contributions of this article are 1. the use and the justication of the inductionless approximation of the MHD, which is more convenient in the blood ow regime than conventional MHD equations, 2. an ecient numerical scheme to solve the MHD equations in this regime, 3. the denition of several benchmarks to verify the numerical procedures, and 4. the full coupling of the three compartments that allowed us to study whether the current induced in the aorta provides a feedback on the heart. Standard 12-lead ECGs have been computed with and without a magnetic eld. A larger Twave was obtained in the presence of a magnetic eld, as expected. The electric potential and the electrical current in the aorta have been computed. According to these computations, we found that the ECG is approximately the superposition of the electrophysiology and of the MHD eects (see [START_REF] Gold | Magnetohydrodynamic pipe ow. part 1[END_REF]). The impact of the aortic MHD eect on the electric activity of the heart has been evaluated and seems limited, even for strong magnetic elds.

Despite the complexity of the model, it still could be enriched. In addition to the usual questions of blood ow modelling (inlet ow, outlet boundary conditions, the length of the artery to be kept, etc.), several specic questions need to be discussed. For example, the blood ow in the pulmonary artery and in the heart might be also considered. We neglected these ows because we focused on the ECG, but they could have an eect on the cardiac electrophysiology. The elasticity of the aorta could be considered as in [START_REF] Moireau | External tissue support and uid-structure simulation in blood ows[END_REF]. A major, but computationally expensive, improvement could be to link directly the aortic ow to an electromechanical model of the myocardium [START_REF] Chapelle | Numerical[END_REF], and in this way the aortic inow pulse would not be a given parameter of the model.

The data for the conductivities of dierent tissues used in the literature [START_REF] Nijm | Inhomogeneous human torso model of magnetohydrodynamic blood ow potentials generated in the mr environment[END_REF][START_REF] Kinouchi | Theoretical analysis of magnetic eld interactions with aortic blood ow[END_REF] can vary by an order of magnitude. It could be interesting to perform a sensitivity analysis in order to evaluate the impact of these conductivities on the MHD eect. In the same direction, here we only took into account a few dierent tissue conductivities in the torso: it would also be interesting to explore a more precise description of the tissues (fat, muscles, more bones, etc.).

The aim of the present study was to demonstrate that the computational tool was able to capture a phenomenon that has actually been observed in MRI machines. In the future, the model could be used, for example to optimize the electrode locations so as to be less sensitive to the MHD eect, or on the contrary to extract more information from the MHD eect about the ow or the electrical properties of the tissues.

Parameters

Value units T time period 700 ms σ i,l intern longitudinal conductivity 3.0 10 -3 S cm -1 σi,t intern transverse conductivity 3.0 10 -4 S cm -1 σ e,l extern longitudinal conductivity 3.0 10 -3 S cm -1 σe,t extern transverse conductivity

4. 1 . 1 .Figure 2 .

 112 Figure 2. Left : the conguration. Middle : velocity (parabola) and potential as a function of y; analytical (solid line) and numerical (crosses) solutions. Right : inlet ux as a function of time, comparison between stationary analytical solution and numerical solutions.

( 10 ) 3 - 2 y 2 + 2 y.

 103222 Solving[START_REF] George | Fully automatic mesh generator for 3d domains of any shape[END_REF], we obtain u = Cy(L y -y), and φ = CB0 3 y CB0Ly G-2Cη σB0 y + φ 0 , where C and φ 0 are xed by the boundary conditions. Two cases are considered. Imposed potentials Assuming that φ(0) = φ 0 and φ(L y ) = φ 0 + δφ, with given constants φ 0 and δφ, we obtain C = 6 η G-σB0 δφ/Ly 12+(σ/η)B 2 0 L Fig. 2 shows the very good agreement of the numerical result and this analytical solution. The parameters are B 0 = -2, σ = 0.5, η = 0.05, φ 0 = 2, δφ = 3, G = 2, L y = 1, H y = 7, τ = 10 -2 (SI units, arbitrary values).

Figure 3 .

 3 Figure 3. Left : the conguration. Middle : velocity (parabola) and potential as a function of y; analytical (solid line) and numerical (crosses) solutions. Right : inlet ux as a function of time, stationary analytical solution and numerical solutions.

( 8 )

 8 at y = L y , we obtain C = 6 η G-σB0 δφ/Hy 12+(σ/η)B 2 0 L 3 y /Hy and φ

4. 1 . 3 .

 13 Comparison with a numerical solution from the literature In[START_REF] Kinouchi | Theoretical analysis of magnetic eld interactions with aortic blood ow[END_REF], a 2D model was proposed to study the electric potential induced by the MHD eects in the aorta. This test case cannot be used for quantitative estimations of the MHD eect in MRI because the magnetic eld B is imposed along the backfront axis and not along the feethead axis. Nevertheless, for validation

Figure 4 .

 4 Figure 4. Gold validation: velocity contours. Left: B = 0 T , colorbar scale [0, 119.9]cm s -1 . Center: B = 10 T , scale [0, 116.3]cm s -1 . Right: B = 100 T , scale [0, 42.6]cm s -1 .

Figure 5 .Figure 6 .

 56 Figure 5. Gold validation: potential and electrical current on a cut plane. Left: potential, B = 10 T , scale [-58.6, 58.6]mV . Middle left: potential, B = 100 T , scale [-274, 274]mV . Middle right: current, B = 10 T , scale [0, 2.97]A m -2 . Right: current, B = 100 T , scale [0, 16.61]A m -2 .

hmax 4 i=1

 4 area(faces) volume (≥ 1). We give the percentage of tetrahedra having a good quality (Qm ≤ 10).

Figure 7 .

 7 Figure 7. Meshes of the three domains: aorta, heart and body with 3 dierent tissues (bones, lungs, other tissues), body surface with some points used to compute the ECG.

Figure 8 .

 8 Figure 8. Left: lead I of the ECG with B = 0, in mV (solid line) and imposed inlet ow Q (dashed-line) in cm 3 s -1 , as a function of time (ms). The total ux is approximately 5L min -1 . Middle: maximum potential value for the lead I (mV ) as a function of B (T ). The peak value in QRS complex is shown (horizontal broken line, constant for all B). Right: lead I ( D1) of the ECG with B = 10T , in mV , for 3 mesh renements. The 3 curves are almost superimposed.

Fig. 9 Figure 9 .

 99 Fig.9shows that the ECG is roughly the superposition of the ECG without a magnetic eld

I(D 1 )

 1 II(D 2 ) III(D 3 ) aVR(D 4 ) aVL(D 5 ) aVF (D 6 ) α i 1 10 -3 0.35 10 -3 -0.65 10 -3 -0.67 10 -3 0.82 10 -3 -0.15 10 -3

Figure 10 .

 10 Figure 10. Potential (mV ) at two points on the heart surface (where the MHD eect is the strongest) as a function of time (ms), for B = 0, 3, 8, 20T .

Figure 11 .

 11 Figure 11. Leads I (Left) and II (Right) with MHD eects (B = 10T ). Comparison between the strong coupling strategy (solid line) and the weak coupling strategy (dashed line).

Figure 12 .

 12 Figure 12. Electric potential. From top to bottom: B = 0T (no MHD), B = 1, 3, 8T . From left to right, time instants: t = 15ms, 20ms, 30ms, 100ms, 200ms, 250ms, 300ms. Scale: [-2, 2]mV .

Figure 13 .

 13 Figure 13. Electric potential on the aorta at t = 210ms (systole): for B = 0T (left, scale: [-0.5, 0]mV ), B = 3T (middle, scale: [-25, 25]mV ) and B = 8T (right, scale: [-70, 70]mV ).

Figure 14 .

 14 Figure 14. Electrical current intensity on the aorta at t = 210ms (systole). Left: B = 3T (scale: [0, 1.84]A m -2 ). Right: B = 8T (scale: [0, 4.9]A m -2 ).

  , and π i is the solution to the ordinary dierential

	equation C i	dπi dt + πi R d,i = Q i . The parameters R d,i , R p,i and C i represent the distal and proximal
	resistances and the capacitance C i respectively (values taken from [21]).

Table 1 .

 1 Gold validation: analytical and numerical results on various mesh renements. The

				188.496	120	0	0
	inductionless, M 1	0	0	187.08	119.37	0	0
	inductionless, M 2	0	0	188.02	119.79	0	0
	inductionless, M 3	0	0	188.29	119.91	0	0
	exact MHD	10	1	184.706	116.400	117.389	2.9546
	inductionless, M 1	10	1	183.30	115.79	117.01	2.973
	inductionless, M 2	10	1	184.24	116.20	117.30	2.970
	inductionless, M 3	10	1	184.50	116.31	117.36	2.965
	exact MHD	20	2	174.758	107.134	221.250	5.6678
	inductionless, M 1	20	2	173.38	106.56	220.36	5.698
	inductionless, M 2	20	2	174.30	106.94	220.94	5.695
	inductionless, M 3	20	2	174.56	107.05	221.06	5.688
	exact MHD	100	10	90.634	42.882	550.171	16.6321
	inductionless, M 1	100	10	89.02	42.21	542.99	16.45
	inductionless, M 2	100	10	90.06	42.65	547.75	16.61
	inductionless, M 3	100	10	90.39	42.78	549.22	16.65

dimensionless Hartmann number is dened as Ha = BR σ/η (= 0.1B here).

Table 2 .

 2 Meshes

for the heart, the aorta and the rest of the body. Mesh quality of a tetrahedron dened as: Qm = 1 6 √ 6

Table 3 .

 3 Potential and current as a function of B: maximum and minimum potentials in the aorta (φ A ) and on the skin (φ S ), peak value for leads I, II and III, peak current modulus in the aorta (j A ). Values at t ≈ 210ms, peak value for the MHD eects (T -wave), and value during the R-wave at t = 27ms (it is identical for all B). The Hartmann number is computed as Ha = B R A σ bl /η, see TableA2.

						6.8	-6.3	2.0	8.3	1.83	1.47	0.66	0.23
		0	0	-0.5	-0.02 0.48 -0.06	0.5	0.56	0.81	0.42	-0.40	0.039
		1	0.12	-9.7	8.5	18.2	-1.0	1.23	2.23	1.05	0.49	-0.57	0.61
		2	0.24 -19.5	17	36.5	-2.0	2.35	4.35	1.30	0.56	-0.76	1.22
	T -	3	0.36	-29	25.6	54.6	-3.0	3.33	6.33	1.59	0.64	-1.02	1.84
	wave	5	0.59	-49	43	92	-4.9	5.7	10.6	2.49	0.87	-1.61	3.06
		8	0.95	-78	68	146	-8.1	8.85	16.9	3.85	1.35	-2.50	4.90
		10	1.19	-97	86	183	-9.8	11	20.8	4.76	1.67	-3.10	6.13
		15	1.78 -146	128	274	-15	16.6	31.6	7.05	2.47	-4.58	9.19
		20	2.38 -195	171	366	-20	22	42	9.34	3.27	-6.08	12.2

Table A1 .

 A1 Main parameter values in the heart. See[START_REF] Boulakia | Mathematical modeling of electrocardiograms: A numerical study[END_REF] for more details. cm -1 s -1 4 10 -3 kg m -1 s -1 σ bl blood conductivity 6.25 10 -4 mA 2 s 3 g -1 cm-3

						1.2 10 -3	S cm -1
		Am Cm φH	membrane area per vol membrane capacitance per area electric potential	200 10 -3 ≈ 10	cm -1 mF mV
	T R A ρ η	Parameters time period aortic arch radius blood density blood viscosity	Value 0.7 0.95 1.06 0.04	{cm/g/s/mA} units Value (SI unit) s cm 9.5 10 -3 m g cm -3 1.06 10 3 kg m -3
							0.625S m	-1
	B umax qmax S Q	magnetic eld mean inlet velocity 98 1 to 20 460 max blood ow 4.66 inlet area blood ow rate ≈ 5.36	g s -2 mA -1 cm s -1 cm 3 s -1 cm 2 L min -1	1 to 20 T 0.98m s -1 0.46 10 -3 m 3 s -1 4.66 10 -4 m 2
	φA	electric potential	≈ 100	g cm 2 s	-3 mA -1	≈ 10mV

g

Table A2 .

 A2 Main parameter values in the aorta expressed in {cm, g, s, mA} (and recalled in SI units). umax denotes the maximum mean blood velocity at the inlet (during systole).

		Parameters	Value	units
	σT σ bl σ bone σ l φT	S cm -1 blood conductivity 6.25 10 -3 S cm -1 tissue conductivity 6. 10 -4 5.65 10 -5 S cm -1 bone conductivity lung conductivity 5.10 -4 S cm -1 electric potential ≈ 10 mV

Table A3 .

 A3 Main parameter values in the body. simulation of the electromechanical activity of the heart. In Functional Imaging and Modeling of the Heart, number 5528 in Lecture Notes in Computer Science, pages 357365. Springer-Verlag, 2009.

Acknowledgment: We wish to thank C.A. Taylor and C.A. Figueroa for the mesh of the aorta.

Appendix A. Data values

In this paper, the MHD computations in the aorta are performed using the {cm, g, s, mA} system, that has the advantage of giving a convenient scaling for the blood ow, while preserving the magnetic eld scaling (1T ←→ 1{cm/g/s/mA}). In the results, the electric potential is given in mV and the current in A m -2 . The computations of Section 4.2 use the following data. See Table A1 for the heart, Table A2 for the aorta, and Table A3 for the body.