X Artru
email: x.artru@ipnl.in2p3.fr

C Ray
email: c.ray@ipnl.in2p3.fr

PARTICLE ACCELERATION IN A HELICAL WAVE GUIDE

The electromagnetic wave field propagating in a helical wave guide is decomposed in an angular momentum basis. Eigenmodes are calculated using a truncation in l and a discretisation of the boundary condition. Modes slightly slower than light are considered for relativistic particle acceleration.

Introduction

An electromagnetic wave can accelerate charged particles over long enough distance if 1) it has a longitudinal electric component 2) its phase velocity is equal or close to the particle velocity, therefore less than c. A free plane wave misses the first condition and a wave guided by a metallic straight tube miss the second condition. The usual solution to satisfy both conditions is a chain of coupled resonating cavities. Another solution is the helically bent tube. Here the "photon path" is longer than the corresponding path in a straight tube, so that some modes can be made slower than light, while having a non-zero longitudinal electric field on the helicity axis.

A helical wave guide is defined by its axis Oz, its period Λ = 2π/q and its base curve r(s) = (X(s), Y (s)) which is the intersection of the wall with the z = 0 plane. The wall surface is generated by helices H s of parametric equation (x, y, z) = (R(qz)r(s), z)

where the 2 × 2 matrix R(qz) describes the rotation of angle qz about 0z. The traveling wave tube, used for generation of radar waves, is a special kind of helical wave guide. It is made of a cylindrical metallic tube inside which a thick conducting helical wire is fixed. Its base is not simply connected : it is a disk with a hole corresponding to the section of the wire at z = 0. Therefore the base curve has two disconnected pieces : the circular base of the cylinder and, inside it, a small closed curve. The same device, with a superconducting helix, has been studied for particle acceleration [START_REF] Klein | Proceedings of the 1972 proton linear accelerator conference[END_REF][START_REF] Aron | Stanford[END_REF].

A tube where the wall itself possesses a helical structure can be made by drawing helical grooves on the internal face of an initially cylindral tube. In the case where the groove section is made of rectangular steps, the propagation modes can be calculated using an approximate field-matching technique [START_REF] Wang | [END_REF], or a 3-dimensional computer code [START_REF] Wallet | Technical characteristics of a novel helical-groove traveling-wave tube structure[END_REF]. In this paper we will consider helical tubes which have smoothly curved walls, associated to smooth base curves. To calculate the modes we will use a truncated partial wave expansion and enforce the boundary conditions only at discrete points of the base curve, called sutture points, the number of sutture points being equal to the number of partial waves. For definiteness, we will take the base cuve to be a ex-centred circle, as shown in Fig. 1. The guide parametres are then

• a = radius of the base circle,

• ε ∈ [0, 1] = relative amplitude of the helical deformation. The centre of the base circle is at distance εa from Oz.

• q = 2π/(helix period).

The average angle of an helix with Oz (i.e., π 2 minus the pitch angle) is ᾱ = arctan(qa). Another significant quantity is the maximum helix angle of the wall with Oz, β = arctan(εqa). The case β = 0 (q = 0 or ε = 0) corresponds to a straight cylindrical guide.

The scalar field case

We first study the simpler case of a scalar field Φ(t, x, y, z), which obeys the Klein-Gordon equation (in units c = 1)

(-∂ 2 t + ∂ 2 x + ∂ 2 y + ∂ 2 z ) Φ(t, x, y, z) = 0 , (2)
with the boundary condition Φ(t, x, y, z) = 0 on the wall .

In a fixed z transverse plane we introduce a "rotating frame" whose basis vectors rotate about Oz by the angle qz when z increases. The fixed-and rotating-frame coordinates, r = (x, y) and R = (X, Y ), are related by x + iy = e iqz (X + iY ).

The corresponding polar coordinates are (r, ϕ) and (r, ϕ -qz) respectively. We look for a solution of the form

Φ(t, x, y, z) = e iP z-iωt Ψ(X, Y ) . (4)
P is the longitudinal pseudo-momentum (in the Bloch wave sense). It is the eigenvalue of the operator

P helic = p z + qJ z = -i∂/∂ z -iq∂/∂ ϕ (5)
which is conserved due to the invariance under the helical displacement {translation by ∆z} × {rotation by q∆z} .

For a given ω the Klein-Gordon equation with the boundary condition has a finite number of solutions or modes. P and ω are linked by the dispersion relation P = P m (ω) or ω = ω m (P ). The phase and group velocities are v ph = ω/P , v gr = dω/dP . Since Φ is an eigenvector of (5) we can decompose it in free cylindrical partial waves of angular momenum l and longitudinal momentum p l = P -lq

Φ(t, x, y, z) = e -iωt +∞ l=-∞ a l e ilϕ e i(P -lq)z ψ l,k(l) (r) , (7)
where

k(l) or k l = (ω 2 -p 2 l ) 1/2 = ω 2 -(P -lq) 2 1/2 (8)
is the transverse momentum of the l component and ψ l,k (r) is a solution of the radial equation

ψ ′′ + ψ ′ /r + (k 2 -l 2 /r 2 )ψ = 0 . (9) k l is real for (P -ω)/q ≤ l ≤ (P + ω)/q , ( 10
)
otherwise it is pure imaginary : k l = iκ l , with κ l = (P -lq) 2 -ω 2 1/2 . The solution of (9) is proportional to the Bessel function J l (kr) or the modified Bessel function I l (κr) ≡ J l (iκr)/i l . We will fix the normalisation of ψ l,k (r) by

ψ l,k (r) = J |l| (kr)/k |l| = I |l| (κr)/κ |l| . (11)
With k = iκ the last two expressions are equivalent. The k |l| denominator allows a finite limit for ψ l,k (r) when k → 0. Let us write the boundary condition for t = z = 0. Il will be satisfied also at all t due to the invariance under time translation and all z due to the invariance under helical displacement. For every point r(s) = (X(s), Y (s)) of the base curve one has the constraint

a l e ilϕ(s) J |l| [kr(s)]/k |l| = 0 , (12)
where (r(s), ϕ(s)) are the polar coordinates of r(s).

The truncation-and-sutture method

One keeps only a finite set L of partial waves and selects an equal number of "sutture points" r 1 , ...r N of the base curve (see Fig. 1). For instance L = [l inf , l sup ] with l sup = l inf + N -1. The a l are approximatively determined by the system of linear equations

l∈L a l e ilϕn J |l| (kr n )/k |l| = 0 , 1 ≤ n ≤ N , (13)
(r n , ϕ n ) being the polar coodinates of r n . A nonzero solution exists if the determinant est zero:

det{C nl } = 0 , where C nl = e ilϕn J |l| (kr n )/k |l| . ( 14
)
This is the dispersion relation in the truncation-and-sutture approximation. The following conditions may be imposed on L and the r n 's :

• the r n 's should give a precise enough definition of the base curve,

• l (ϕ n+1 -ϕ n ) should not be too large compared to unity,

• k l (r n+1 -r n ) should not be too large compared to unity.

The ultimate criterium of the quality of the approximation is the stability of the roots of ( 14) when N is changed by one or two units. The coefficients a l should be stable also, at least for l far enough from the border of L.

Case of a symmetric base. If the base is symmetrical about Oz, the modes are either even or odd under this symmetry. One can study the even and odd modes separately by taking either even or odd l's in L. Similarly, if the base has a ternary symmetry (invariance under rotation by 2π/3), L can be restricted to either l = 0 or l = 1 or l = 2 modulo 3. To avoid redundant equations in (14), the points r n should be all different modulo the symmetry. One can choose them on one part of the base curve only, the other parts being deduced by the symmetry.

Numerical example in the scalar case

We took a circular base with parameters ε = 0.4, q = 1 and imposed the phase velocity v ph = ω/P = 1 -10 -5 . Plotting ln |det{C nl }| as a function of ωa for different sets L and equidistant r n 's, with N ranging from 5 to 21, we found zeros at ωa = 5.0, 7.0 which are stable at the 0.1 precision when L contains at least the set [-1, +5]. Another root, ωa = 8.5, seems to stabilise for L ⊃ [-1, +7]. Higher modes are seen, but would require more partial waves to get stable. Fig. 2 shows the scalar field Ψ(X, Y ) = Φ(0, X, Y, 0) of the ωa = 5.0 mode. Due to the low number of sutture points, the boundary condition Ψ(X, Y ) = 0 at the boundary is only appoximately satisfied. One observes that Ψ(X, -Y ) = Ψ * (X, Y ). This property comes from the symmetry of the base about the X axis. The maximum of |ℜΨ| is not at the centre of the base, but farther from the helix axis. This is a manifestation of the centrifugal effect. One can also observe that the phase arg Ψ = arctan(ℑΨ/ℜΨ) increases with the azimuth, which means that the average angular momentum is positive, as expected.

The Maxwell field

We replace the Klein-Gordon equation by the Maxwell equations in vacuum,

∇ × E = -∂ t B ∇ • B = 0 ∇ × B = ∂ t E ∇ • E = 0 . ( 15
)
These equations are invariant under the duality B → E, E → -B. In the following we gather the fields E and B in a 6-component vector F = {E, B}. In analogy with (4) we look for solutions of the form

F(t, x, y, z) = e iP z-iωt R(qz)F (X, Y ) , (16)
in the complex notation (the physical field is the real part of it). F = {E, B} is the field in the z = 0 plane and R(qz) is the rotation matrix of (1) applied to the vectors E and B :

E x = E X cos qz -E Y sin qz , E y = E X sin qz + E Y cos qz (idem for B). ( 17
)
To write the boundary conditions we introduce the vectors

σ(s) = r ′ (s) =   X ′ (s) Y ′ (s) 0   , τ (s) = ẑ + q ẑ × r(s) =   -qY (s) qX(s) 1   , (18)
which are tangent to the wall at z = 0, and the normal vector

ν(s) = σ(s) × τ (s) =   Y ′ (s) -X ′ (s) q(XX ′ + Y Y ′ )   . ( 19
)
The boundary conditions are E • σ = 0, E • τ = 0 and B • ν = 0. In cylindrical coordinates, they read

E r σ r + E ϕ σ ϕ = 0 , ( 20
) qr E ϕ + E z = 0 , ( 21
) B r σ ϕ + (qrB z -B ϕ )σ r = 0 . ( 22
)
These 3 conditions are not independent. Only two suffice. They break the duality B → E, E → -B.

F can be decomposed in a manner similar to (7). However, since the photon has two polarisation states, there are two free partial waves for each l, which is now the total (spin + orbital) angular momentum.

Study of a free partial wave

Let us first study separately one partial wave of definite ω, l and p (here p = p l = P -lq). In cylindrical coordinates it takes the form

F ω,p,l (t, r, ϕ, z) =        E r E ϕ E z B r B ϕ B z        = e -iωt e ipz e ilϕ        E r (r) E ϕ (r) E z (r) B r (r) B ϕ (r) B z (r)        . (23)
In this subsection the quantities written in caligraphic depend only on r. The Maxwell equations become

(l/r) E z -p E ϕ = ω B r , -∂ r E z + ip E r = iω B ϕ , r -1 + ∂ r E ϕ -i(l/r) E r = iω B z , (24)
and the dual equations obtained by B → E, E → -B. These equations are redundant and one can omit the third one and its dual.

The TM/TE basis. A first basis is made of the so-called transverse magnetic (TM) and transverse electric (TE) states. In cylindrical coordinates,

F (TM) ω,p,l (r) =        E r (r) E ϕ (r) E z (r) B r (r) B ϕ (r) B z (r)        (TM) = k -|l|-2δ l,0         +ip k J ′ |l| (kr) -p l/r J |l| (kr) +k 2 J |l| (kr) +ω l/r J |l| (kr) +iω k J ′ |l| (kr) 0         , (25)
with k = (ω 2 -p 2 ) 1/2 . The TM and TE fields are dual, up to a factor i : B (TE) = -iE (TM) , E (TE) = iB (TM) . For the l = 0 waves,

F (TM) ω,p,0 (r) =        -ip ψ 1 0 ψ 0 0 -iω ψ 1 0        , F (TE) ω,p,0 (r) =        0 +ω ψ 1 0 -p ψ 1 0 -iψ 0        , (26)
ψ being defined in (11).

The chiral basis. When k → 0 with l = 0, F (TM) and F (TE) become equal up to a sign and they no more form a complete basis. We avoid this degeneracy with the combinations F + = F (TE) + F (TM) , F -= F (TE) -F (TM) , which are respectively the states of right and left chirality, or helicity χ = +1 and χ = -1. They are such that E = iχB. For real k a right-handed (resp. lefthanded) solution is a superposition of plane waves of positive (resp. negative) helicity. More precisely, we normalise these states as :

F ω,p,l,χ = -2 F (TE) + χF (TM) χp + s l ω =         -i ψ λ-1 + (p -χ ′ ω) 2 ψ λ+1 s l ψ λ-1 -(p -χ ′ ω) 2 ψ λ+1 2(p -χ ′ ω) ψ λ -χ ψ λ-1 + (p -χ ′ ω) 2 ψ λ+1 -iχ ′ ψ λ-1 -(p -χ ′ ω) 2 ψ λ+1 -2iχ (p -χ ′ ω) ψ λ         . (27)
Here s l = sign(l), χ ′ = s l χ, λ = |l| and ψ λ = J λ (kr)/k λ . Use was made of the recurrence relations

(l/x) J l (x) = [J l-1 (x) + J l+1 (x)]/2 , J ′ l (x) = [J l-1 (x) -J l+1 (x)]/2 .
The denominator (χp + s l ω) eliminates a kinematical zero, so that these states have finite and linearly independent limits when p → -χ ′ ω. We will take the chiral basis except for the l = 0 states, where we keep the TM-TE basis (26).

The partial wave expansion

According to (26-27), the decomposition ( 7) is replaced by

F(t, r, ϕ, z) = e -iωt e iP z a 0 F (TM) ω,P,0 (r) + b 0 F (TE) ω,P,0 (r) + e -iωt l =0 e ilϕ e ip(l)z a l F ω,p(l),l,χ ′ =+1 (r) + b l F ω,p(l),l,χ ′ =-1 (r) , (28)
with p(l) = P -lq. For l = 0 the a l (resp. b l ) coefficients multiply the states where the chirality has the same sign as l (resp. opposite sign to l).

The boundary conditions (20-21) read

0 = (pσ r a 0 + iωσ ϕ b 0 ) ψ 1 + l =0 e ilϕ { a l (σ r + is l σ ϕ ) ψ λ-1 + (p -ω) 2 (σ r -is l σ ϕ ) ψ λ+1 + b l (σ r + is l σ ϕ ) ψ λ-1 + (p + ω) 2 (σ r -is l σ ϕ ) ψ λ+1 } , (29) 0 = a 0 ψ 0 + ωqr ψ 1 b 0 + l =0 e ilϕ { a l s l qr (ψ λ-1 -(p -ω) 2 ψ λ+1 ) + 2(p -ω) ψ λ + b l s l qr (ψ λ-1 -(p + ω) 2 ψ λ+1 ) + 2(p + ω) ψ λ } . ( 30
)
In these equations, λ = |l|, p = p l = P -lq, k = k l = ω 2 -p 2 l 1/2 , the argument of ψ λ , ψ λ-1 , and ψ λ+1 is k l r(s), ϕ = ϕ(s) and σ = σ(s) is given by (18). Generalising the truncation-and-sutture method of Section 1.2 we write conditions (29-30) for N points r n of polar coordinates (r n , ϕ n ), keeping only N values of l. It gives 2N linear equations linking N coefficients a l and N coefficients b l . The dispersion relation is the vanishing of the 2N × 2N determinant of this system, the coefficients of which are given in Appendix (35-36). L and r 1 , ...r N are chosen like in the scalar case.

We have tested this method with a straight (q = 0) cylindrical wave guide. The well-known TM and TE lower modes were found using no more than ten points, even when the Oz axis was not the guide axis, but shifted by about half a radius.

Application to the acceleration of relativistic particles

A particle moving along the axis of the helical wave guide is accelerated or decelerated by the E z component at r = 0, whose physical value is ℜ{a 0 e iP z-iωt }.

Besides, the l = ±1 partial waves give a transverse field at r = 0, the components of which are, in the rotating frame

E X = ℑ{(a 1 + b 1 + a -1 + b -1 ) e iP z-iωt } E Y = ℜ{(a 1 + b 1 -a -1 -b -1 ) e iP z-iωt } B X = ℜ{(-a 1 + b 1 + a -1 -b -1 ) e iP z-iωt } B Y = ℑ{(a 1 -b 1 + a -1 -b -1 ) e iP z-iωt } . (31)
If an electron moves along the axis at the same velocity as the wave (z = v(t-t 0 ) with v = v ph ), it will experience a Lorentz force f = (-e) (E + v × B) which is constant in the rotating frame. In the ultrarelativistic case

(v ph ≃ 1), this force is f X f Y ≃ (-e) × 2ℜ e -iωt0 b 1 -i 1 -b -1 i 1 . (32)
In the laboratory frame, it is the force excerted by a helical undulator having the same period as the helix. We therefore get a parasitic undulator radiation. It can be avoided with a base which has a symmetry of order S ≥ 2 about Oz, choosing a mode involving only the waves l = 0 modulo S.

Numerical results

Like in the scalar case (Section 1.3) we took a circular base with parametres ε = 0.4, q = 1 and imposed the phase velocity v ph = 1 -10 -3 . Plotting ln |det{C nl }| or, better, arg |det{C nl }| we located modes at ωa = 2.28, 3.6, 4.8 and 5.0 which are stable at least at the 0.1 precision for N ∼ 15 (they were the same for L = [-6, +8], [-5, +9] and [-3, +11]). The first three modes seem to be mostly of the TE type : |E z (0)/B z (0)| = 0.3, 0.25 and 0.4. We give below more details about the fourth one (ωa = 5.0) which is mainly TM : |E z (0)/B z (0)| = 5.4. Fig. 3 shows the cross section of the fields along the segment x ∈ [-0.6, +1.4], y = 0, which is a diametre of the base circle (in units a = 1). At the extremities the boundary condition B r = 0, which was not directly imposed by (29-30) is well satisfied. A figure of merit of the helical wave tube as an accelerating device is the ratio of the accelerating field over the strongest electric field at the boundary,

f.o.m. = E z (0)/ max(|E(r(s))|) . (33)
Indeed a too strong field at the boundary may trigger an electron emission. Fig. 4 indicates that the strongest boundary field is at r = r max and Fig. 3 gives f.o.m. ≃ 1/8. Another relevant quantity is the group velocity v g = dω/dP , which can be easily calculated from (14) when the phase velocity v ph = ω/P is not imposed. For the ωa = 5.0 mode studied above, v g = 0.68. The guide can therefore be quickly filled with the wave and quickly emptied. To save energy, one may use a recirculated wave.

We also have preliminary results for a waveguide with a stronger helical bending: ε = 0.6, q = 1.5. Imposing v ph = 1 -10 -3 , three first modes are seen at ωa = 2.2, 3.3 and 4.0 with |E z (0)/B z (0)| = 0.35, 0.3 and 0.8 respectively. The number of sutture points, N ∼ 10 -20, was enough to obtain stable roots but not yet stable field amplitudes. This is like in the variational method in quantum mechanics, where the energy converges faster than the wave function.

Conclusion

This study has shown that a linear accelerator for ultrarelativistic particles may work with a helical wave guide. Until now we have made detailed numerical calculations only for a guide with circular base, one set of parametres {ε, q} and the lower TM-like mode. The figure of merit (accelerating field ) / (strongest electric field at the boundary) of this mode was only about 1/8, but we hope that different guide parametres or different base shapes can give higher values.

The study has also validated the truncation-and-sutture method. This economical method seems to be well suited to calculate the modes of a scalar or vector field in a 2-or 3-dimensional homogeneous cavity of arbitrary, but smooth enough shape.

Figure 1 :

 1 Figure 1: Base curve (b.c.) and sections of the helical wave guide at successive values of z. The rotating frame {OX, OY } of the last section is shown. In this example the base curve is an ex-centered circle of radius a. The balls (s.p.) represent the sutture points.

Figure 2 : 3 -

 23 Figure 2: 3-dimensional representation of the scalar field Ψ(x, y) = Φ(0, x, y, 0) of the ωa = 5.0 mode, for a circular base with ε = 0.4 and q = 1. Left: real part; right: imaginary part.

Fig. 4

 4 shows the 3-dimensional plot of |E r (x, y)|, |E ϕ (x, y)|, |E z (x, y)| and |B z (x, y)|. We did not plot |B r | and |B ϕ |, which look like |E ϕ | and |E r |.

Figure 3 :Figure 4 : 3 -

 343 Figure 3: Cross section of the electromagnetic field F (x, y) = F(0, x, y, 0) along the line y = 0 for the TM-like mode at ωa = 5.0 (circular base, ε = 0.4 and q = 1). Left curves, from top to bottom : |E r (x, y)|, |E ϕ (x, y)| and |E z (x, y)| at y = 0. Right curves : same for B r , B φ and B z .

6

  APPENDIX : The M matrix Defining C n,l = coefficient of a l in the condition E • σ = 0 at r n , C n,l+N = coefficient of b l ≡ a l+N in the condition E • σ = 0 au point r n , C n+N,l = coefficient of a l in the condition E • τ = 0 at r n , C n+N,l+N = coefficient of b l in the condition E • τ = 0 at r n , the 2N × 2N matrixof the linear system of equations expressing the boundary conditions is decomposed in four N × N submatrices : C n+N,l ... | ... C n+N,l+N ... are, for l = 0,C n,l = ψ λ-1 σs l + (p -ω) 2 ψ λ+1 σ-s l e ilϕn C n,l+N = ψ λ-1 σs l + (p + ω) 2 ψ λ+1 σ-s l e ilϕn C n+N,l = s l qr (ψ λ-1 -(p -ω) 2 ψ λ+1 ) + 2(p -ω) ψ λ e ilϕn C n+N,l+N = s l qr (ψ λ-1 -(p + ω) 2 ψ λ+1 ) + 2(p + ω) ψ λ e ilϕn ,(35)and for l = 0,C n,0 = P ψ 1 σ r C n,N = iω ψ 1 σ ϕ C n+N,0 = ψ 0 C n+N,N = ωqr ψ 1 .(36)In these formulas p = p l = P -lq, k = k l = ω 2 -(P -lq) 2 1/2 , λ = |l|, s l = sign(l), ψ λ = J λ (k l r n )/k λ l and we have introduced the unitary complex number σn = σ r (r n ) + iσ ϕ (r n ) | σ(r n )| = e -iϕn σ x (r n ) + iσ y (r n ) | σ(r n )| .


Antoine Bulteau
email: antoinebulteau@g2elab.grenoble-inp.fr

Yves Lembeye

David Frey

Design and limit of a 4 legs inverter with unbalanced grid injection operation connected to a Quad Active Bridge converter

come

Design and limit of a 4 legs inverter with unbalanced grid injection operation connected to a Quad Active

Bridge converter

Introduction

Currently, a converter network architecture is proposed based on a multiport converter (MC). This MC is composed of two converters showed in Fig 1 , the Quad Active Bridge, connecting the DC components and the 4 legs 2 levels inverter, allowing the interconnection to the AC grid making it possible to exchange energy between these equipments.

In this article, the focus will be done on the 4 legs inverter part which is represented in Fig 2 . First, we will discuss the constraints and limitations of connecting this converter to a grid and exchange power in an unbalanced way. Next, we will propose an adapted control and a more efficient design for this converter.

The study case we propose to address here is a three-phase network with single-phase domestic loads. This network is therefore potentially * Institute of Engineering Univ. Grenoble Alpes unbalanced as shown in Fig 2 . In this application, the aim of the inverter is to furnish ancillary services to the network by acting as an energy router. Therefore, it helps to increase the load balancing by transferring power from one phase to another one. As we will see, this functioning will add new constraints in the inverter design. This configuration brings constraints that will be explained in section 2. Section 3 will validate these constraints by simulation. Section 4 will present a validation of the simulation by experimentation. Finally, we will conclude in section 5. The design hypothesis is that the neutral current limit is the same as the one of the phases. The neutral current can be expressed by one equation with 6 unknowns (3 active Pi and reactive Qi powers). This equation ( 1) is used in an algorithm to determine the cases that can fulfill the requirements for the inverter. To obtain these values, it was chosen to vary the active and reactive power independently for each phase with a linear increment from -Pmax to +Pmax and identically for the reactive power per phase from -Qmax to +Qmax. For each calculation, a value of the neutral current is obtained Only 49% of the cases have its neutral current acceptable considering this nominal current limitation. If the criterion is to accept all configurations, then the neutral part needs an oversizing and specific protection to accept these cases. This possibility is not in line with this converter which induces an extra cost.

Inverter Design

Neutral current limitation

Attended those multiple converters have to be implemented simultaneously, this solution should not be economically sustainable. This constraint can be deflected in using centralized control [START_REF] Han | Review of Active and Reactive Power Sharing Strategies in Hierarchical Controlled Microgrids[END_REF] which enables to cover any power request by using several converters. The current has a counterclockwise, whereas the fluctuated power is in opposite rotation. In addition, both systems are not the same frequency (50Hz for the current frame, 100Hz for the fluctuating power). Therefore, it indicates that there is no direct link between the amount of neutral current and the fluctuating power in the DC bus.

Table -1 gives the result obtained for one configuration having 3680VA per phase corresponding to a current per phase of 16A.

2,29A

In this specific case, the maximum 𝑃 𝑓 can be equal to 3 times the apparent power for 1 phase. In the worst case 𝑃 𝑓 maximum value can reach 11kW, which is near to the value obtained in

Table -1.

For the neutral current, in that case, even if the fluctuating power is high, its value remains low. This current can reach 3 times the nominal current (48A) in the worst case for the neutral current. This configuration shows the independence between the value of the fluctuating power and that of the neutral current

Classical design capacitor

The classical method to design the DC bus capacitor is first to determine the maximum voltage fluctuation which enables the system to function properly. Then by using equation ( 3) the value of the capacitance can be determined:

𝐶 𝑑𝑐 = max(𝑝 𝑓 (𝑡- 𝑇 𝑔 8
))

2𝜔 𝑔 𝑉 𝑑𝑐 max(𝑣 𝑎𝑐 (𝑡))

(

With 𝑝 𝑓 (𝑡) the fluctuating power at a frequency of 2𝜔 𝑔 , 𝑉 𝑑𝑐 the DC bus voltage and 𝑣 𝑎𝑐 (t) the allowed fluctuating voltage.

Minimal design capacitor

In the literature, there are many ways to design the bus capacitor depending on the hypothesis [START_REF] Hava | A DC bus capacitor design method for various inverter applications[END_REF]. In this part, the only constraint used to minimize the capacitance value, is the inverter saturation corresponding to minimal voltage needed at its DC input.

𝑣 𝑐 (𝑡) = 𝑉 𝑑𝑐 + 𝑝 𝑓 (𝑡 - 𝑇 𝑔 8 ) 2𝐶 𝑑𝑐 𝑉 𝑑𝑐 𝜔 𝑔 (4)
By definition, the modulation depth can be expressed as the ratio of the amplitude of the DC input to the peak-to-peak amplitude of the AC output of the inverter. Thus, per phase, the modulation depth is:

α x (t) = 2𝑣 𝑜𝑛𝑑,𝑥 (𝑡) 𝑣 𝑐 (𝑡) (5)
To simplify the equation and the resolution, the 𝑣 𝑜𝑛𝑑,𝑥 (𝑡) magnitude is expressed by 𝑉 𝑜𝑛𝑑,𝑥 ̂ in the following equations and the variable 𝑣 𝑐 (𝑡) is replaced by its expression (3):

𝛼 𝑥 (𝑡) = 2𝑉 𝑜𝑛𝑑,𝑥 (𝑉 𝑑𝑐 + 𝑝 𝑓 (𝑡 - 𝑇 𝑔 8 ) 2𝐶 𝑑𝑐 𝑉 𝑑𝑐 𝜔 𝑔 ) (6)
Using this equation, the minimum capacitance value can be determined having in mind that the modulation depth must be kept under 1 [START_REF] Eren | A Digital Current Control Technique for Grid-Connected AC/DC Converters Used for Energy Storage Systems[END_REF] to avoid the saturation. The expression of the minimum capacitance value is given by equation ( 8).

1 ≥ max(𝛼 𝑥 (𝑡)) (7)
C dc ≥ max(𝑝 finv (𝑡)) 2𝜔 𝑔 𝑉 𝑑𝑐 . 1

𝑉 𝑑𝑐 -2 𝑉 𝑜𝑛𝑑,𝑥 ̂ (8)
This equation could be more precise by taking account the phase shift between the fluctuation voltage and the AC voltage, but the improvement would be negligible. Furthermore, the determination of the capacity value does not need to be so precise due to the accuracy of these components.

In this application, we will consider an inverter with a DC bus voltage of 800V and a maximum fluctuating power of 10.8kW. In this case, to choose a constraint on the voltage fluctuation, the peak-to-peak voltage amplitude on the AC side is known, which is about 650V. Therefore, one can choose to use a minimum DC bus voltage that allows to be above this value. For example, the fluctuating DC bus value should not be below 680V, which allows a fluctuating voltage of 120V amplitude. The result of the dimensioning of this configuration is represented on This capacity value cannot cover all the cases given by the green curve. It only covers 95% of the cases. This means that 5% of the noncovered cases can lead to saturation of the inverter. This method shows that the minimum capacitance value cannot be deduced simply by imposing a voltage fluctuation, but needs to take also into account the minimum voltage linked to the operating points and unbalanced power exchanged.

Average model

It is desired to control the inverter so that it can be used in unbalanced three-phase power configurations without error. For this purpose, a corrector must be implemented to avoid any static error. Before dimensioning the corrector, it is necessary to know the transfer function of the system under study. For this purpose, an average model is used and demonstrated. The average model used is presented in Fig. 7.

Afterwards, the equations from average model are synthetized in state space model represented in equation ( 9 Besides, the inductance and resistance matrix have for expression: For the application, only 𝑖 𝑎 /𝑉 𝐴𝑁 , 𝑖 𝑏 /𝑉 𝐵𝑁 and 𝑖 𝑐 /𝑉 𝐶𝑁 are used. To reduce the error, feedback is used. In the case of an unbalanced three phase control, in the Park domain, it is necessary to use proportional integrator and resonant (PI-R) corrector [START_REF] Eren | A Digital Current Control Technique for Grid-Connected AC/DC Converters Used for Energy Storage Systems[END_REF]. Therefore, the controller computation is not simplified versus a standard time domain. To simplify the design, a PR corrector [START_REF] Abdulbaqi | Modeling and Analysis of a Four-Leg Inverter Using Space Vector Pulse Width Modulation Technique[END_REF] [9], [START_REF] Demirkutlu | A scalar resonant-filter-bank-based output-voltage control method and a scalar minimumswitching-loss discontinuous PWM method for the four-leg-inverter-based three-phase four-wire power supply[END_REF] is directly used with the following transfer function:

𝐴 𝐿 = [ 𝐿 𝑛 +
𝑃𝑅(𝑠) = 𝐾 𝑝 + 2𝐾 𝑟 𝜔 𝑐 𝑠 2 + 2𝜔 𝑐 + 𝜔 50 2 (18)
The value used for this transfer function is summarized in the Table -2:

Simulation

The implementation of the controller and its effectiveness has first been tested on Matlab/Simulink with a commutated model to represent the converter.

The DC capacitor takes have the value of 340𝜇𝐹. The Fig shows the result for depth modulation. This modulation doesn't fluctuate above 1, which is the limit of saturation. Furthermore, the error in power of PR control is below 1%. It shows that this control remains very precise, add confirms its use.

The Table -3 shows the difference between order and measure. It is desired to control the inverter so that it can be used in unbalanced three-phase configurations without error. For this purpose, it is seen in the literature that a PR corrector can compensate for these errors. In addition, in order to dimension this PR corrector, it is first necessary to know the transfer function of the system under study. For this, a modelling of the mean model is used and demonstrated. The average model used is shown in In order to measure accurately the impact of the fluctuating power induced by the unbalance operation on the DC voltage bus and avoid any interaction with the power supply, an additional inductance (LDC see Figure 12) has been added between the power supply and the DC bus. Its value is given in Table -4. The inductance has been designed in order to present an impedance almost 100 times higher than the capacitive DC bus at the 100Hz frequency. Thus, allowing to ensure that the fluctuating power is effectively passing through the capacitance and therefore have accurate results. In order to maximize the ratio of the fluctuating voltage versus the DC one, the latter has been adjusted in order to have a depth modulation close to 1. The 𝑉 𝑑𝑐 value is predetermined from simulation and applied on the bench test. The approach remains the same.

Exploitation of the results

Several operation points have been tested. The Table -5 shows the results obtained theoretically and experimentally. The Table -5 shows the theorical model accuracy for extreme operating point. The operating points with a low fluctuated power have a lower accuracy than the other points. In addition, for the same active power transferred the fluctuated power can be above or below the test 1 which should be one the most unbalanced case for the inverter. Test 1 corresponds to the single-phase mode for a three-phase inverter, which should correspond to one of the most unbalanced cases. Finally, configuration 2 with the same active power shows that the fluctuating power is much higher. Thus, this configuration is more unbalanced than the single-phase case. Test 3 shows that by having the same powers as case 2 but reversing 2 phases. The fluctuating power becomes almost non-existent. Finally, case 4 shows a case providing the same active power as the other tests, with a distribution of this power between two phases, with the addition of reactive power exchange. It is obtained a fluctuating power 2 times lower than in the single-phase case.

Conclusion

The proposal of this dimensioning of the capacity on the DC voltage bus, allows to predict, according to the use of the inverter, the minimum capacity to be used to avoid the saturation of the inverter.

The simulated and experimental study validated this approach. However, the inaccuracy error is due to the use of a simple model. This simplification does not take into account certain phenomena such as the voltage drops induced by the conduction of the transistors. This reduction in capacitance makes it possible to reduce the size of the components, their cost, and even to change the capacitor technology.

The only condition to allow such a fluctuation on the DC bus is that the equipment connected to this bus is able to operate with this disturbance, using a control law adaptation and/or an active filter.

The capacitance value chosen, at the end of the sizing study, must take into account many more constraints, such as EMC, filtering or the voltage withstand of components such as transistors. By superimposing all these constraints, the most important value of the capacitance must be kept. It has been seen that the energy router operation of an inverter cannot take all configurations because of the limit on the neutral current. Knowing that the application of energy routing is in the context of an inverter grid, it is possible to use a higher level control on the principle of hierarchical control [START_REF] Han | Review of Active and Reactive Power Sharing Strategies in Hierarchical Controlled Microgrids[END_REF] .The purpose of this control would be to allocate viable configurations to the different inverters while responding to grid services. This control ensures that the neutral current does not exceed its nominal value.

Fig 1 :Fig 2 :

 12 Fig 1 : Application of a multiport converter, splited in two converters; the QAB and the inverter, for energy routing in micro-grid application

  Therefore, this algorithm allows to determine the different possible cases and to evaluate the worst configuration. Thanks to the limitation in neutral current value (Inmax = Inom), it is possible to reduce the number of cases to design the capacitor. The Fig 3 shows the percentage of configurations which respect a current neutral value under its nominal value.

Fig 3 :

 3 Fig 3: The percentage of cases below the value of the neutral current divided by the rated current

Furthermore, this neutral

  current constraint is not linked to the fluctuated power constraint which will allow to design the capacitive DC bus. Indeed, the two Fresnel Diagram shown in Fig 4 indicate the current and power of each phase.

Fig 4 :

 4 Fig 4: Fresnel Diagram for current and fluctuated power On the Fig 4, the active power can be projected on direct and quadrature components, which allows to determine the magnitude and phase of fluctuated power that can be expressed with the following equation: 𝑝 𝑓 (𝑡) = 𝑃 𝑓 cos(2𝜔 𝑔 𝑡 -𝜙 𝑖 ) (2)

Fig 5

 5 shows the limit voltage shape wave base on this hypothesis It shows the worst configuration possible. Taking into account the fluctuating power, it induces a fluctuating voltage on the DC bus. The expression of the capacitive bus voltage is represented in equation (4):

Fig 6

 6 by the purple dot.

Fig 5 :

 5 Fig 5: Voltage wave form from DC side and AC side. Case saturation limit when 𝑚𝑖𝑛 (𝑉 𝑐 (𝑡)/2) = 𝑚𝑎𝑥 (𝑉 𝑜𝑛𝑑1 (𝑡))

Fig 6 :

 6 Fig 6: Comparison between standard method and specific method to design the capacitance for several value of DC voltages, in function of percentage of possible cases from energy router configuration.

Fig 7 :

 7 Fig 7: Average model from 4 legs and 2 level inverter

  Fig 7.

Fig 8 :Fig 10 :Fig 9 : 4 Experimentation 4 . 1

 8109441 Fig 8: Transfer function of inverter with and without PR corrector

Fig 12 :

 12 Fig 12: Electrical circuit of Bench test

Table - 2

 - : Coefficient values for PR corrector

		VALUE
	𝑲 𝒓	20
	𝑲 𝒑	0.1
	𝝎 𝟓𝟎	314,5 𝑟𝑎𝑑/𝑠
	𝝎 𝒄	10𝑟𝑎𝑑/𝑠

This enables to have the transfer function represented in Fig 8. Both transfer function are represented, without and with PR(s) corrector.

Table - 4

 - : Capacitance and inductance value for LC filter

	𝐿 𝑑𝑐	720mH
	𝐶 𝑑𝑐	0.340µF
	On the AC connection, an inductive filter with a
	value of 5mH is added to cancel current
	harmonic on each of the four legs. Finally, the
	test bench corresponds to the one presented on
	Fig 12.

Table - 5

 - : Test results from bench test with different unbalanced configuration


Damien Lemaitre
email: damienlemaitre@grenoble-inp.fr

Benoît Sarrazin

Yohan Wanderoild

Alexis Derbey

Thierry Brincourt

Yves Lembeye

Method to Increase WPTS Robustness to Frequency Splitting and Bifurcation Phenomena

Wireless Power Transfer System (WPTS) for electrical vehicle (EV) charging presents many advantages, but the classical values of coupling and load resistance can lead to frequency splitting phenomenon (FSP) and bifurcation phenomenon (BP). The latter affect the system performances due to a potential lower power factor, loss of the zero-voltage switching (ZVS), high sensitivity to over coupling. The present paper provides a method for selecting capacitor values to avoid/reduce these phenomena, specifically allowing to preserve ZVS until an over coupling 67 % higher than classical case.

Introduction

The development of electric vehicles (EVs) is leading to a growing need for chargers. EVs are largely charged with plug-in chargers. The safety of no exposed conductor, the galvanic isolation, the resistance to vandalism, weathering and the user comfort make wireless power transfer system (WPTS) attractive. The principle of the WPTS is presented in Fig. 1. The electrical circuit of the system is given in Fig. 2 considering a seriesseries (SS) resonant topology. Maximising efficiency is one of the main objectives of WPTS to compete with plug in chargers. One of the main concerns is how to design the coils for a desired power with maximum efficiency considering the all physical parameters, i.e. airgap, compensation topology, coil shape, etc.

In addition, the design procedure becomes more complicated considering Zero Voltage Switching (ZVS), bifurcation and splitting frequency phenomena. The latter are defined in part 2, where their influence on the system operation is studied. Part 3 proposes a study to measure the accentuation of bifurcation and frequency splitting according to the values of the coupling coefficient and the load resistance. Considerations used later for the design of the system are defined. Part 4 proposes a new method to reduce these phenomena and to increase the robustness of the system to coupling and load resistance variations. These results are applied and validated experimentally in part 5. Looking for highly efficient WPTS lead to consider bifurcation phenomenon (BP) and frequency splitting phenomenon (FSP). These phenomena are however defined differently in the literature, sometimes assimilated, sometimes confused.

The section below provides the definitions of BP and FSP considered in this article. The influence of these phenomena on the system will then be presented.

Definitions

To begin with, it is necessary to introduce the input phase angle (IPA) φ between the voltage V1 and the current I1 defined in Fig. 2. The frequency where the phase shift φ is equal to zero is called "zero phase angle" (ZPA) frequency. The power factor on the primary side of the coupler PF is the cosine of φ as defined in equation (1). The system therefore has a maximum power factor of 1 when operating at the ZPA frequency. The impedance seen from the source is in fact purely resistive.

PF = cos (φ )

(

Few articles address BP and FSP simultaneously.

Several articles in the literature present the bifurcation phenomenon as the occurrence of several ZPA frequencies [START_REF] Namadmalan | Design Optimization of Inductive Power Transfer Systems Considering Bifurcation and Equivalent AC Resistance for Spiral Coils[END_REF]- [START_REF] Moon | Analysis and Design of a Wireless Power Transfer System With an Intermediate Coil for High Efficiency[END_REF]. The term "bifurcation" is sometimes used to refer to another phenomenon describing several maximums of the output power as a function of the switching frequency [START_REF] Niu | Exact Analysis of Frequency Splitting Phenomena of Contactless Power Transfer Systems[END_REF], [START_REF] Kim | Mitigation of Frequency Splitting Phenomena Using a Matching Capacitor in Wireless Power Transfer System for Automated Guided Vehicle[END_REF].

The occurrence of multiple ZPA frequencies and multiple power peaks is very often linked as we will see in part 3. However, one phenomenon can clearly appear without the other as illustrated in [START_REF] Truong | Further Results on "Design Guidelines to Avoid Bifurcation in a Series-Series Compensated IPTS": Theoretical Analysis and Experimental Validations[END_REF]. It is therefore necessary to distinguish these two phenomena.

Following the definition proposed in [START_REF] Truong | Further Results on "Design Guidelines to Avoid Bifurcation in a Series-Series Compensated IPTS": Theoretical Analysis and Experimental Validations[END_REF], bifurcation will here refer to the appearance of several output power peaks as a function of the switching frequency, see Fig. 3 (a). The term "frequency splitting" will therefore refer to the occurrence of several ZPA frequencies for the IPA, see Fig. 3 (b).

Influence on System Performances

The BP and FSP can be more or less pronounced as shown in Fig. 3 (a) Fig. 3 (b). The curves in these figures show the output power and the IPA φ as a function of switching frequency for the same system, sized at a resonant frequency f of 81 kHz, for a rated power of 20 kW. The coupling coefficient is different in the three cases shown. The grey, orange and blue curves are obtained for coupling coefficients k of 0.1, 0.135 and 0.25 respectively.

If the BP is sufficiently pronounced, it is possible that the maximums of the output power no longer reach the rated power as illustrated with the blue curve (k=0.25) in Fig. 3 (a). This is a major drawback.

The FSP is problematic because it can cause the loss of the ZVS as illustrated in Fig. 3 (b) by the blue curve (k=0.25). The ZVS mode for the primary side inverter is ensured when operating above the resonant frequency [START_REF] Ferrieux | Maximising transferred power and preserving zero voltage switching in grid to vehicle and vehicle to grid modes of a wireless charging system[END_REF]. This requires an inductive impedance seen by the inverter that can be expressed by the condition from equation ( 2).

(𝑎) (𝑏)

FSP induced several zero-phase angle (ZPA) frequencies and reduced the input phase angle (IPA) value at the nominal point, that can become lower than the minimum value to guarantee ZVS mode. ZVS improves WPTS efficiency and protect the switch from hard switching that can bring over heat if designed for ZVS only. If ZVS can be obtained under FSP and BP, these last make the IPA much more sensitive to coupling variations. The ZVS is then only possible for a very small range of coupling coefficient k. The FSP also induces a decrease in efficiency through a higher IPA at rated power. In addition, several control strategies can become inefficient or less effective due to FSP [START_REF] Aditya | Design Guidelines to Avoid Bifurcation in a Series-Series Compensated Inductive Power Transfer System[END_REF], [START_REF] Wei | A SiC-Based 22kW Bi-directional CLLC Resonant Converter with Flexible Voltage Gain Control Scheme for EV On-Board Charger[END_REF].

The solutions provided in the state of the art to avoid these phenomena are described in the following section.

Criteria to avoid frequency

The electrical diagram in Fig. 4 is equivalent to the diagram in Fig. 2, considering the expression of Rs defined in [START_REF] Chopra | Analysis and design considerations for a contactless power transfer system[END_REF]. The primary quality factor Q and secondary quality factor Q of the system according to the definitions in [START_REF] Wang | Power transfer capability and bifurcation phenomena of loosely coupled inductive power transfer systems[END_REF] are given by the equations ( 4) and ( 5) respectively.

R = 8 π . R (3)
Fig. 4: Equivalent electrical circuit of WPTS.

Q = L . R w . M (4)
Q = w . L R (5)
While the FSP occurrence can be clearly defined through the transition from one to several zero IPA frequencies, BP apparition is harder to define due to a progressive appearance. This is why the criteria for avoiding the FSP as shown in Fig. 3 (b) are predominant in the literature. These are summarised in table 1 of [START_REF] Hatchavanich | Secondary-Side Voltage Control via Primary-Side Controller for Wireless EV Chargers[END_REF] for systems with a SS resonant compensation topology. It should be noted that some of these criteria for avoiding the FSP, contained in equations ( 6), ( 7) and ( 8), are equivalent. The coupling coefficient between the primary and secondary windings of the system is noted k.

The critical coupling coefficient above which the FSP occurs is denoted k .

Q < 1 2. 1 -√1 -k (6) k < k = 1 Q 1 - 1 4. Q (7) Q > 4. Q 4. Q -1 (8)
These criteria ensure that the operation of the switches in ZVS mode is maintained at the nominal operating point for which the system has been designed. However, they severely restrict the coupling coefficient in the coupler design, which impacts the system efficiency. Furthermore, choosing to maximize the coupling by taking the coupling coefficient k equal to k defined in equation ( 7) can be a trade-off between maximizing the coupling coefficient and limiting the accentuation of the BP and FSP.

Influential Parameters

The influence of the coupling coefficient and load resistance on the BP and FSP is carried out on a 20 kW WPTS, dimensioned for a resonant frequency of 81 kHz with an air gap of 21 cm. This coupler is presented here [START_REF] Sarrazin | Bidirectional Wireless Power Transfer System with Wireless Control for Electrical Vehicle[END_REF] and in part 5. The latter has been designed for a coupling coefficient k = 0.155, considering primary L1 and secondary L2 inductances of 338.74 µH and 80.66 µH respectively. The theoretical capacitances C1 and C2 calculated according to equation ( 9) at the resonant frequency are therefore 11.40 nF and 47.86 nF respectively.

C = 1 4. π . f . L (9)
The coupling coefficient of the system obtained under the design conditions, i.e. an air gap of 21 cm and with the windings aligned is noted kd.

It should be noted that the system has been designed such that k = k defined in equation ( 7) in order to maximize the coupling and limit the BP and FSP.

Coupling Coefficient

The influence of the coupling variation on the power transfer and the input phase angle is clearly illustrated in Fig. 3 (a) Fig. 3 (b). Over-coupling defines the case where k > kd. The larger the coupling k is in front of kd, the more pronounced the BP and FSP are. In contrast, a coupling k < kd reduces these phenomena.

Charge Resistance

The load resistance is a second parameter with a strong influence on the BP and FSP. The voltage across the EV battery represented by the load resistance Rch varies during charging. Considering that the power transferred to the battery Pout remains constant during charging, this implies that the resistance Rch varies with the charging voltage. The variation of Rch does not change the power factor PF . However, a charging voltage of 350 V, 400 V and 450 V implies charging resistances of 6.1 Ω , 8.0 Ω and 10.1 Ω respectively for a constant charging power of 20 kW. Fig. 5

Method proposed

The choice of the capacitance of the primary and secondary resonant capacitors is classically made from equation [START_REF] Truong | Further Results on "Design Guidelines to Avoid Bifurcation in a Series-Series Compensated IPTS": Theoretical Analysis and Experimental Validations[END_REF]. It is often difficult to obtain in practice the exact values of C1th and C2th due to the limited values of capacitors on the market. The output power vs. frequency curve can vary greatly depending on the choice made for a few percent difference in capacitance as will be detailed later. The following section proposes to add an additional capacitor selection criterion that provides a significant gain in system robustness to the BP and FSP. The robustness of the system gained by the proposed method will be studied analytically on the 20 kW system introduced in part 3.

Description

The proposed sizing method is divided into two steps. In order to maximize the efficiency while avoiding having pronounced BP and FSP, the first step is to perform the sizing of the WPTS with a coupling not lower than but equal to the critical coupling kFSP given by equation [START_REF] Niu | Exact Analysis of Frequency Splitting Phenomena of Contactless Power Transfer Systems[END_REF]. Inductances and resonant capacitances values L1, L2, C1, C2 are thus obtained on primary and secondary side at the desired resonant frequency fr, equal here to 81 kHz. The second step consists in a slight change of resonant capacitance values C1th, C2th. The new capacitance values thus obtained are noted C1c, C2c and are selected according to the classical criterion n°1, to which is added criterion n°2 below. The parameter Xc must be close to 1.

 Criterion n°1: C1 ≈ C1 & C2 ≈ C2  Criterion n°2: C1 = X . . C2 , with X ≠ 1
Criterion 1 ensures a good power factor, necessary to maximise efficiency. Criterion n°2 comes from inequation (10) below and allows the correction of the double resonance seen by the reflected impedance at the primary in the case of BP and FSP.

The influence of the parameter Xc is detailed in the following part, noting that Xc = 1 is the classical case. In the under-coupling case where k < kd, the BP and FSP attenuate when the coupling coefficient decreases, as illustrated in the case k = 0.11 Fig. 8 (a) and Fig. 8 (b). The inconvenience of reduced coupling is related to the lower power factor, see table 2. (𝑎) (𝑏) (𝑎) (𝑏) The previous theoretical study is validated in this section through the realisation of a 20 kW WPTS in the case Xc > 1. After presenting the system configuration, the output power versus frequency curve will be compared to the curve with Xc > 1, presented in the theoretical analysis in Fig. 7 (a).

Setup

The WPTS tested and shown in Fig. 9 has been sized for a nominal power of 20 kW, at a resonant frequency of 81 kHz. Full system characteristics are available in paper [START_REF] Sarrazin | Bidirectional Wireless Power Transfer System with Wireless Control for Electrical Vehicle[END_REF], [START_REF] Lemaitre | Efficiency-Oriented Design of Litz Wire for Several kW Power Experimented on 20 kW Prototype[END_REF]. As indicated in section 4.

Results

The relative error presented in Table 3 is the one between the measured output power of the prototype presented on Fig. 9 and the one calculated from the analytical model. Xc = 1.041 in both cases.

The theoretical and measured output power versus switching frequency curves are shown in Fig. 10. This graph is a zoom between 76 kHz and 83 kHz of Fig. 6 (a) of the blue theoretical curve (with Xc > 1). The relative error between the theoretical and measured values is a maximum of 4.2 % as shown in Table 3. The experimental results thus validate the analytical results presented above.

Conclusion

Bifurcation phenomenon (BP) and frequency splitting phenomenon (FSP) frequently appear in WPTS in the case of EV chargers, and are synonymous with additional losses, a control method rendered inoperative, and high sensitivity to the over coupling phenomenon. The impact of an increase in coupling and a decrease in load resistance on the accentuation of bifurcation and frequency splitting phenomena was presented. An easy to follow proposed method allows to avoid FSP and correct BP by slightly modifying the values of the resonant capacitances during the post-design. In particular, the proposed method makes it possible to reach the nominal power for a coupling coefficient of + 115 % and to maintain the ZVS for an over-coupling of + 67 % compared to the conventional case. The robustness of the system during coupling and load resistance variations is thus greatly increased without degrading the efficiency thanks to the preservation of a good power factor.

Fig. 1 : 2 :

 12 Fig. 1: WPTS principle of an EV charging Fig. 2: Schematic of a series-series resonant DC-DC converter

Fig. 3 :

 3 Fig. 3: Illustration of the bifurcation phenomenon (BP) (a) and of the frequency splitting phenomenon (FSP) (b) as function of switching frequency, resulting respectively in two maximums regarding output power and in several ZPA frequencies. These curves are obtained considering a designed WPTS of 20 kW at a resonant frequency fr equal to 81 kHz, with inductances and resonant capacitances values on primary and secondary side of L1, L2, C1, C2 respectively equal to 338.74 µH, 80.66 µH, 11.49 nF and 47.86 nF. Coupling coefficient considered for the grey, orange and blue curves are respectively 0.1, 0.135 and 0.25.

  (a) and (b) show respectively for a coupling coefficient kd = 0.155 the output power and IPA obtained for these three resistance values. The trend shows that the BP and FSP increase with decreasing Rch, i.e. at the beginning of the battery charge.

Fig. 5 :

 5 Fig. 5: Output power transferred (a) and Input phase angle (b) for the load resistances Rch equal to 6.1Ω, 8.0Ω, and 10.1Ω corresponding to grey, orange, and blue curves respectively. Power factor cos(φ ) is the same for all three cases and illustrated in brown. The coupling is k = 0.155 for the three cases. The WPTS studied has a nominal power of 20 kW, a resonance frequency of 81 kHz.

Fig. 6 :

 6 Fig. 6: Output power (a) and IPA (b) for a coupling coefficient k = 0.155 = kd. Three cases are presented: Xc > 1 (blue), Xc < 1 (orange), and the classical case (grey).

Fig. 7 :

 7 Fig. 7: Output power (a) and IPA (b) for a coupling coefficient k = 0.2 > kd. Three cases are presented: Xc > 1 (blue), Xc < 1 (orange), and the classical case (grey).

Fig. 8 :

 8 Fig. 8: Output power (a) and IPA (b) for a coupling coefficient k = 0.11 < kd. Three cases are presented: Xc > 1 (blue), Xc < 1 (orange), and the classical case (grey).

  2, the theoretical capacitance values are C1th = 11.40 nF and C2th = 48.86 nF, the selected capacitances are C1c = 12.5 nF and C2c = 50 nF, which corresponds to Xc = 1.041 > 1. This slight difference is sufficient to change the operating point at a given frequency. The electronics devices used for the measurements are current probes LeCroy CP150, voltage probe LeCroy HVD3605, WaveRunner RO 66Zi 600 MHz oscilloscope and thermal camera.

Fig. 9 :

 9 Fig. 9: 20 kW prototype realized with 21 cm air gap.

Fig. 10 :

 10 Fig. 10: Comparison of the measured and theoretical output power values between 12 kW and the 20 kW nominal output power in case of Xc = 1.041, for a coupling coefficient k = kd = 0.155.

Table 1 :

 1 Maximum coupling to guarantee the ZVS mode at Poutd and to be able to reach Poutd in the classical case and with the proposed method.

		Classical case	Proposed method	Difference
	k	0.18	0.3	+ 67 %
	k	0.2	0.43	+ 115 %

Table 2 :

 2 Comparison of WPTS performances in the classical case and with the case A of the proposed method. The capacity to transfer the nominal power P , as well as to operate in ZVS at P are investigated. The power factor cos (φ ) is given at P .

			Classical Case		Proposed Method : Case Xc > 1
		P	ZVS	cos (φ )	P	ZVS	cos (φ )
	k = 0.11	Yes	Yes	0.66	Yes	Yes	0.72
	k = 0.13	Yes	Yes	0.75	Yes	Yes	0.84
	k = k = 0.155	Yes	Yes	0.84	Yes	Yes	0.975
	k = 0.18	Yes	Boundary	0.99	Yes	Yes	0.986
	k = 0.2	Boundary	No	0.975	Yes	Yes	0.99
	k = 0.3	No	No	/	Yes	Boundary	0.99
	k = 0.43	No	No	/	Boundary	No	0.99

Table 3 :

 3 Comparison of the measured and theoretical output power values between 12 kW and the 20 kW nominal output power in case of Xc = 1.041, for a coupling coefficient k = kd = 0.155. The efficiency measured is presented.

	Frequency	Efficiency Measured Pout Measured	Pout Th	Relative error
	81 688 Hz	97,1 %	12 030 W	12 530 W	4,2 %
	81 448 Hz	97,0 %	13 037 W	13 310 W	2,1 %
	81 098 Hz	97,0 %	14 058 W	14 140 W	0,6 %
	80 568 Hz	97,0 %	15 059 W	14 990 W	0,5 %
	79 718 Hz	97,0 %	16 062 W	15 920 W	0,9 %
	78 618 Hz	96,8 %	17 064 W	16 940 W	0,7 %
	77 638 Hz	96,6 %	18 084 W	18 040 W	0,2 %
	76 908 Hz	96,4 %	19 091 W	19 140 W	0,3 %
	76 348 Hz	96,2 %	20 104 W	20 260 W	0,8 %

Study on 20 kW WPTS

The robustness of the system gained by the proposed method is studied based on the 20 kW system introduced in part 3. The comparison between the classical case corresponding to Xc = 1 and the proposed method is made according to three criteria: the ability of the system to transfer its nominal output power Pout , its ability to operate in ZVS mode, and its power factor. These three parameters are studied for different coupling coefficients k. The ZVS operation and power factor are calculated at the nominal output power for each value of k. The coupling coefficient varies according to the distance between the windings as well as their misalignment.

Capacitance values C1c and C2c of case Xc > 1 of the proposed solution is the one used during the experimental validation of the method in part 5, hence this choice for the theoretical comparison. The values chosen are C1c = 12.5 nF and C2c = 50 nF, which corresponds to Xc = 1.041. The second proposed case Xc < is performed for Xc = 0.94. In the classical case and the case with Xc > 1, the attainment of the Poutd power, the possibility of operating in ZVS and the power factor are shown for different coupling coefficients in Table 2. The results for Xc > 1 and Xc < 1. A first comparison is made under design conditions, with k = kd = 0.155. The output power and IPA φ as a function of switching frequency for k = kd are shown in Fig. 6 (a) and Fig. 6 (b). Two peaks exceed the rated power in the classical case, which is avoided in the Xc ≠ 1 cases. The arrows indicate the points at nominal operation. When several points allow the Pout specification power of 20 kW to be reached, the point that allows operation in ZVS, guarantees the best power factor and avoids passing a power peak is selected. In the case of over-coupling, i.e. k > kd, the increase of k in front of kd makes it increasingly difficult to meet two important conditions, i.e. operating in ZVS mode and being able to reach Poutd, see Table 2. This is illustrated in particular by Fig. 7 (a) and Fig. 7 (b) for which the coupling coefficient is 0.2. The system is in the classical case at the limit of being able to transfer Poutd while the ZVS mode is no longer possible while transferring Poutd. This is however still possible in the cases Xc < 1 and Xc > 1 as shown in Fig. 7 (a) and Fig. 7 (b


Kossigan Bernard Dedey

David Grenier
email: david.grenier@inrae.fr

Laurent Blondel

Yves Diascorn

Antoine Lejeune

Tiphaine Lucas
email: tiphaine.lucas@inrae.fr

Qualification of a small-sized bubble inflation device: gas tightness, uniformity of temperature, effect of the time initially left for dough to relax and numerical sensitivity analysis of stress-strain relationship to the initial disc diameter and thickness

This work was adpated

Introduction

Bread dough is often described as a dispersion of gas cells in a continuous hydrated gluten/starch matrix [START_REF] Gan | The microstructure and gas retention of bread dough[END_REF]. The crumb structure (specific volume and local gas fraction) of bread greatly depends on when and where gas cell walls (GCWs) separating these gas cells rupture during baking (Hayman et al., 1998a;[START_REF] Hayman | Factors controlling gas cell failure in bread dough[END_REF]. Crumb collapse with large cells is obtained for early GCW rupture. Shrinkage due to the lowering in pressure at the cooling stage (just after baking) occurs where a few or even none GCWs rupture. Rupture of GCWs at intermediate temperatures leads to a well-distributed cell crumb structure. At the end of proofing and during baking, part of the thinnest GCWs is reduced to about the size of a few starch granules [START_REF] Bloksma | Rheology of wheat flour doughs[END_REF][START_REF] Dobraszczyk | The rheological basis of dough stickiness[END_REF][START_REF] Sandstedt | The microscopic structure of bread and dough[END_REF]. These are the GCWs that are the most likely to rupture first. Although believed of great importance in controlling the expansion of bread, analysis of rupture at the scale of a thousand of starch granules has little received attention. Attempts to dynamically visualize dough extension at that scale on such thin films in order to identify where GCWs rupture are extremely scarce. The smallsized bubble inflation device that was qualified in this study makes it possible the control of the atmosphere surrounding a thin dough film, its temperature and the observation of the inflating thin dough film from the above when associated with a low-to-mid range digital microscope. Such microscope was not associated to the bubble inflation device in this study which only dealt with the analysis of the possible experimental and calculation biases relevant to the measurement of the stressstrain relationship. The objective of the work was to complement published studies dealing with the observation of the inflating bubble and the measurement of the stress-strain relationship on this device. In the published studies that refer to that work it was not possible to show all the preliminary characterization required for good measurement of stress and strain that accompany the visualization of the deformation of dough films. The preliminary characterization dealt with the analysis of the following biases which were experimentally or analytically investigated:

 Gas tightness of the small-sized bubble inflation device at 20°C and 70°C.  Uniformity in temperature.  Experimental and analytical analysis of the influence of:

o The time left to the dough film to relax its stress prior to inflation. o Diameter and thickness of the dough film.  Propagation of uncertainty on the calculation of stress, strain, thickness of the film, bubble height and radius.

2 Materials and methods

Sample preparation

245 g wheat flour (F1: Type 65, 0.53% ash, 10.58% proteins, from Moulins Soufflet Pantin, France), 144.6 g water, and 4.9 g salt (La Baleine, France) were mixed together using the kneader of the Alveolab ® (Chopin, France) at 60 rpm for 12.5 min. The initial water content (0.448 ± 0.002 wb) was determined from three dough samples after oven-drying at 104 ± 1°C for at least 24 h. The dough was removed from the kneader, stored in a sealed box, and allowed to rest for 15 min at least and 4 hours at most before testing.

Experimental procedure

A piece of dough of about 20 g was removed from the sealed box before each bubble inflation experiment. Thin films of dough were then produced by stretching the piece of dough by hands until the wanted thickness was reached. Two initial thicknesses (see sections 2.3.1 and 3.3.2.1) were experimentally investigated. The bi-extended piece of dough was then placed in the receiving clamp of the small-sized bubble inflation device. Two inner radii (𝒂) of the clamp (see sections 2.2.4 and 3.3.3.1) were experimentaly tested. The device was closed and the sample left to relax before bubble inflation. Two times of relaxation (see sections 2.2.4 and 3.3.1) were also experimentaly investigated.

Inflations were carried out on dough samples not systematically coming from the same preparation. Three replicates were carried out for each condition.

Bubble inflation

Dough bubbles were inflated within a home-made and small-sized bubble inflation device that consisted of two similar chambers placed one on top of the other (Fig. 1 #1 and #2) and separated by the inflating dough film. The dough film and the two chambers were surrounded by a double envelope in which thermoregulated water flows (Fig. 1 #3). Water is also flowing inside the double-glazed windows that were managed at the top and the bottom of the two chambers (Fig. 1 #6) for possible (not used in this study) dough lightening from the bottom and observation from above. The chambers were designed the same size (35 mm in height and 30 mm in diameter) in order to limit slight temperature inhomogeneities that would inevitably lead to evapo-condensation-diffusion and uncontrolled or biased hydration or drying of the dough film. A receiving clamp (Fig. 1 #5) consisting of two removable flat rings insured the holding of the dough film and gas tightness between the upper and lower parts of the dough film. Gas tightness between the clamp and the two chambers was insured using grease. A O-ring silicone flat gasket was placed between the chambers to ensure a tight seal between the two chambers using four wing nuts (Fig. 1 #7) that insured squeezing of the O-ring silicone flat gasket. A space of 100 μm was managed in both the upper and lower rings in order to only squeeze the piece of dough close to the outer rim of the sample. The shape of the rings of the clamp was designed in order not to expose sharp angles to the dough that could initiate undesired rupture of the dough bubble at the clamp level. Four dough tanks (Fig. 1 #4) were managed all around the device to be filled with dough in order to set the gaseous atmosphere and prevent the dough film from drying as much as possible.

Verification of gas tightness

No dough film inflation or dough tank were involved in this experiment. The empty small-sized bubble inflation device was first stabilized at 20°C or 70°C for 15 minutes and the two chambers were then sealed before increasing the pressure to 2 kPa using a syringe connected to the dedicated injection port on the side of the device (Fig. 1 #8). Both pressure and temperature were recorded. Three runs were carried out at each temperature.

Verification of the uniformity of temperature

For the study of the uniformity of temperature, the stretched dough film made by hand was placed within the clamp, a hole (< 1 mm in diameter) was managed by puncturing the middle of the dough film to insure gas to communicate chamber to chamber. The four dough tanks (Fig. 1 #4) were filled with extra dough (about 8 g of dough each) in order to stabilize the gaseous composition around the dough film. No gas was injected and the chambers remained in contact with the surrounding atmosphere through the small metallic ducts of the ports. The temperature was measured during a heating from 20°C to 70°C (at a rate of 5°C/min) that was then followed by a plateau at 70°C. Three runs were carried out at each temperature.

Experimental sensitivity study

The stretched dough film was placed within the clamp. The four dough tanks (Fig. 1 #4) were also filled with extra dough and the device was sealed and connected to the motorized syringe (see section 2.4.1) trough the injection port. The bubble was then inflated at constant strain rate at two extreme values relevant of bread dough during proving, baking and usual mechanical tests performed on dough (0.0006 s -1 and 0.09 s -1 ). The experiments were carried out either with the clamp of 7.5 or 10 mm in diameter and with either the initial thickness of 445 or 621 μm on average (see section 2.3.1).

Measurements

Initial thickness of the dough film

The initial thickness 𝑇(𝑡 = 0) of the dough film was measured using a digital caliper on dedicated samples (6 replicates) and dried for 24 hours at 104°C before measurement. The initial dry thicknesses were found to be 445 ± 143 μm and 621 ± 122 μm.

Temperature

When gas tightness was under study, the temperature was measured using four calibrated thermocouples (T-type, 0.2 mm in diameter, TC SA, France). T1, T2 and T3, T4 were located in the upper and lower chambers respectively (Fig2. a). The mean (Tm) of the four temperature measurements was then calculated in order to take into account the changes in gas pressure due to possible temperature variations in the calculation of the measured pressure (see section 2.4.2).

For the study of temperature uniformity, T1 was placed within one of the dough tank used for atmosphere regulation at the location where the dough was the most exposed to the gas in the chambers (Fig2. c). T4 was placed at the coldest spot of the dough tank i.e. the spot that was not directly in contact with the regulating bath, at the bottom of the dough tank (Fig2. c). T2 and T3 were placed at about 2 mm away from the dough film within both the upper and lower chambers (Fig2. b).

Gas pressure

A calibrated miniaturized pressure transducer of 9.5 mm in length (2BXCQ-093-25A, Kulite, USA) was connected to the metallic duct of the syringe at the lower chamber (Fig. 2 a). The upper chamber was connected to the outer atmosphere through a small metallic duct except for the study of gas tightness where it was closed. The pressure P used for stress calculation (section 2.4.3) was that measured in the lower chamber, it was expressed relatively to the atmospheric pressure. That pressure was corrected from temperature variation bias (section 2.4.2) when gas tightness of the device was under consideration.

Data acquisition

A data logger was used (Almemo®, Ahlborn, Germany) and the acquisition frequencies were 10 Hz or 1 Hz respectively when theoretical strain rates of 0.0006 and 0.09 s -1 were applied. For the study of gas tightness, the data were recorded at the frequency of 1 Hz for 1h30 which was about the maximum time of experiment for the very low strain rate, 0.0006 s -1 .

Alveolab ®

Bubble inflation was also carried out using the reference protocol AACC54-30 (Alveolab ® , Chopin technologies, France). The initial thickness of the dough sample used in the Alveolab ® was close to 5 mm. The stress-strain relationship was automatically calculated by the Alveolab ® software [START_REF] Dubois | AlveoConsistograph Handbook[END_REF]. The results were used for the comparison with the results obtained using the small-sized bubble inflation device.

Calculations

2.4.1 Control of constant strain rate: advance of the syringe plunger

The calculations of the stress, the height of the bubble (ℎ), the thickness of the bubble at the top (𝑇) (at any time t) at non-constant strain rate were previously presented by [START_REF] Bloksma | A calculation of the shape of the alveograms of some rheological model substances[END_REF] and [START_REF] Charalambides | Large deformation extensional rheology of bread dough[END_REF]. These equations were adapted in the present study for constant strain rate.

Let's consider the Hencky strain 𝜀 (Eq.( 1)).

𝜀 = 𝑙𝑛 1 + (1)
The strain rate 𝜀̇ Eq.( 2) is obtained by time derivating Eq.( 1).

𝜀̇= = (2)
By reorganizing Eq.( 2), it comes the differential equation Eq.( 3).

ℎℎ ̇-̇ ℎ -𝜀̇ = 0 (3)
By making a change of variable in Eq.( 3) and assuming 𝑦 = ℎ with 𝑦̇= 2ℎℎ ̇, Eq.(3) then writes Eq.(4).

̇ - ̇ 𝑦 -𝜀̇= 0  𝑦̇- ̇ 𝑦 -𝜀̇= 0  𝑦̇-𝜀ẏ = 𝑎 𝜀̇ (4)
Eq.( 4) is a differential equation of order 2 with a second member. The solution of such an equation is the sum of the solution of two sub-equations: the homogeneous equation i.e. without second member and the affine equation. The solution of the homogenous equation Eq.( 5) writes Eq.( 6).

𝑦̇-𝜀ẏ = 0 (5)

𝑦 = 𝐾𝑒 ̇ (6)
where K is a constant which will be determined later on using the initial conditions.

The solution of the affine equation writes Eq.( 7).

𝑦 = -𝑎

The solution to Eq.( 4) is therefore Eq.( 8).

𝑦 = 𝐾𝑒 ̇ -𝑎 (8)
Eq.( 8) then becomes Eq.( 9) using the initial conditions. At t=0, y=0 and so 𝐾 = 𝑎 .

𝑦 = 𝑦 + 𝑦 = 𝑎 (𝑒 ̇ -1) (9)
The solution of Eq. ( 4) then writes Eq.( 10).

𝑦 = ℎ => ℎ = 𝑎 𝑒 ̇-1 (10)
As already mentioned the knowledge of the height of the bubble (ℎ) as a function of time for a given strain rate makes it possible to calculate the plunger advance (𝑑(𝑡)). Eq.( 11) can be deduced using the assumption that temperature and pressure in the bubble inflation device and the syringe are equal.

𝑑 = 𝑎 𝑒 ̇-1 2 + 𝑒 ̇ /𝑆 (11)
Where 𝑆 was the surface of the plunger of the syringe. The advance speed (𝑣) of the plunger of the syringe (Fig. 3 #1) was used for programming the rotation speed of the stepper motor (3 axis CNC control board, stepper motor, 2-phase 4-wire, Huizhou Bachin Electronic Technology Co., Ltd, China) (Fig. 3 #2) Eq.( 12).

𝑣 = (12)

The use of two motor reductions (Fig. 3 #3) made it possible to reach a large span of strain rates. The programming was managed in Arduino language using the software dedicated for programming the Arduino Uno card (Uno, Arduino, Italy) and its associated motorshield interface (Motor Shield, Arduino, Italy) (Fig. 3 #4).

Correction of the measured pressure

The measured pressure was corrected from the induced pressure due to variation in temperature of the air within the device when the verification of gas tightness was carried out. Without this correction a decrease in pressure could have been interpreted as a gas leak when it was only due to a slight decrease in temperature. Under the assumption of constant volume (V) within the two chambers and no gas leaks, the perfect gas law reduces to Eq.( 13).

= 𝑐𝑡𝑒 (13)

The variation of absolute pressure 𝑃 _ due to the variation in temperature Tm writes Eq.( 14).

𝑃 _ = 𝑃 (14)
where 𝑃 was the absolute initial pressure in the chambers and 𝑇 the initial temperature (K). The variation in relative pressure (𝑃 ) due to the variation in temperature then writes Eq.( 15).

𝑃 = 𝑃 _ -𝑃 (15)
The pressure relative to the initial pressure without the effect of temperature variation in the chambers 𝑃 is given by Eq.( 16).

𝑃 = 𝑃 -𝑃 -𝑃(𝑡 = 0) (16)
where 𝑃(𝑡 = 0) was the relative pressure at t = 0 s when the device was sealed.

Stress

The pressure 𝑃 measured in the lower chamber was used for the calculation of the stress 𝜎 at the top of the bubble Eq. ( 17).

𝜎 =

.

(

) 17
where 𝑅 was the radius of the bubble (Eq. ( 18)) and where 𝑇 was the thickness at the top of the bubble at any time (Eq.( 19)) [START_REF] Bloksma | A calculation of the shape of the alveograms of some rheological model substances[END_REF][START_REF] Charalambides | Biaxial deformation of dough using the bubble inflation technique[END_REF].

𝑅 = (18) 𝑇 = 𝑇(𝑡 = 0) 1 + (19)
where 𝑇(𝑡 = 0) was the initial thickness of the film.

Analytical sensitivity study

A numerical sensitivity study to the initial thickness of the bubble was first managed assuming that the gas pressure was constant and set at 0.65 kPa, a pressure close to the average of the pressure measured over all the duration of a bubble inflation experiment. The strain rate at the top of the bubble was set at 0.09 s -1 . The initial thickness of the film was varied from 100 to 1,500 µm, with a step of 100 µm, directly in the set of equations. The two average values of film thickness measured in the experiments 445 ± 143 μm, 621 ± 122 μm were also considered. The radius of the receiving clamp was varied (7.5 and 10 mm).

Then, the uncertainty analysis on the parameters used in the calculations was performed at the light of the uncertainty evaluated on each parameter. The Monte-Carlo method was applied and 500 inputs for each parameter were randomly selected within their uncertainty range. The uncertainties on the experimental parameters used for calculation, i.e. initial bubble radius (7.5 mm) and initial thickness (404 µm), gas pressure (set constant at 0.65 kPa) and the strain rate (0.09 s -1 ), were estimated at 5%, 30%, 5%, 30% respectively. At each step-time, each outputs were calculated for the random 500 inputs and the mean, the standard deviation and the relative error (standard deviation to mean ratio) were calculated.

3 Results and discussions 3.1 Gas tightness Fig. 4 shows the evolution of the pressure relative to the initial pressure (𝑃 ) against time at room temperature (Fig. 4 a) and at 70°C (Fig. 4 b).

At 20°C, the device was gas tight during at least 1.5 hour (Fig. 4 a). The decrease in pressure was lower than 5.10 -4 kPa/min. This slight leakage was judged not to affect the analysis of stress-strain relationship at 20°C as its impact was lower than 1% of the measured stress.

At 70°C (Fig. 4 b), a rise in pressure relative to the initial pressure during the first 1,000 seconds (blue rectangle, Fig. 4 b) was noticed after the closure of the injection device. It was followed by a pressure stabilization to a plateau (red line almost horizontal, Fig. 4 b). The rise in pressure during the first 1,000 s was due to the thermal expansion of the gas injected in the chamber which was initially at room temperature within the syringe. The variation in temperature of the gas injected at room temperature was not measured by the thermocouples and could not be taken into account in the correction of the measured pressure (section 2.4.2). For the next generation of devices, the syringe should also be thermoregulated. The device looked to be more gas tight on average at 70°C than at 20°C (See linear fit). This was probably due to the thermal dilation of the device that increased squeezing of the O-ring silicone flat gasket. The discrepancies between trials was however higher at 70°C (±1 kPa vs 0.5 kPa). Again, the impact of this leakage do not influence the measured stress in an extent greater than 1%. 3.2 Uniformity of temperature Fig. 5 presents the temperature evolution at the coldest spots (within the dough in a dough tank) and close to the punctured thin dough film maintained in between the lower and upper chambers (Fig. 2. b and c).

During heating and once stabilized at 70°C, a superimposition of temperatures measured by thermocouples T2 and T3 was noticed, which warrants a temperature homogeneity around the dough film better than 0.8°C.

The temperatures measured at the core of the dough tank (T4) was about 10-15°C lower than those close to the dough film (T2 and T3) for the first 1,000 s of heating (Fig. 5a). Then this difference decreased with the progressive establishment of the thermodynamic equilibrium within the device following the stabilization of temperature at 70°C or about. The difference of temperature between the core of the dough tank (T4) with the temperature close to the dough film was lower than 0.5°C (Fig. 5 b). The difference in temperature between the surfaces of the dough film and the dough tank is indeed the most determinant of the vapour migration within the device and measurements of temperature at the surface of the dough tank could usefully complete this analysis and demonstrate that the difference is even less than the above mentioned value. After 1,000 s of stabilization, the drying of the dough film due to evapo-condensation from hot spots where the vapour pressure was the greatest to the cold spots where it was the lowest was thought to be minimized. It was considered that this remaining difference in temperature was non-significant and did not allow intensive further drying of the dough film. Mechanical testing of the film was started after this time of stabilization (> 1,700 s).

Sensitivity of the stress to different factors

In the followings, when not specified, the reference values for the radius of the receiving clamp, the time left to the film sample to relax and the initial thickness of the film were a= 7.5 mm, τ = 45 min and T(t=0) = 445 𝜇𝑚 respectively.

Effect of the time initially left for dough to relax

Fig. 6 presents the stress at the top of the bubble versus the Hencky strain for the two resting times left to the dough sample to relax before bubble inflation within the small-sized bubble inflation device at 0.0006 s -1 and 0.09 s -1 . No significant effect of the resting time was evidenced up to 1.5 Hencky strain at the two extreme settings of strain rates experienced in this work. 45 min of rest is the time commonly used in the literature before mechanical testing of bread dough film [START_REF] Charalambides | Biaxial deformation of dough using the bubble inflation technique[END_REF][START_REF] Dobraszczyk | Extensional rheology and stability of gas cell walls in bread doughs at elevated temperatures in relation to breadmaking performance[END_REF][START_REF] Shewry | Measurement of biaxial extensional rheological properties using bubble inflation and the stability of bubble expansion in bread doughs and glutens[END_REF] and we have retained this value as a reference for comparison to other works, even if 10 min seem to be sufficient for complete relaxation in these experiments.

Fig. 6. Experimental stress vs. Henchy strain using the sealed small-sized bubble inflation device; effect of two resting times of dough (τ = 10 min and 45 min) at strain rates of 0.09 s -1 and 0.0006 s -1 .

3.3.2 Effect of the initial thickness of the film 3.3.2.1 Experimental Fig. 7 shows stress against Hencky strain at 0.09 s -1 for the two initial average thicknesses of the film, 445 μm (blue line) and 621 μm (red line). The greatest initial thickness was 1.4-fold that of the lowest. The greatest the initial thickness of the dough film, the lower the increase in stress with Hencky strain. At 1.5 Hencky strain, stress decreased about 2-fold from 240 to 120 kPa when the initial thickness of the film increased 1.4-fold (from 445 to 621μm). This makes hard to compare stress-strain curves obtained between reports in the literature as the slightest difference in initial thickness can change a lot the value of stress. The increase in stress observed on the small-sized device in comparison to that obtained by [START_REF] Dobraszczyk | Extensional rheology and stability of gas cell walls in bread doughs at elevated temperatures in relation to breadmaking performance[END_REF] was possibly explained by the extent of the initial strain. The current result supports the idea according to which the difference in stress could come from the difference in initial thickness. The lower the initial thickness, the greater the initial strain and the greater the ability of the dough for strain hardening. The results coming from the Alveolab ® for an initial film thickness close to 5 mm (solid red line) also showed that the same dough could exhibit Hencky strain, ε H 0.0006 s -1 0.09 s -1 different stress-strain relationship according to its initial thickness. Note however that the strain rate applied at the top of the bubble was not constant during the Alveolab ® test and not of the same order of magnitude as that applied in this study within the small-sized bubble inflation device.

Analytical

Fig. 7 presents the analytical stress vs Hencky strain under Bloksma's assumptions [START_REF] Bloksma | A calculation of the shape of the alveograms of some rheological model substances[END_REF] and the assumption of a constant pressure set at 0.65 kPa throughout the duration of the bubble inflation for an initial film thickness varying from 100 μm (the smallest thickness expectable by hand stretching) to 1.5 mm (the approximate thickness used in the experiments of [START_REF] Dobraszczyk | Extensional rheology and stability of gas cell walls in bread doughs at elevated temperatures in relation to breadmaking performance[END_REF] and [START_REF] Charalambides | Large deformation extensional rheology of bread dough[END_REF]). Consistently with the experimental trends discussed above, these results exhibited the fact that at great Hencky strain, the stress was very sensitive to the initial film thickness, in particular when one tries to work with dough films lower than 300 μm in thickness i.e. close to the average thickness of gas cell walls at the end of proving. The stress at the top of the bubble for a film of 100 μm in initial thickness was found about 4-fold that calculated for an initial thickness of 445 μm. This difference decreased as the initial film thickness increased with reference to 100 µm. This analysis shows how much it is important to directly measure strain or bubble thickness if one wants not to make too much error on the calculation of stress. Fig. 7. Effect of the initial film thickness varying from 100 μm to 1.5 mm on the analytical stress-strain (0.09 s -1 strain rate), 45 min resting time and a = 7.5 mm ; comparison to the experimental stress-strain curves obtained with initial film thicknesses (445 ± 143 μm and 621 ± 122 μm) using the small-sized bubble inflation device at 0.09 s -1 . The stress-strain curve obtained with the Alveolab ® using a dough film of about 5 mm thick initially and non-constant strain rate.

3.3.3 Effect of the radius of the receiving clamp 3.3.3.1 Experimental Fig. 8 a shows the stress vs Hencky strain for the two initial radii of the receiving clamp of 7.5 and 10 mm at the two settings of constant strain rates at the top of the bubble of 0.09 s -1 and 0.0006 s -1 . The effect of the initial radius on the mechanical behaviour of the dough film was low at extremely low strain rate of 0.0006 s -1 ; the curve at the greatest radius (10 mm) was 1.09-fold that at the smallest radius at Hencky strain of 1.5. The effect increased with the increasing strain rate at 0.09 s -1 , the stress was about 1.5-fold using the receiving clamp of the smallest radius (7.5 mm) compared to the greatest radius (10 mm). The greater the setting of the strain rate, the greater the effect of the initial radii of the dough film on the stress.

Analytical

Fig. 8 b presents the analytical stress vs Hencky strain for an initial thickness of 445 μm and the two radii of the receiving clamp (7.5 and 10 mm) for two settings of constant strain rates at the top of the bubble of 0.09 s -1 and 0.0006 s -1 . The stress for a radius of 10 mm was found to be 1.3-fold that when using a receiving clamp of 7.5 mm in radius at 1.5 Hencky strain. This result was in agreement with the experimental data at low strain rate (0.0006 s -1 ). However, it was the opposite of what was experimentally observed at 0.09 s -1 with a variable pressure (Fig. 8 a). This observation was probably due to the constant pressure assumption we considered in the case of analytical estimation of the stress while the real pressure measured was non-constant and greater when using the small diameter receiving clamp than when using the largest diameter. These results once again support the point that the knowledge of the experimental parameters such as pressure, bubble height and especially the initial thickness is crucial for a good estimation of the stress at the top of the bubble. 3.4 Sensitivity study to the uncertainty on all the entry parameters used in the analytical equations Fig. 9 a presents the mean analytical height and radius of the bubble (black lines), calculated over the 500 combinations of input parameters at each time step. The associated uncertainty is presented in dash red lines. The bubble height first decreased from the infinite (the radius of the initial flat dough film) to the value of 7.5 mm when the bubble turned from an 'open' to a 'closed' spherical cap where the radius, the height of the bubble and the radius of the receiving clamp were all equal. From that point on, the radius increased. In the meantime, the height first quickly increased to the value of 3 mm due to the constant strain rate condition imposed at the top of the bubble. Then it linearly increased up to 14 mm at 1.5 Hencky. The impact of the uncertainties on the bubble height was of the order of 2.5 mm for a mean height of 14 mm. The relative error on the bubble height reached about 18% (Fig. 10 green) and the relative error on the bubble radius reached about 10% (Fig. 10 black dot). Fig. 9 b presents the analytical thickness at the top of the bubble against time. It decreased from 404 μm to about 30 μm at 1.5 Hencky strain where the divergence with the hypothesis of non-uniformity of the bubble thickness started to be very significant. The associated relative error was found to be of the order of 52% at 1.5 Hencky strain (Fig. 10 red line). Such an error was quite significant and revealed the importance of being accurate on the input parameters if one still wants to estimate the thickness of the bubble in the domain were the equation are supposed to remain valid (< 1.5 Hencky strain). Fig. 9 c presents the analytical stress at the top of the bubble against time. Let's first remind that the applied pressure was uniform throughout all this analytical bubble inflation, even at the beginning of the inflation. This explains why the stress was not null at the beginning of the inflation. It then increased exponentially. The relative error reached a value of about 60% at 1.5 Hencky strain (Fig. 10 black solid), showing again how the calculation of the stress amplifies the errors on the input parameters. . Relative error on analytical outputs issued from the calculation under Bloksma's assumption with the assumption of uncertainties of 5% on the initial bubble radius, 30% on the initial thickness of the dough film, 30% on the strain rate at the top and 5% on the pressure measurement.

Conclusions and perspectives

Gas leak was observed at 20°C and 70°C but found not to be influencing the results significantly. One might stay confident about the tightness of the device and the reliability of the measured pressure and strain though some improvements in gas tightness should be made in the next generation of devices.

The greatest difference in temperature within the device was reduced as much as possible to 0.8°C once the temperature has stabilized to 70°C within the device. This give quite good confidence on the ability of the device to minimize water content variations of the dough film during long lasting mechanical tests involved at low strain rate. Without this precaution, any cold spots inside the device would condense water vapor within the device and bring a progressive dehydration of the sample. Remember that it is difficult to check the water content of the dough film after mechanical testing because of its extremely low thinness; this is true at 20°C and even more at 70°C.

The small-sized bubble inflation device remains a delicate tool to handle when one wants to estimate the stress-strain relationship. The equations used for the calculation of the stress are based on assumptions that quickly (Hencky strain > 1.5) become not accurate enough to allow precise identification of the stress. The structure of the equations used highly amplifies the slightest error on each input parameter. The stress-strain relationship was found to be very sensitive to the film thickness and especially when it became thinner than 300 µm.

Fig. 1 .

 1 Fig. 1. The small-sized two-chamber inflation device. 1 and 2: Upper and lower chambers 3: Connexions for the thermoregulated water circulation 4: Dough tanks (two per chamber) 5: O-ring Silicone flat gasket 6: Double glazing (one per chamber), 7: wing nuts.

Fig. 2 :

 2 Fig. 2: Location of the thermocouples for the verification of gas tightness (a), for the verification of the uniformity in temperature (b) and (c): (b) in the vicinity of the dough film and (c) inside one of the four dough tanks.

Fig. 3 .

 3 Fig. 3. (a) (1) The motorized syringe used for bubble inflation. (2) The stepper motor and its associated (3) motor reductor and the coupled controller and motorshield cards used for controlling the stepper motor and programmed for exponential movement of the plunger (constant strain rate at the top of the bubble). (b) Image of a bubble without the upper chamber. (c) Schematic of the inflating bubble (from Charalambides et al. (2002)).

Fig. 4 .

 4 Fig. 4. Relative pressure Pr within the chambers vs time, (a) at room temperature (b) at 70°C. All for 5,000s (about 1 hour 20 minutes). Mean and standard deviation between three runs.

Fig. 5 .

 5 Fig. 5. (a) Temperature vs time for thermocouples T2, T3 near the dough film, T1, T4 in a dough tank, (b) zoom in on the blue boxed area in Fig. 5 a. The temperature was set to 70°C. Mean and standard deviation between three runs.

Fig. 8 .

 8 Fig. 8. Stress vs Hencky strain for the two radii of the receiving clamp a = 7.5 mm and 10 mm and two strain rates of 0.09 s -1 and 0.0006 s -1 . (a) experimental data (45 min resting time and initial film thickness of 0.445 mm), (b) analytical data.

Fig. 9 d

 9 Fig. 9 d presents the analytical strain at the top of the bubble against time. It was directly impacted by the uncertainty on strain rate which could have a real value different from that set. The relative error was constant and of the order of 18% (Fig. 10 blue).

Fig. 9 .

 9 Fig. 9. Analytical (a) mean height and radius of the bubble vs time, (b) mean thickness of the bubble vs time, (c) mean stress vs Hencky strain and (d) Hencky strain vs time, black lines. The range of uncertainty (mean ± standard deviation) is shown in red dash lines.

  Fig.10. Relative error on analytical outputs issued from the calculation under Bloksma's assumption with the assumption of uncertainties of 5% on the initial bubble radius, 30% on the initial thickness of the dough film, 30% on the strain rate at the top and 5% on the pressure measurement.