A Pierquin
email: antoine.pierquin@ec-lille.fr

S Brisset

T Henneron

S Clénet

Benefits of Waveform Relaxation Method and Output Space Mapping for the Optimization of Multirate Systems

Keywords: Differential algebraic equations, Electromagnetic coupling, Finite element methods, Optimization methods

We present an optimization problem that requires to model a multirate system, composed of subsystems with different time constants. We use waveform relaxation method in order to simulate such a system. But computation time can be penalizing in an optimization context. Thus we apply output space mapping which uses several models of the system to accelerate optimization. Waveform relaxation method is one of the models used in output space mapping.

I. INTRODUCTION

I N the framework of the optimization of a multi-physics system, it is necessary to model the whole system and to perform a coupling of different numerical models. But modeling of a system including components with very different time constants is particularly problematic. On one hand, a strong coupling involves a time discretization according to the smallest time constant, and thus a large numerical system to solve and a long computation time. On the other hand, a weak coupling implies a lack of consistency of the results. In fact, a model used in optimization has to be as precise as possible but not too long to be computed because of the huge number of evaluations which the optimization process requires.

Waveform relaxation method (WRM) [START_REF] Lelarasmee | The waveform relaxation method for time-domain analysis of large scale integrated circuits[END_REF], [START_REF] Schöps | A cosimulation framework for multirate time integration of field/circuit coupled problems[END_REF] is an iterative process which allows to model each component of the multiphysics system with respect to its own time constant. This can reduce computation time, while keeping a good precision since the method converges to the exact solution [START_REF] Arnold | Preconditioned dynamic iteration for coupled differential-algebraic systems[END_REF], [START_REF] Crow | Convergence properties of the waveform relaxation method as applied to electric power systems[END_REF]. However, even WRM optimization can be extremely long to execute. With the aim of reducing more optimization time, an output space mapping (OSM) strategy [START_REF] Guofo | An output space-mapping algorithm to optimize the dimensional parameter of electromagnetic relay[END_REF], [START_REF] Bandler | Em-based surrogate modeling and design exploiting implicit, frequency and output space mappings[END_REF], [START_REF] Encica | Aggressive output space-mapping optimization for electromagnetic actuators[END_REF], [START_REF] Tran | Adapted output spacemapping technique for a bi-objective optimization[END_REF] can be set up. This is still an iterative process which requires at least two models of the same device, but with different accuracy and computation time. A coarse model, the fastest one, but the less accurate, is used during optimization. A fine model, the most time consuming, more precise, is evaluated once per iteration to correct the other model. Thus, the WRM can be used to produce the most precise model of the OSM.

The first two parts of this article present the waveform relaxation method and the output space mapping technique. In the last part, these methods are applied to the minimization of a transformer mass, using a finite element model (FEM).

II. WAVEFORM RELAXATION METHOD

WRM allows to model a multi-physics system of which components have heterogeneous time constants. The system is split to model each component with respect to its own time constant. Then the communication between subsystems is done by an exchange of waveforms. A loop is performed, with a relaxation at each iteration, until convergence.

Let a Differential Algebraic Equation (DAE) represent a system on the time domain T = [t 0 , t f ]:

ẏ(t) = h(y(t), z(t)), (1
) 0 = g (y(t), z(t)) .

(2)

y : [t 0 , t f ] -→ R m , z : [t 0 , t f ] -→ R p , f : (R m , R p ) -→ R m , g : (R m , R p ) -→ R p .
The system is decomposed into r subsystems. Each subsystem i satisfies:

ẏi (t) = h i (y(t), z(t)), (3) 0
= g i (y(t), z(t)), (4)
with

y(t) = [y 1 (t), . . . , y i-1 (t), y i (t), y i+1 (t), . . . , y r (t)] T , z(t) = [z 1 (t), . . . , z i-1 (t), z i (t), z i+1 (t), . . . , z r (t)] T . y i : [t 0 , t f ] -→ R mi , z i : [t 0 , t f ] -→ R pi , f i : (R mi , R pi ) -→ R mi , g i : (R mi , R pi ) -→ R pi , m = r i=1 m i , p = r i=1 p i .
Equation ( 3) is the differential equation of subsystem i, and (4) is the algebraic equation. In these equations, y are state variables and z coupling variables.

The WRM produces iteratively an approximation (ỹ k (t), zk (t)) of the solution (y(t), z(t)), where k is the iteration index.

The initial iteration is fixed using the known values of y and z at time t 0 . This is the extrapolation step:

ỹ0 (t) = y(t 0 ), z0 (t) = z(t 0 ), ∀ t ∈ [t 0 , t f ]. (5)
Then, at iteration k, the DAE (3)-( 4) are solved successively from subsystem 1 to r using Gauss-Seidel relaxation process:

ẏk i (t) = f i ( Ỹ k i (t), Zk i (t)), (6)
0 = g i ( Ỹ k i (t), Zk i (t)), (7)
with

Ỹ k i (t) = [ỹ k 1 (t), . . . , ỹk i-1 (t), ỹk i (t), ỹk-1 i+1 (t), . . . , ỹk-1 r (t)] T , Zk i (t) = [z k 1 (t), . . . , zk i-1 (t), zk i (t), zk-1 i+1 (t), . . . , zk-1 r (t)] T .
The algorithm stops when the norm of the difference between two successive iterations is smaller than a given tolerance.

In the case of an exact resolution of ( 6)-( 7), convergence of the WRM to the exact solution is proven [START_REF] Lelarasmee | The waveform relaxation method for time-domain analysis of large scale integrated circuits[END_REF], [START_REF] Arnold | Preconditioned dynamic iteration for coupled differential-algebraic systems[END_REF], [START_REF] Crow | Convergence properties of the waveform relaxation method as applied to electric power systems[END_REF]. In the case of numerical resolution, the convergence is always effective, but an error is introduced by the discretization schemes. Equations ( 6)-( 7) are solved numerically by using the most adapted time-scheme to the subsystem i, so each subsystem is solved using its own time discretization. We choose to use a linear interpolation between two discrete values of a waveform to obtain the value of this waveform at any time of T .

III. OUTPUT SPACE MAPPING

Computation time of an optimization process depends on the complexity of the model to be evaluated during the process. A precise model is often long to simulate, and conversely a fast model is less accurate. Space mapping techniques allow to perform a fast and precise optimization by using the advantages of both models.

The following optimization problem has to be solved:

x ⋆ f = arg min x f (x) such that k f (x) ≤ 0, (8)
with

f : X → R, k f : X → R q , X ⊂ R n .
Objective function f and constraints k f form the fine model, with both high precision and computation time. A second model of the same problem is considered: c and k c are the objective function and constraints of the coarse model, faster but less accurate. The coarse optimization problem associated is:

x ⋆ c = arg min x c(x) such that k c (x) ≤ 0, (9)
with c : X → R, k f : X → R q .
The principle of OSM is to optimize with the coarse problem, then to evaluate the fine model at the solution found to obtain correction of the coarse model. This process is applied iteratively. At the j-th iteration of the OSM procedure, we consider a corrector O j ∈ R for the objective function and a corrector Õj ∈ R q for the constraints. The corrected problem is:

x j = arg min x O j .c(x) such that Õj .k c (x) ≤ 0. ( 10
)
From the solution x j , new correctors are computed by evaluating the fine model, so the number of evaluations of the fine model is equal to the number of iterations. The algorithm stops when a convergence criterion is satisfied: for example, when

x j -x j-1

x j-1 is less than a given tolerance ε. We can also use a criterion on the difference between fine and coarse models. The algorithm is given:

1. j = 0 2. Initial point x 0 3. O 0 = f (x 0 ) c(x 0 ) , Õ0 i = k f,i (x 0 ) k c,i (x 0 ) , i = 1, . . . , q 4
. while

x j -x j-1 x j-1 > ε 4.1. x j = arg min x O j .c(x) s. t. Õj .k c (x) ≤ 0 4.2. O j+1 = f (x j ) c(x j ) , Õj+1 i = k f,i (x j ) k c,i (x j ) , i = 1, . . . , q 4.3. j = j + 1
OSM implies choosing two models: one coarse and one fine. In a system of components with heterogeneous time constants, the WRM is an adapted way to obtain a fine model with a shorter computing time than a strongly coupled model.

IV. APPLICATION

The OSM strategy is applied to the optimization of a transformer. We consider a device composed of a circuit supplying a transformer (Fig. 1(a)): a pulse width modulation (PWM) voltage source, an LC filter, a resistor and a transformer. Two models of this device are necessary to apply the space mapping.

A. Coarse and fine models

The coarse model is a circuit model of the device (Fig. 1(b)), where the transformer is represented by an inductance L = µ0µrN 2 S ℓ , with N the number of turns in the primary coil, S the section and ℓ the length of the magnetic core.

The fine model is a simulation by WRM where the system is decomposed into two subsystems (Fig. 1(c)): the circuit and the transformer. The circuit consists of the PWM voltage source with the LC filter and a current source. The transformer is modeled by 3D FEM (only one eighth of the transformer is modeled, Fig. 2) in vector potential formulation with a voltage coupling [START_REF] Bossavit | A rationale for edge-elements in 3-d fields computations[END_REF], [START_REF] Piriou | Finite element analysis in electromagnetic systems-accounting for electric circuits[END_REF]. The source of the circuit is the current i R in the resistor, and the source of the transformer is the voltage v c of the capacitor. In the WRM loop, at the k-th iteration, circuit part is solved with i k-1 R as a source. Voltage v k c is obtained and used as source for the transformer. This gives i k R , and the loop goes to the next iteration. For FEM, the magnetic flux density B is expressed with respect to the vector potential A:

B(χ, t) = curlA(χ, t), (11) χ ∈ D ⊂ R 3 , t ∈ T.
By Ampere's law and the coupling equation we obtain the following system:

curl 1 µ curlA(χ, t) -N (χ)i R (t) = 0, ( 12
) d dt D A(χ, t).N (χ)dv + R.i R (t) = v c (t). ( 13
)
where N = J /i R , with J the current density.

Last, time discretizations are different in the two subsystems. Because of the PWM, time-step for the first subsystem is dt 1 = 5 × 10 -7 seconds, whereas in the transformer part, time-step is dt 2 = 10 -3 seconds.

i L v C i R v PWM (a) Complete device. i L v C i R v PWM L (b) Coarse model. v C i R i L i R v C dt=5.0e-7 dt=1e-3 v PWM i R v C (c) Fine WRM model.

B. Optimization problem

The aim is to minimize the transformer mass m, and to impose RMS current value into the transformer. The design variables are width L and height H of the transformer: all other dimensions are deduced from this two length (Fig. 3). These two length form the optimization variables: x = [H, L]. We denote by i c and i f the RMS values of current i R , obtained respectively with the coarse and the fine model. The optimization problem is: OSM algorithm is applied to the optimization problem of the transformer. Objective function is the same for the two models, the correction is applied only on the constraint on the current in the primary coil. At iteration j of the OSM process, a corrector Õj = i f (x j ) i c (x j ) is applied to the value i c such that the following corrected problem is solved:

             min H,L m(H, L), 20cm ≤ H ≤ 40cm, 12cm ≤ L ≤ 24cm, H -2L 3 > 0, i f = 3 A. (14)
             min H,L m(H, L), 20cm ≤ H ≤ 40cm, 12cm ≤ L ≤ 24cm, H -2L 3 > 0, Õj .i c = 3 A. (15)
All the optimizations are executed with the corrected coarse model. The FEM is evaluated once per iteration to compute the corrector. The algorithm stops when the difference between two iterations is small enough. Optimizations use sequential quadratic programming (SQP) algorithm. This algorithm needs an initial point to start with. For the first iteration, a random initial point is used. For the following iterations, the solution of the previous iteration is used. Five OSM procedures are performed, with five random initial points, and the best solution is kept. The algorithm used is the following:

1. random point x 0 2. while Error between two iterates (a) Evolution of the error between two iterations vs space-mapping iterations. We compare in Table I the best solution obtained by using: 1. The OSM strategy. 2. An optimization with the WRM model. 3. An optimization with the analytical model. Reference solution is the solution obtained by the optimization with the WRM model. The OSM process stops at iteration j when the criterion

x j -x j-1 x j-1 > ε 2.1. Õj = i f (x j-1 ) i c (x j-1 ) 2.2. x j = arg min x m(x) s. t. Õj .i c (x) = 3 2.4. j = j + 1
x j -x j-1 x j-1 is inferior to a given tolerance ε. The OSM algorithm converges quickly to a solution which minimize the transformer mass: on the five trials, three iterations maximum are enough to obtain the optimum (Fig. 4(a)), so FEM is evaluated three times, and the objective function decreases at each iteration (Fig. 4(b)). The computation time of the optimization process is considerably reduced due to the few evaluations of the FEM, but the solution is close to the reference solution. Compared to the reference solution, the error on the objective function is 13.77% with the analytical model but 0.13% with OSM.

V. CONCLUSION

A multirate system is composed of components with very different time constants that are evaluated many times during the optimization process. In order to reduce optimization time, an OSM strategy is applied to solve the problem. Two models of the system are chosen: an analytical one and a FEM. FEM is computed using the WRM, which allows to simulate each subsystem with respect to its own time constant and guarantees the consistency of the result with a reasonable computation time. The joint action of OSM and WRM allows to obtain a solution as accurate as WRM but in a shorter time.

Fig. 1 .

 1 Fig.1. The coarse and fine models of the device.

Fig. 2 .

 2 Fig. 2. Mesh of the eighth of transformer used for the FEM.

3 Fig. 3 .

 33 Fig. 3. Transformer's geometry.

Fig. 4 .

 4 Fig. 4. Evolution of the error and of the mass vs space mapping iterations.


T Henneron

A Benabou

S Clénet

Non Linear Proper Generalized Decomposition method applied to the magnetic simulation of a SMC microstructure

Keywords:

Improvement of the magnetic performances of Soft Magnetic Composites (SMC) materials requires to link the microstructures to the macroscopic magnetic behavior law. This can be achieved with the FE method using the geometry reconstruction from images of the microstructure. Nevertheless, it can lead to large computational times. In that context, the Proper Generalized Decomposition (PGD), that is an approximation method originally developed in mechanics, and based on a finite sum of separable functions, can be of interest in our case. In this work, we propose to apply the PGD method to the SMC microstructure magnetic simulation. A non-linear magnetostatic problem with the scalar potential formulation is then solved.

I. INTRODUCTION

oft Magnetic Composites (SMC) are magnetic materials offering an interesting alternative to the use of classical laminated iron steel in electromagnetic energy conversion [START_REF] Cyr | Design and Optimization of Soft Magnetic Composite Machines With Finite Element Methods[END_REF]. Nevertheless, their magnetic and mechanical performances still remain below those of laminated magnetic steel. In order to improve the SMC performances, it is necessary to investigate their properties at the microstructure scale and to link them with the macroscopic behavior that is of interest for the designer of electromagnetic devices. Previous investigations have been realized [START_REF] Cyr | Methodology to Study the Influence of the Microscopic Structure of Soft Magnetic Composites on Their Global Magnetization Curve[END_REF] in order to reconstruct the macroscopic magnetic behavior law from the microstructure geometry. Interesting results have been obtained by the use of a classical 2D Finite Element (FE) calculation approach based on the geometry reconstruction using imaging obtained from the Scanning Electron Microscope (SEM). Nevertheless, the procedure is strongly linked to a rather delicate image processing step. In fact, the contours of the particles are interpolated between the pixels of the image. This step can introduce some error on the geometry model. It is then interesting to investigate other numerical techniques such as the Proper Generalized Decomposition (PGD) method [START_REF] Chinesta | Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models[END_REF][START_REF] Henneron | Proper Generalized Decomposition Method to Solve Quasi Static Field Problems[END_REF][START_REF] Lamari | Non-conventional numerical strategies in the advanced simulation of materials and processes[END_REF]. This approach allows representing the solution separately on each geometric axis with functions approximated by a 1D FE discretization [START_REF] Chinesta | Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models[END_REF][START_REF] Lamari | Non-conventional numerical strategies in the advanced simulation of materials and processes[END_REF]. The main interest is that we have to solve several 1D FE systems instead of a 2D FE system, but it also avoids the geometry reconstruction as required in [START_REF] Cyr | Methodology to Study the Influence of the Microscopic Structure of Soft Magnetic Composites on Their Global Magnetization Curve[END_REF] for the 2D FE approach.

In this work, the PGD approach is applied to a simple, but representative, image of a SMC microstructure. The equivalent non-linear magnetic behavior of the studied microstructure is then calculated. In a first part, the problematic and the magnetostatic problem are briefly presented. In the second part, the non-linear PGD approach is developed. Finally, the results obtained from the PGD are presented and compared with a reference 2D FE model.

II. POSITIONING OF THE PROBLEM

A. Problematic

Whatever the chosen numerical approach, it is first necessary to perform the microstructure imaging in order to identify the distribution of the different materials. For the magnetic calculation, the porosities and the resin, that is used to agglomerate the iron particles, are considered with the air properties whereas the particles are considered with pure iron magnetic properties.

In the classical FE approach, it is necessary to extract the geometry in order to build a numerical model based on a mesh. This step is rather delicate as it requires an image processing procedure that must reflects with accuracy the real particle shapes and air gaps thicknesses. In fact, at the microstructure level, the material is strongly inhomogeneous and simulations have proven that the way the image processing is carried out has an impact on the magnetic behavior that is simulated. An alternative is to treat the image using each pixel as an element of the mesh (quadrangle elements in 2D). Nevertheless, in both cases, the bigger the image is, the higher the memory and computational cost will be. In this work, we focus on the accuracy of the method regarding the global behavior representation of a SMC microstructure simulation. Figure 1 shows the example of the microstructure that is actually studied.

B. Mathematical model

Let us consider a domain D of boundary Γ (Γ=Γ B ∪Γ H1 ∪Γ H2 and Γ B ∩Γ H1 ∩Γ H2 =0) (Fig. 2). We denote by ε the magnetomotive force between Γ H1 and Γ H2 and Φ the magnetic flux flowing through the domain. The magnetostatic problem can be described by the following Maxwell's equations,

curl H = 0 div B = 0 B = µ (H) H (1) (2)
(3) with B the magnetic flux density, H the magnetic field, µ (H) the non-linear magnetic permeability that depends on the magnetic field. The unicity of the solution is obtained by boundary conditions such that, B.n=0 on Γ B and H×n=0 on Γ H1 and Γ H2 [START_REF] Henneron | Proper Generalized Decomposition Method to Solve Quasi Static Field Problems[END_REF] with n the outward unit normal vector. To solve the previous problem, the magnetic scalar formulation can be used by defining a magnetic scalar potential Ω in the whole domain.

The magnetomotive force is introduced by a scalar function α [START_REF] Dular | Dual Finite Element Formulations for Lumped Reluctances Coupling[END_REF] and the magnetic field can be expressed such that, H = -grad Ωε grad α with Ω = cst on Γ H1 and Γ H2 and α = 1 on Γ H1 and α = 0 on Γ -Γ H1 [START_REF] Lamari | Non-conventional numerical strategies in the advanced simulation of materials and processes[END_REF] Combining ( 3) and ( 5) in relation ( 2), we obtain the magnetic scalar potential formulation of the problem. The weak form to be solved is then written,

dD ' α. ε dD ' . D (H) D (H) ∫ ∫ - = grad grad grad grad (6)
with Ω' a test function defined in the same space as Ω.

III. PROPER GENERALIZED DECOMPOSITION

A. Separated representation

In order to solve (6), a method based on the PGD approach can be used [START_REF] Chinesta | Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models[END_REF][START_REF] Henneron | Proper Generalized Decomposition Method to Solve Quasi Static Field Problems[END_REF][START_REF] Lamari | Non-conventional numerical strategies in the advanced simulation of materials and processes[END_REF]. Then, the scalar potential is approximated by a separated representation in the spatial dimensions x and y:

∑ = ⋅ ≈ Ω M 1 n n n (y) S (x) R y) (x, (7)
with x∈D x , y∈D y and M the number of modes. To approximate R n (x) and S n (y), functional spaces of finite dimension are introduced, i.e. 1D nodal shape function space.

To compute the functions R n (x) and S n (y), an iterative enrichment method is used. The couple (R n (x), S n (y)) is calculated regarding the previous couples (R i (x), S i (y)) with i∈[1,n-1]. The function Ω' (see [START_REF] Dular | Dual Finite Element Formulations for Lumped Reluctances Coupling[END_REF]) can be written such that:

(y)' S (x) R (y) S (x)' R y) (x, ' n n n n ⋅ + ⋅ = Ω (8)
with R n (x)' and S n (y)' the test functions defined in the same spaces of R n (x) and S n (y) respectively.

The scalar function α used to impose the magnetomotive force is also expressed by a separated representation such that:

(y) (x).α α y) α(x, y x = (9)
The 1D scalar functions α x (x) and α y (y) are discretized in the same way as R n (x) and S n (y). The magnetic behavior law of the microstructure is represented by a matrix M µ which terms correspond to the magnetic permeability associated to each pixel of the image (either iron or air). The magnetic behavior law is also expressed by a separated representation that uses functions defined on each axis. These ones are deduced from a Singular Value Decomposition (SVD) method [START_REF] Lamari | Non-conventional numerical strategies in the advanced simulation of materials and processes[END_REF], i.e. the magnetic permeability matrix is written such that,

∑ = = = N 1 j t j jj j t V D U UDV M (10)
with D the diagonal matrix of the singular values, U and V the matrixes of the left and right singular vectors, N the number of non zero singular values and "j" stands for the j th vector. The number of lines of U and V depend on the spatial discretization of the image, i.e. the number of pixels np x and np y along x and y axes. From this decomposition, it is possible to express the magnetic permeability µ(x,y) for a given position (x,y) in the domain such that,

∑ ∑ ∑ = = = = = = y x np 1 e' je' e' j np 1 e je e j N 1 j j j jj (y).V w (y) v and (x).U w (x) u with (y) (x)v u D y) (x, (11)
where w px (x), resp. w py (y), is the scalar interpolation function associated to the pixel p x , resp. p y . This function is equal to 1 if x belongs to the pixel position and 0 otherwise. The permeability in (11) corresponds to a given magnetic state.

B. Determination of (R n (x),S n (y))

Each couple (R n (x),S n (y)) is calculated by solving iteratively two equations determined from [START_REF] Dular | Dual Finite Element Formulations for Lumped Reluctances Coupling[END_REF]. First, we suppose that S n (y) is known. Then, the function S n (y)' vanishes in (8) and the test function Ω' is equal to R n (x)'.S n (y). Equation ( 6) is then solved in order to determine the function R n (x) using the following expression, In a second step, we compute the function S n (y) assuming that the function R n (x) is known. In this case, the function R n (x)' vanishes in (8) and the test function Ω' is equal to R n (x).S n (y)'. To determine S n (y), the relation ( 6) is solved such that, The weak forms (12) and ( 13) are solved successively, using an iterative procedure, until convergence of the couple (R n (x),S n (y)). A criterion of convergence is defined by using two successive iterations. If we denote by (R n (x) j ,S n (y) j ) and (R n (x) j-1 ,S n (y) j-1 ) the couples of functions obtained for the iterations j and j-1, the error ε can be expressed in the following form,

    - - -     - =     +     ∑ ∫ ∫ ∑ ∫ ∫ ∫ ∫ ∑ ∫ ∫ ∫ ∫ ∫ ∫ ∑ - = - = = = 1 n 1 D D x ' n i j y n i j 1 n 1 D D x ' i j y n i j D D x ' n x
    - - -     - =     ∑ ∫ ∫ ∑ ∫ ∫ ∫ ∫ ∫ ∫ ∑ ∫ ∫ ∫ ∫ ∑ - = - = = = 1 n 1 i D D y ' i j x n i j 1 n 1 i D D y ' i j x n i
2 1 - j n 1 - j n j n j n ) (y) S ( ) x ( R ) (y) S ( ) x ( R t t ⋅ - ⋅ = ε (14)
with || || 2 the L 2 -norm. The stop criterion is defined when the error lower than a given tolerance.

C. Non linear problem

In order to determine the equivalent magnetic behavior law of the SMC microstructure, the non-linear behavior of the iron must be introduced in the numerical model. Then, to solve the non-linear problem, a fixed point method is used to compute Ω n q the scalar potential approximated with n modes and q denotes the non-linear iteration for a given source term ε. For each non-linear iteration q, the magnetic permeability associated to each pixel (11) of the 2D image is computed from Ω n q-1 . The SVD method (10) is applied with the new distribution of the magnetic permeability of the SMC microstructure. Ω n q is computed with the help of the relations (12) and ( 13). These steps are repeated until convergence defined by the condition ε nl lower than the tolerance. The error is calculated such that,

2 1 - q n q n nl ε Ω - Ω = (15) IV. APPLICATION

A. Presentation of the studied microstructure

We consider the microstructure given in Fig. 1. The resolution of the image is 101x101 pixels. For the discretization of the functions associated with each axis, the native resolution is kept. Then, for each axis, the 1D mesh is composed by 101 elements. The source term is introduced by imposing the magnetomotive force ε between the opposite boundaries of the system (bottom and top of the image). The behavior of the iron particles is modeled by the non linear law given in Fig. 3.

B. Influence of the number of modes

In order to evaluate the PGD approach, we compare the magnetic energy obtained from this method with the reference 2D FE approach. In this approach, the 2D mesh is directly obtained from the pixels of the image that is composed of 10201 quadrangles. In Fig. 4, the evolution of the magnetic energy W versus the number of modes is presented for a given magnetomotive force. The 2D reference result is presented in dashed line. One can observe that, with the PGD approach, the magnetic energy converges towards the one of the reference model with increasing the number of modes. Therefore, the stop criterion can be obtained from the evaluation of the energy variation between two successive results.

In term of computation time, in that case, the PGD approach is more time consuming than the 2D FE approach. In fact, the studied system is rather simple and has a small number of unknowns. The PGD is penalized by the important number of modes that is required to obtain acceptable results. Nevertheless, in term of memory resources, the PGD solves two 1D FE systems whereas the reference solves a 2D FE system. For high resolution images, the PGD method will be more interesting as its number of unknowns increases linearly. At the opposite, the 2D FE model has a number of unknowns that increases with a quadratic law.

C. Determination of the equivalent magnetic behavior

In order to determine the equivalent behavior law µ eff (H) of the microstructure image, we impose different values of the magnetomotive force. For each calculation, the magnetic coenergy Wc(ε) is computed. Then, as the global magnetic permeance P(ε) can be expressed in function of the derivative of Wc(ε) with respect to ε, the effective magnetic permeability µ eff (H) can be written such that µ eff (H) =(L/S)P(ε) with H=ε/L where L is the length of the microstructure and S the section that is seen by the magnetic flux. In Fig. 5, the magnetic energy and coenergy are presented versus H. The equivalent magnetic behavior law can then be deduced as shown in Fig. 6. Due to air-gaps and porosities between the iron particles, the equivalent magnetic behavior law is as expected beneath the one of pure iron. Due to the inhomogeneity, the behavior law is rather uncommon. In fact, the reconstruction of the macroscopic behavior law of a SMC will require larger micro-structural images, typically with more than 15 iron particles as observed in [START_REF] Cyr | Methodology to Study the Influence of the Microscopic Structure of Soft Magnetic Composites on Their Global Magnetization Curve[END_REF]. As illustration of magnetic states for two levels of H, maps of the effective magnetic permeability are given in Fig. 7. The darker the pixel is, the lower the magnetic permeability is (black is µ 0 ). One can observe that for low level of H, the local saturation occurs on the left of the image. For the higher level of H, multiple local saturations appear. The non-linear PGD approach with scalar potential formulation has been developed and applied to the study of a SMC microstructure. To illustrate the feasibility to use this method in such case, a simple image of the microstructure is used. The PGD approach has been compared with the 2D FE and the obtained results are in good agreement. Nevertheless, the reconstruction of the macroscopic behavior law of a SMC requires several larger micro-structural images [START_REF] Cyr | Methodology to Study the Influence of the Microscopic Structure of Soft Magnetic Composites on Their Global Magnetization Curve[END_REF]. In that case, the PGD method could be more interesting in terms of computational effort.

Fig. 1 .

 1 Fig. 1. Studied microstructure obtained from SEM

Fig. 2 .

 2 Fig.2. Magnetostatic problem

Fig. 3 .

 3 Fig.3. Non linear behavior law of the iron

Fig. 4 .

 4 Fig.4. Magnetic energy versus the number of modes (dashed line is the 2D reference)

Fig. 5 .

 5 Fig.5. Magnetic energy and coenergy

Fig. 6 .

 6 Fig.6. µgl(H)


S L Gavrilyuk
email: sergey.gavrilyuk@univ-amu.fr

M Yu Kazakova
email: m.u.kazakova@gmail.com

Hydraulic jumps in two-layer flows with a free surface

A closure relation is proposed for description of hydraulic jumps in twolayer flows with a free surface over a horizontal bottom. Such a relation comes from the momentum equations of each layer which become in a sense conservative when the total momentum and the masses in each layer are conserved. It is also shown that this relation guarantees that the energy flux decreases through the jump.

Governing equations

Consider a flow of two immiscible heavy fluid layers of depths h 1 , h 2 over a horizontal bottom (see Figure 1). In the long wave approximation the following model can be derived [START_REF] Ovsyannikov | Nonlinear problems in the theory of surface and internal waves[END_REF][START_REF] Baines | Topographic effects in stratified flows[END_REF]:

(γ 1 h 1 ) t + (γ 1 h 1 u 1 ) x = 0, (1)
(γ 2 h 2 ) t + (γ 2 h 2 u 2 ) x = 0, (2)
γ 1 h 1 (u 1t + u 1 u 1x ) + gh 1 (γ 1 h 1 + γ 2 h 2 ) x = 0, (3)
γ 2 h 2 (u 2t + u 2 u 2x ) + gh 2 γ 2 (h 1 + h 2 ) x = 0. (4)
Here γ 1 , γ 2 (γ 1 > γ 2 ) are constant densities, and u 1 , u 2 are velocities in each layer. The gravity acceleration is denoted by g. It was recently shown [START_REF] Montgomery | On the number of conserved quantities for the two-layer shallow water equations[END_REF]Moodie, 2001, Barros, 2006) that, in addition to the mass conservation laws, the only conservation laws admitted by the system are the conservation of the total momentum

(γ 1 h 1 u 1 + γ 2 h 2 u 2 ) t + (γ 1 h 1 u 2 1 + γ 2 h 2 u 2 2 + E) x = 0, (5)
γ 1 h 1 u 2 1 2 + γ 2 h 2 u 2 2 2 + E t (6)
+ u 1 h 1 γ 1 u 2 1 2 + g(γ 1 h 1 + γ 2 h 2 ) + u 2 h 2 γ 2 u 2 2 2 + gγ 2 (h 1 + h 2 ) x = 0,
and the Bernoulli integrals

γ 1 u 1t + γ 1 u 2 1 2 + g (γ 1 h 1 + γ 2 h 2 ) x = 0, (7)
γ 2 u 2t + γ 2 u 2 2 2 + gγ 2 (h 1 + h 2 ) x = 0. ( 8
)
Here the "internal energy" E is defined as

E = gγ 1 h 2 1 /2 + gγ 2 h 1 h 2 + gγ 2 h 2 2 /2.
The corresponding characteristic polynomial can easily be derived from (1) -( 4) :

(u 1 -λ) 2 -gh 1 (u 2 -λ) 2 -gh 2 -g 2 γ 2 γ 1 h 1 h 2 = 0. ( 9
)
The system is conditionally hyperbolic provided the relative velocity w = u 1 -u 2 is small [START_REF] Ovsyannikov | Two-layer "shallow water" model[END_REF][START_REF] Baines | Topographic effects in stratified flows[END_REF][START_REF] Barros | Conservation laws for one-dimensional shallow water models for one and two-layer flows[END_REF][START_REF] Abgrall | Two-layer shallow water systems : a relaxation approach[END_REF]. The hyperbolicity implies a possibility of shock formation (corresponding to hydraulic jumps) even if the initial data are smooth.

A dispersive regularization of (1), ( 2), ( 3) and (4) has been derived in [START_REF] Liska | Nonhydrostatic two-layer models of incompressible flow[END_REF] and studied in [START_REF] Barros | Dispersive Nonlinear Waves in Two-Layer Flows with Free Surface[END_REF] and Barros and Gavrilyuk (2007). In particular, solitary wave solutions of the dispersive equations have been studied in the last two references. For the case where the velocities at infinity in each layer are equal, it has been found that dispersive equations admit "table-top" (and "table-down") internal solitary waves for some singular values of Froude numbers. The amplitude of the free surface waves is, in general, much smaller than the amplitude of internal waves. The half of such a "table-top" soliton can be considered as a regular "hydraulic jump" where the total momentum and the total energy are conserved.

The singular hydraulic jumps considered for the two -layer shallow water system (1), ( 2), ( 3) and ( 4) are a symplification of ondular bores having nonstationary oscillating tail after the first smooth front representing a solitary wave. Such ondular bores have been intensively studied in one-layer shallow water flows [START_REF] El | Unsteady undular bores in fully nonlinear shallow-water theory[END_REF]Smyth, 2006, LeMetayer, Gavrilyuk and[START_REF] Metayer | A numerical scheme for the Green-Naghdi model[END_REF], however their nonstationary analysis is always absent for the two-layer case.

An old question concerning singular hydraulic jumps in two-layer flows is : what are the Rankine-Hugoniot relations determining such a jump? The mass conservation laws (1), ( 2) and the conservation of the total momentum (5) are obvious candidates. The energy conservation law plays the role of the entropy inequality through the jumps (see the corresponding discussion in the case of one-layer flows in [START_REF] Stoker | Water Waves[END_REF]. Finally, one of the Bernoulli conservation laws can be used, but which one? The same problem appears in two-layer flows between rigid lids (see Baines, 1995, §3.5) where one of the possible closing hypotheses could be to privilege the Bernoulli equation in the contracting layer. However, as noticed in [START_REF] Baines | Topographic effects in stratified flows[END_REF]this assumption can not be strictly correct, as observations show that there is some dissipation occurring in each layer". The use of the Bernoulli laws for the two-layer flows with free surface is thus also questionable.

Studying shear flows in homogeneous fluids, [START_REF] Liapidevskii | Mathematical models of propagation of long waves in non-homogeneous fluids[END_REF] proposed to use the local vorticity to closure the jump relations. In the case of two-layer flows an analogue of such a vorticity is the velocity difference. [START_REF] Ostapenko | Stable shock waves in two layer "shallow water[END_REF] has been proved that if we choose the layer depths h 1 , h 2 , the total momentum γ 1 h 1 u 1 + γ 2 h 2 u 2 and the velocity difference u 1 -u 2 as the basic variables, the domain of convexity of the total energy as a function of such variables is maximal compared with other choices of conservative variables (for example, h 1 , h 2 , u 1 , u 2 ). He considered it as a mathematical argument justifying the choice of the conservation law for the relative velocity. Also, in the case of very weak stratification (γ 1 ≈ γ 2 ), in [START_REF] Liapidevskii | Mathematical models of propagation of long waves in non-homogeneous fluids[END_REF] the choice of the conservation law for the relative velocity has been used to solve the Riemann problem for two-layer flows under a rigid lid. In particular, such a choice allowed them to construct a solution of the Riemann problem in the hyperbolicity region satisfying the Lax stability condition on the hydraulic jumps.

Abgrall and Karni (2009) used the rigid-lid assumption at shocks to study numerically the two-layer flows over topography. The same problem with jump conditions resulting from the classical choice of a linear path through the space state was studied in [START_REF] Mandli | Finite volume methods for the multilayer shallow water equations with applications to storm surges[END_REF]. The jump relations in the last two cases are not, a priori, compatible with the entropy inequality. However, this inequality was always checked a posteriori.

In this paper, we exploit another idea proposed by V. Teshukov in 2007. He remarked that the non-conservative momentum equations for each layer become conservative at the variety where the mass of each layer and the total momentum are conserved. In particular, we will show that such a closure is compatible with the entropy inequality.

Basic relations

We suppose the hydraulic jump is stationary in the reference system moving with the jump. The mass and the total momentum equations ( 1), ( 2), ( 5) can then be integrated :

γ 1 h 1 u 1 = Q 1 = const, (10)
γ 2 h 2 u 2 = Q 2 = const, (11)
γ 1 h 1 u 2 1 + γ 2 h 2 u 2 2 + gγ 1 h 2 1 /2 + gγ 2 h 1 h 2 + gγ 2 h 2 2 /2 = M = const. ( 12
)
Eliminating the velocities, the momentum equation becomes

Q 2 1 γ 1 h 1 + Q 2 2 γ 2 h 2 + gγ 1 h 2 1 /2 + gγ 2 h 1 h 2 + gγ 2 h 2 2 /2 = M = const. (13)
For given γ 1 , γ 2 (γ 1 > γ 2 ) and small Q 1 and Q 2 , the function

M (h 1 , h 2 ) = Q 2 1 γ 1 h 1 + Q 2 2 γ 2 h 2 + gγ 1 h 2 1 /2 + gγ 2 h 1 h 2 + gγ 2 h 2 2 /2 (14)
has only one non-degenerate critical point which is the minimum point. Indeed, consider the function M (h 1 , h 2 ) defined by ( 14). Its Hessian matrix (denoted by M ′′ (h 1 , h 2 )) is positive definite

M ′′ (h 1 , h 2 ) =   2Q 2 1 γ1h 3 1 + gγ 1 gγ 2 gγ 2 2Q 2 2 γ2h 3 2 + gγ 2   > 0 because γ 1 > γ 2 . The minimum point of M (h 1 , h 2
) is defined from the equations :

h 2 1 h 1 + γ 2 γ 1 h 2 = Q 2 1 gγ 2 1 , (15)
h 2 2 (h 2 + h 1 ) = Q 2 2 gγ 2 2 .
The level curves of M (h 1 , h 2 ) in the vicinity of the critical point are then closed. Consider the total energy

E = γ 1 h 1 u 2 1 2 + γ 2 h 2 u 2 2 2 + gγ 1 h 2 1 /2 + gγ 2 h 1 h 2 + gγ 2 h 2 2 /2.
The partial derivatives of the energy with respect to velocities u i and depths h i define the fluxes of the Bernoulli integrals ( 7), ( 8) and the mass fluxes(1), ( 2) [START_REF] Bridges | Criticality manifolds and their role in the generation of solitary waves for two-layer flow with a free surface[END_REF] :

∂E ∂h 1 = γ 1 u 2 1 2 + g(γ 1 h 1 + γ 2 h 2 ) = R 1 , ∂E ∂u 1 = γ 1 h 1 u 1 = Q 1 , ∂E ∂h 2 = γ 2 u 2 2 2 + gγ 2 (h 1 + h 2 ) = R 2 , ∂E ∂u 2 = γ 2 h 2 u 2 = Q 2 .
The energy flux is given by

D (h 1 , h 2 ) = Q 1 γ 1 R 1 + Q 2 γ 2 R 2 = Q 1 γ 1 Q 2 1 2γ 1 h 2 1 + g(γ 1 h 1 + γ 2 h 2 ) + Q 2 γ 2 Q 2 2 2γ 2 h 2 2 + gγ 2 (h 1 + h 2 ) . ( 16
)
The Hessian of D is :

D ′′ (h 1 , h 2 ) =   3Q 3 1 γ 2 1 h 4 1 0 0 3Q 3 2 γ 2 2 h 4 2   .
The extremum point of D is given by :

Q 1 γ 1 - Q 2 1 γ 1 h 3 1 + gγ 1 + Q 2 γ 2 gγ 2 = 0, Q 1 γ 1 gγ 2 + Q 2 γ 2 - Q 2 2 γ 2 h 3 2 + gγ 2 = 0.
Or

h 3 1 = Q 3 1 gγ 2 1 (Q 1 + Q 2 ) , ( 17
)
h 3 2 = Q 3 2 gγ 2 2 γ2 γ1 Q 1 + Q 2 .
Two cases should be distinguished.

• the flow in both layers is in the same direction, i.e. the signs of the mass fluxes are the same. For definiteness, we suppose that the fluxes are positive, i.e. the flow is from the left to the right :

Q 1 > 0, Q 2 > 0. ( 18
)
• The signs of the mass fluxes are opposite. One can choose, for example,

Q 1 > 0, Q 2 < 0. ( 19
)
The case ( 19) is more difficult because some additional singularities will appear in (17). In the following we will consider only the case (18). In this case, the extremum point is a minimum point, and the level curves D (h 1 , h 2 ) are also closed.

At the minimum points of M (h 1 , h 2 ) and D (h 1 , h 2 )( denoted by M M and D M , respectively), the velocities of each layer u 1 , u 2 are related to the corresponding depths h 1 , h 2 by ( 15) and ( 17). One can prove that the points M M and D M always belong to the hyperbolicity region (the correponding characteristic polynomial (9) has three positive and one zero root at these points). The fact that one of the roots is always zero, is a very important property which is worth to be discussed in details. Consider the critical curve :

C (h 1 , h 2 ) = Q 2 1 γ 2 1 h 2 1 -gh 1 Q 2 2 γ 2 2 h 2 2 -gh 2 -g 2 γ 2 γ 1 h 1 h 2 = 0. ( 20
)
obtained from the characteristic polynomial ( 9) by taking λ = 0. Theorem 1.If (h 1 , h 2 ) are extrema of M (h 1 , h 2 ) and D (h 1 , h 2 ) defined, respectivement, by ( 15) and ( 17) then C = 0, i.e. they belong to the critical curve.

Proof. The minimum point of M is defined by ( 15) and satisfies the relations

Q 2 1 γ 2 1 h 2 1 = g h 1 + γ 2 γ 1 h 2 , Q 2 2 γ 2 2 h 2 2 = g (h 2 + h 1 ) . Hence C = Q 2 1 γ 2 1 h 2 1 -gh 1 Q 2 2 γ 2 2 h 2 2 -gh 2 -g 2 γ 2 γ 1 h 1 h 2 = g h 1 + γ 2 γ 1 h 2 -gh 1 (g (h 2 + h 1 ) -gh 2 ) -g 2 γ 2 γ 1 h 1 h 2 = 0.
The minimum point of D is defined by ( 17) and satisfies the relations

Q 2 1 γ 2 1 h 2 1 = g (Q 1 + Q 2 ) Q 1 h 1 , Q 2 2 γ 2 2 h 2 2 = g γ2 γ1 Q 1 + Q 2 Q 2 h 2 . Hence C = Q 2 1 γ 2 1 h 2 1 -gh 1 Q 2 2 γ 2 2 h 2 2 -gh 2 -g 2 γ 2 γ 1 h 1 h 2 = g (Q 1 + Q 2 ) Q 1 h 1 -gh 1   g γ2 γ1 Q 1 + Q 2 Q 2 h 2 -gh 2   -g 2 γ 2 γ 1 h 1 h 2 = 0.
The theorem is proved.

Remark 1.The line of equal velocities u 1 = u 2 which can be written as :

Q 1 γ 1 h 1 - Q 2 γ 2 h 2 = 0
is always above the minimum points of M (h 1 , h 2 ) and D(h 1 , h 2 )in the plane (h 1 , h 2 ). Indeed, consider the minimum point of M (h 1 , h 2 ). We obtain

Q 2 1 γ 2 1 h 2 1 = u 2 1 = g h 1 + γ 2 γ 1 h 2 , Q 2 2 γ 2 2 h 2 2 = u 2 2 = g (h 2 + h 1 ) . Hence Q 2 1 γ 2 1 h 2 1 - Q 2 2 γ 2 2 h 2 2 = u 2 1 -u 2 2 = gh 2 γ 2 γ 1 -1 < 0, i.e. u 1 -u 2 < 0 at the minimum point of M (h 1 , h 2 ).
Consider now the minimum point of D(h 1 , h 2 ). We get

Q 3 1 γ 3 1 h 3 1 = g (Q 1 + Q 2 ) γ 1 , Q 3 2 γ 3 2 h 3 2 = g γ2 γ1 Q 1 + Q 2 γ 2 . Hence Q 3 1 γ 3 1 h 3 1 - Q 3 2 γ 3 2 h 3 2 = u 3 1 -u 3 2 = g (Q 1 + Q 2 ) γ 1 - g γ2 γ1 Q 1 + Q 2 γ 2 = g Q 2 γ 2 γ 2 γ 1 -1 < 0,
i.e. u 1 -u 2 < 0 at the minimum point of D(h 1 , h 2 ). One can also show that the point D M is situated above the point M M . In Figure 2, we show a typical relative position of the curves C =0, M = const, D = const, u 1 = u 2 and the minimum points D M and M M at the plane (h 1 , h 2 ). Theorem 2. The curves M (h 1 , h 2 ) = const and D(h 1 , h 2 ) = const are tangent at a point of the plane (h 1 , h 2 )if and only if at this point u 1 = u 2 or C =0.

Proof. Consider the curves M (h 1 , h 2 ) = const and D(h 1 , h 2 ) = const defined by ( 14), ( 16)

M (h 1 , h 2 ) = Q 2 1 γ 1 h 1 + Q 2 2 γ 2 h 2 + gγ 1 h 2 1 /2 + gγ 2 h 1 h 2 + gγ 2 h 2 2 /2, D(h 1 , h 2 ) = Q 1 γ 1 Q 2 1 2γ 1 h 2 1 + g(γ 1 h 1 + γ 2 h 2 ) + Q 2 γ 2 Q 2 2 2γ 2 h 2 2 + gγ 2 (h 1 + h 2 ) .
The curves are tangent if and only if

0 = ∂M ∂h 1 ∂D ∂h 2 - ∂D ∂h 1 ∂M ∂h 2 = - Q 2 1 γ 1 h 2 1 + gγ 1 h 1 + gγ 2 h 2 g Q 2 + Q 1 γ 2 γ 1 - Q 3 2 γ 2 2 h 3 2 - M = const С 0 = h 1 h 2 D = const M M D M

8

--

Q 2 2 γ 2 h 2 2 + gγ 2 h 1 + gγ 2 h 2 g (Q 2 + Q 1 ) - Q 3 1 γ 2 1 h 3 1 = = - (Q 2 γ 1 h 1 -Q 1 γ 2 h 2 ) h 3 1 h 3 2 γ 2 1 γ 2 2 Q 2 1 -Q 2 2 + gh 3 2 γ 2 2 + gh 3 1 γ 1 Q 2 2 γ 1 + gh 3 2 γ 2 2 (γ 2 -γ 1 ) = = g 2 (Q 2 γ 1 h 1 -Q 1 γ 2 h 2 ) Q 2 2 gh 3 2 γ 2 2 -1 Q 2 1 gh 3 1 γ 2 1 -1 - γ 2 γ 1
Hence, the Jacobian matrix is degenerate if and only if

Q 2 γ 1 h 1 -Q 1 γ 2 h 2 = 0 or Q 2 2 gh 3 2 γ 2 2 -1 Q 2 1 gh 3 1 γ 2 1 -1 - γ 2 γ 1 = 0.
This is equivalent to relations u 2 = u 1 or C = 0.

Supercritical-subcritical transition

Consider a discontinuous piecewise constant stationary solultion of ( 1) -( 4)

(h 1 (x) , h 2 (x) , u 1 (x) , u 2 (x)) T = (h 10 , h 20 , u 10 , u 20 ) T , if x < 0 (h 1 , h 2 , u 1 , u 2 ) T , if x > 0
For convenience, we will use the same notations for the right state constant variables as for the basic variables (without additional indices). The critical curve C (h 1 , h 2 ) = 0 is defined by ( 20). The flow is supercritical if C > 0 and subcritical if C < 0. Supppose that the flow satisfies Q 1 > 0, Q 2 > 0 and the state "0" at the left is supercritical. The constant states satisfy the relations (10), ( 11) and ( 12). An additional relation is needed satisfying the natural dissipation inequality (very often called in the context of hyperbolic equations "entropy inequality") :

[D (h 1 , h 2 )] = Q 1 γ 1 R 1 + Q 2 γ 2 R 2 = Q 1 γ 1 Q 2 1 2γ 1 h 2 1 + g(γ 1 h 1 + γ 2 h 2 ) + Q 2 γ 2 Q 2 2 2γ 2 h 2 2 + gγ 2 (h 1 + h 2 ) ≤ 0.
where

[D] = D (h 1 , h 2 ) -D (h 10 , h 20 )
Since the curve D (h 1 , h 2 ) = D (h 10 , h 20 ) is closed, the "entropy inequality" means that the subcritical state should be inside this closed curve. In particular, the entropy inequality must imply that the state at the right is subcritical.

For given positive Q 1 , Q 2 , we choose (h 10 , h 20 ) such that the whole closed curve M (h 1 , h 2 ) = M (h 10 , h 20 ) belongs to the hyperbolicity domain (see Figure 3). This allows us to avoid additional mathematical difficulties relating to the interpretation of the solution in the elliptic domain.

The curve M (h 1 , h 2 ) = M (h 10 , h 20 ) can be parametrized, if we introduce the polar coordinates (r, θ) :

h 1 = r (θ) cos θ, h 2 = r (θ) sin θ, θ ∈ θ -, θ + .
Here r (θ) is determined from the third order polynomial coming from the relation

Q 2 1 γ 1 h 1 + Q 2 2 γ 2 h 2 + gγ 1 h 2 1 /2 + gγ 2 h 1 h 2 + gγ 2 h 2 2 /2 = M =const.
It has the form

r 3 g γ 1 (cos θ) 2 /2 + γ 2 cos θ sin θ + γ 2 (sin θ) 2 /2 -M r + Q 2 1 γ 1 cos θ + Q 2 2 γ 2 sin θ = 0.
(21) One of the roots of this polynomial is always negative (hence, non-physical), the two other roots are real and positive : r 1 (θ) < r 2 (θ), for θ ∈ (θ -, θ + ). The union of these two brunches form a smooth closed curve M = const in the plane (r, θ) : 4). Consider now the stationary non-conservative momentum equations ( 1), (2) for each layer :

r 1 (θ -) = r 2 (θ -), r 1 (θ + ) = r 2 (θ + ) (see Figure
γ 1 h 1 u 1 u 1x + gh 1 (γ 1 h 1 + γ 2 h 2 ) x = 0, γ 2 h 2 u 2 u 2x + gh 2 γ 2 (h 1 + h 2 ) x = 0.
They become conservative along the curve M = const. Indeed, the nonconservative product h 1 (h 2 ) x is now :

h 1 (h 2 ) x = r(θ) cos θ dr(θ) dθ θ x sin θ + r(θ)θ x cos θ =   r 2 (θ) 2 sin θ cos θ + 1 2 θ θ * r 2 (θ)dθ   x = =   h 1 h 2 2 + 1 2 θ θ * r 2 (θ)dθ   x .
Analogously,

h 2 (h 1 ) x = r(θ) sin θ dr(θ) dθ θ x cos θ -r(θ)θ x sin θ =   r 2 (θ) 2 sin θ cos θ - 1 2 θ θ * r 2 (θ)dθ   x = =   h 1 h 2 2 - 1 2 θ θ * r 2 (θ)dθ   x . M = onst с С 0 = h 1 h 2 q + q - q + с q - с - Figure 4:
The momentum curve M = const can be parametrized in polar coordinates in the plane (h 1 , h 2 ).

Here θ * is any constant belonging to (θ -, θ + ). Finally, one of the following conserved quantities can be chosen:

Q 2 1 γ 1 h 1 +g   γ 1 h 2 1 2 + γ 2   h 1 h 2 + 1 2 θ θ * r 2 (θ)dθ     = Q 2 1 γ 1 h 10 +g   γ 1 h 2 10 2 + γ 2   h 10 h 20 2 + 1 2 θ0 θ * r 2 (θ)dθ     = M 1 , Q 2 2 γ 2 h 2 +g   γ 2 h 2 2 2 + γ 2   h 1 h 2 2 - 1 2 θ θ * r 2 (θ)dθ     = Q 2 2 γ 2 h 20 +g   γ 2 h 2 20 2 + γ 2   h 10 h 20 2 - 1 2 θ0 θ * r 2 (θ)dθ     = M 2 .
This system is compatible with the total momentum. Indeed, summing them we obtain

M 1 + M 2 = M.
We can choose the conservation of the local momentum M 1 (or M 2 ) as an additional closure relation. We have now to check if the corresponding state is in the subcritical region, and that the "entropy inequality" is satisfied. For definiteness, we choose the local momentum equation of the first layer. Obviously, this relation does not depend on the choice of θ * . Finally, for a given supercritical state h 10 , h 20 we have to find a subcritical state h 1 , h 2 from the following system of equations : 4 Solution algorithm to find the subcritical state

Q 2 1 γ 1 h 1 + gγ 1 h 2 1 2 + gγ 2 h 1 h 2 2 + gγ 2 2 θ θ0 r 2 (θ)dθ = Q 2 1 γ 1 h 10 + gγ 1 h 2 10 2 + gγ 2 h 10 h 20 2 , (22)
Q 2 1 γ 1 h 1 + Q 2 2 γ 2 h 2 +gγ 1 h 2 1 /2+gγ 2 h 1 h 2 +gγ 2 h 2 2 /2 = Q 2 1 γ 1 h 10 + Q 2 2 γ 2 h
The following particular parameters were chosen to calculate the constant M (h 10 , h 20 ):

h 10 = 0.45, h 20 = 0.045, Q 1 = 10, Q 2 = 1.1, γ 1 = 10, γ 2 = 9.5.
They determine the curve M (h 1 , h 2 ) = M (h 10 , h 20 ) entirely belonging to the hyperbolicity region. In the following, (h 10 , h 20 ) will mean any generic point at this fixed curve, and not only these chosen particular values. It is more simpler to see the solution structure in the plane (r, θ). The curve C = 0 cuts the curve M = const at points corresponding to angles θ

- c < θ + c such that r q q + q - С 0 = M = onst с с q + q - с Figure 5: A typical behaviour of the curves C = 0, M = const in (r, θ) -plane is shown.
and look for a non-trivial solution θ with r = r 2 (θ). It was found that this value of θ was unique and always satisfied the "entropy inequality". A typical picture is shown in (Figure 6 where S is an initial supercritical state, and J is the corresponding subcritical state. The following property was also established : if the initial supercritical data vary monotonically with respect to the angle θ 0 (they are denoted by points enumerated from 1 to 8 in Figure 7), the same monotonicity is observed for the subcritical states (they are denoted by squares enumerated from 1 to 8) (see Figure 7).

Comparison with other approach based on the Bernoulli integrals

An alternative approach consists in choosing the conservation law for the velocity difference as a closing jump relation : For stationary hydraulic jumps, it corresponds to the conservation of the quantity :

(u 1 -u 2 ) t + u 2 1 2 - u 2 2 2 + gh 2 γ 2 γ 1 -1 x = 0. s J С 0 = D = onst с r M = onst с q q + q - с q + q - с
u 2 1 2 - u 2 2 2 + gh 2 γ 2 γ 1 -1 = const
It can be shown that in both approaches, the "entropy inequality" is satisfied.

For ten points distributed uniformly at the subcritical part of the curve M = const and enumerated in increasing order with respect to θ, the results for the corresponding subcritical states are shown in Figures 8 and9. For both, internal and surface jumps, the amplitudes of jumps obtained by a new method are larger than those obtained by using the difference of the Bernoulli integrals.

One can also numerically prove that in the new approach the following inequality is always satisfied :

h 1 h 10 > h 1 + h 2 h 10 + h 20 (24)
i.e. the amplitude of the internal jump is always larger than the amplitude of the jump corresponding to surface waves. When the Bernoulli approach is used, the inequality ( 24) is not always satisfied (Figure 10). The inequality (24) corroborates a well-known fact concerning two-layer flows [START_REF] Armi | The hydraulics of two flowing layers with different densities[END_REF] : "two-layer flows with a free surface and small non-dimensional-density difference behave exactly as bounded two-layer flows with the upper boundary level".

Figure 1 :

 1 Figure 1: A physical picture of a two-layer flow with a free surface over a flat bottom

Figure 2 :

 2 Figure 2: A typical behaviour of the curves C = 0, M = const, D = const, and u 1 = u 2 is shown. The minimum points D M and M M are always below the line u 1 = u 2 , and M M is always below D M .

Figure 6 :

 6 Figure 6: The supercritical state S and the subcritical state J belong to the same closed solid curve M = const. The state J satisfies the "entropy inequality" : it is inside the dashed closed curve D = const passing by S.

Figure 7

 7 Figure 7: 8 initial states were uniformly distributed at the supercritical part of the curve M = const (shown by points). The corresponding subcritical states are shown by squares.

Figure 8 :

 8 Figure8: The deviation of the internal jump for an array of ten points is shown for a new approach (by dots ) and the Bernoulli approach (by empty dots).

Acknowledgments The authors acknowledge support from the Ministry of Education and Science of the Russian Federation (agreement 2012-1.5-8503) and from the President Program for Support of Leading Scientific Schools (grant No. NSh-6706.2012.1). The authors would like to thank S. V. Golovin for pertinent discussions, and R. Leveque for useful references.

The hyperbolicity domain is non-colored, while the ellipticity domain is in gray (we have there two real and two complex conjugate eigenvalues of the polynomial ( 9)). The curve M = const is chosen to be entirely in the hyperbolicity region.

N

Figure 10: For the new approach, the amplitude of the internal jump is always larger than the amplitude of the jump corresponding to surface waves (shown by triangles), while for the Bernoulli approach it is not always the case (shown by empty triangles).


D Chablat

R Jha

F Rouillier
email: fabrice.rouillier@inria.fr

G Moroz

Non-singular assembly mode changing trajectories in the workspace for the 3-RPS parallel robot

Keywords: Parallel robot, 3-RPS, Singularity, Operation mode, Aspect, Cylindrical algebraic decomposition

Having non-singular assembly modes changing trajectories for the 3-RPS parallel robot is a well-known feature. The only known solution for defining such trajectory is to encircle a cusp point in the joint space. In this paper, the aspects and the characteristic surfaces are computed for each operation mode to define the uniqueness of the domains. Thus, we can easily see in the workspace that at least three assembly modes can be reached for each operation mode. To validate this property, the mathematical analysis of the determinant of the Jacobian is done. The image of these trajectories in the joint space is depicted with the curves associated with the cusp points.

Introduction

When designing a robot, the last step is the trajectory planning. The task of the robot is generally defined in the workspace whereas the control loop depends on the joint space parameters. While defining the home pose of the robot, the Cartesian pose and the Joint values of the actuators are known. If the trajectory planning is done in the workspace by analyzing only the determinant of the Jacobian, we can reach a Cartesian pose different from the home pose but with the same joint value. This feature is called a nonsingular assembly mode changing trajectory and stands only for the parallel robot.

For such robots, the inverse and direct kinematic problem (DKP) can have several solutions. To cope up with this problem, the notion of aspects was introduced for the serial robot in [START_REF] Borrel | A study of manipulator inverse kinematic solutions with application to trajectory planning and workspace determination[END_REF] and for the parallel robot in [START_REF] Wenger Ph | Definition Sets for the Direct Kinematics of Parallel Manipulators[END_REF] and [START_REF] Chablat | Working Modes and Aspects in Fully-Parallel Manipulator[END_REF]. For the serial robots, the aspects are defined as the maximal singularity-free sets in the joint space whereas in case of parallel robots, the aspects are defined as the maximal singularity-free sets in the workspace or the cross-product of the joint space by the workspace. However, there exists robots, referred as cuspidal robots, which are able to change the inverse kinematic solution without passing through a singularity for serial robots or direct kinematic solution without passing through a singularity for parallel robots [START_REF] Innocenti | Singularity-free evolution from one configuration to another in serial and fully-parallel manipulators[END_REF][START_REF] Chablat | Séparation des solutions aux modèles géométriques direct et inverse pour les manipulateurs pleinement parallèles[END_REF][START_REF] Caro | Non-Singular Assembly Mode Changing Trajectories of a 6-DOF Parallel Robot[END_REF][START_REF] Husty | Non-Singular Assembly Mode Change of 3-RPS Manipulators[END_REF][START_REF] Macho | Transitions between Multiple Solutions of the Direct Kinematic Problem[END_REF]. The uniqueness domains are the connected subsets of the aspects induced by the characteristic surface. These notions are defined more precisely in sections 2.3 and 2.4 .

The paper elucidates the non-singular assembly mode changing trajectories in the workspace for the 3-RPS parallel robot. In Section 2.1 we describe the 3-RPS parallel robot, in section 2.2 we set the related kinematic equations while in section 2.3 we define the aspects for an operation mode. In section 2.4 we analyze the characteristic surfaces for an operation mode, and in section 2.5 we report the non-singular assembly modes changing trajectory between the two basic regions.

Kinematics

Mechanism under study

The robot under study is the 3-RPS parallel

A 1 A 2 A 3 B 1 B 2 B 3 P O g h x y z u v w ρ 1 ρ 2 ρ 3
Fig. 1 3-RPS parallel robot robot with three degrees of freedom. It been studied by many researchers [START_REF] Husty | Non-Singular Assembly Mode Change of 3-RPS Manipulators[END_REF][START_REF] Babu | Design for optimal performance of 3-RPS parallel manipulator using evolutionary algorithms[END_REF]. It is the assembly of two equilateral triangles (the base and the moving platform) by three identical RPS legs where R is a revolute passive joint, P an prismatic joint and S a passive spherical joint. Thus, the revolute joint is connected to the fixed base and the spherical joint to the mobile platform.

Considering the 3-RPS parallel manipulator, as shown in figure 1, the fixed base consists of an equilateral triangle with vertices A 1 , A 2 and A 3 , and circumradius g. The moving platform is another equilateral triangle with vertices B 1 , B 2 and B 3 , circumradius h and circumcenter P . The two design parameters g and h are positive numbers. Connecting each of the pairs of vertices of A i , B i (i = 1, 2, 3) by a limb, a rotational joint lies at A i and a spherical joint lies at B i . ρ i denotes the length of each limb and their adjustment is done through an actuated prismatic joint. Thus we get five parameters, namely g, h, ρ 1 , ρ 2 and ρ 3 . g and h are the two design parameters determine the design of the manipulator whereas the joint parameters ρ 1 , ρ 2 and ρ 3 determine the motion of the robot. To simplify the equations, we will study a unit robot with g = h = 1.

Kinematic equations

The transformation from the moving frame to the fixed frame can be described by a position vector p = OP and a 3 × 3 rotation matrix R. Let u, v and w be the three unit vectors defined along the axes of the moving frame, then the rotation matrix can be expressed in terms of the coordinates of u, v and w as:

R =   u x v x w x u y v y w y u z v z w z   (1)
The vertices of the base triangle and mobile platform triangle are

A 1 =   g 0 0   A 2 =   -g/2 g √ 3/2 0   A 3 =   -g/2 -g √ 3/2 0   (2)
b 1 =   h 0 0   b 2 =   -h/2 h √ 3/2 0   b 3 =   -h/2 -h √ 3/2 0   (3)
The coordinates of b i with respect to fixed frame reference are obtained by B i = P + Rb i for i = 1, 2, 3. Also the coordinates of the centre of the mobile platform in the fixed reference is P = [x y z] T . The distance constraints yields:

||A i -B i || = ρ 2 i with i = 1, 2, 3 (4)
As A i are revolute joints, the motion of the B i are constrained in planes. This leads to the three constraint equations:

u y h + y = 0 (5) y -u y h/2 + √ 3v y h/2 + √ 3x - √ 3u x h/2 + 3v x h/2 = 0 (6) y -u y h/2 - √ 3v y h/2 - √ 3x + √ 3u x h/2 + 3v x h/2 = 0 (7)
Solving with respect to x and y we get:

y = -hu y (8) x = h √ 3u x - √ 3v y -3u y + 3v x √ 3/6 (9)
The coefficients of the rotation matrix can be represented by quaternions. The quaternion representation is used for modeling the orientation as quaternions do not suffer from singularities as Euler angles do. The quaternion rotation matrix for the parallel robot is then

R =   2q 1 2 + 2q 2 2 -1 -2q 1 q 4 + 2q 2 q 3 2q 1 q 3 + 2q 2 q 4 2q 1 q 4 + 2q 2 q 3 2q 1 2 + 2q 3 2 -1 -2q 1 q 2 + 2q 3 q 4 -2q 1 q 3 + 2q 2 q 4 2q 1 q 2 + 2q 3 q 4 2q 1 2 + 2q 4 2 -1   (10)
with q 2 1 + q 2 2 + q 2 3 + q 2 4 = 1. In Equations 4, 6, 7, we substitute x, y using relations 8 and 9, and u, v, w by quaternion expressions using 10. Then ( 6) and ( 7) become q 1 q 4 = 0. Thus, we have either q 1 = 0 or q 4 = 0. This property is associated with the notion of operation mode [START_REF] Schadlbauer | Operation Modes in Lower-Mobility Parallel Manipulators[END_REF].

The notion of operation mode (OM) was introduced in [START_REF] Zlatanov | A Singularities as C-Space Singularities[END_REF] to explain the behavior of the DYMO robot. An operation mode is associated with a specific type of motion. For the DYMO, we have 5 operation modes: translational, rotational, planar (2 types) and mixed motions. In the workspace W , for each motion type, the W OMj is defined such that

• W OMj ⊂ W • ∀X ∈ W OMj , OM is constant
For a parallel robot with several operating modes, the pose can be defined by fixing the control parameters. For an operation mode OM j , if we have a single inverse kinematic solution, we can then define an application that maps X onto q: g j (X) = q (11)

Then, the images in W OMj of a posture q in the joint space Q is defined by:

g -1 j (q) = X | (X, q) ∈ OM j ( 12
)
where g -1

j
is the direct kinematic problem restricted to the operation mode j. Differentiating with respect to time the constraint equations leads to the velocity model:

A ṫ + B q = 0 ( 13
)
where A and B are the parallel and serial Jacobian matrices respectively, ṫ is the velocity of P and q is the joints velocity. The parallel singularities occur whenever det(A) = 0. Let OM 1 (reps. OM 2 ) be the operation mode where q 1 = 0 (reps. q 4 = 0). Then S OM1 and S OM2 are the loci of the parallel singularities and are characterized by: S OM1 : q 4 (8q 2 q 2 3 q 6 4 + 2q 2 q 8 4 -64zq 6 3 q 4 -96zq 4 3 q 3 4 -36zq 2 3 q 5 4 -6zq 7

The serial singularities occur whenever ρ 1 ρ 2 ρ 3 = 0. The common coordinates for both operation modes are z, q 2 and q 3 . Due to the redundancy of the quaternion representation, there exists two triplets defined by these three coordinates that represent the same pose in the same operation mode. To overcome this problem, we set q 1 > 0 and q 4 > 0. We can then depict a slice of this hypersurface by fixing one parameter as shown in Figure 2.

q 2 q 3 q 4 (a) q 1 q 2 q 3 (b)
Fig. 2 Singularity curves for z = 3, q 1 = 0 (a) and q 4 = 0 (b)

Aspect for an operation mode

In [START_REF] Wenger Ph | Definition Sets for the Direct Kinematics of Parallel Manipulators[END_REF], the notion of aspect is defined for parallel robots with only one inverse kinematic solution. An aspect W A i is a maximal singularity free set defined such that:

• W A i ⊂ W • W A i is connected • ∀X ∈ W A i , det(A) = 0 and det(B) = 0
This notion is now extended for a parallel robot with several operation modes such that:

• W A ij ⊂ W OMj • W A ij is connected • ∀X ∈ W A ij , det(A) = 0 and det(B) = 0
In other words, an aspect W A ij is the maximum connected region without any singularity of the OM j . The analysis of the workspace is done in the projection space (z, q 2 , q 3 ), and shows the existence of four aspects as shown in Fig. 3. However, no further analysis is done to prove this feature in the four dimension space. As there are several solutions for the DKP in the same aspect, non-singular assembly mode trajectories are possible. The cylindrical algebraic decomposition (CAD) implemented in the SIROPA library has been used to decompose an aspect into a set of cells where algebraic equations define its boundaries [START_REF] Chablat | Uniqueness domains and non singular assembly mode changing trajectories[END_REF]. The CAD provides a formal decomposition of the parameter space in cells where the polynomials det(A) and det(B) have a constant sign [START_REF] Collins | Quantifier Elimination for Real Closed Fields by Cylindrical Algebraic Decomposition[END_REF] and the number of solutions for the DKP is constant.

z q 2 q 3 (a) z q 2 q 3 (b) z q 2 q 3 (c) z q 2 q 3 (d)

Characteristic surfaces for an operation mode

The notion of characteristic surface was introduced in [START_REF] Wenger | A new general formalism for the kinematic analysis of all nonredundant manipulators[END_REF] to define the uniqueness domains for serial Cuspidal robots. This definition was extended to parallel robots with one inverse kinematic solution in [START_REF] Wenger Ph | Definition Sets for the Direct Kinematics of Parallel Manipulators[END_REF] and with several inverse kinematic solutions in [START_REF] Chablat | Séparation des solutions aux modèles géométriques direct et inverse pour les manipulateurs pleinement parallèles[END_REF]. In this paper, we introduce this notion for a parallel robot with several operating modes. Let W A ij be one aspect for the operation mode j. The characteristic surfaces, denoted by S C (W A ij ), are defined as the preimage in W A ij of the boundary

W A ij of W A ij . S C (W A ij ) = g -1 j g(W A ij ) ∩ W A ij (16)
These characteristic surfaces are the images in the workspace of the singularity surfaces. By using the singularity and characteristic surfaces, we can compute the basic regions as defined in [START_REF] Wenger Ph | Definition Sets for the Direct Kinematics of Parallel Manipulators[END_REF]. The joint space is divided by the singularity surfaces in regions where the number of solutions for the DKP is constant. We also name these regions the basic components as in [START_REF] Wenger Ph | Definition Sets for the Direct Kinematics of Parallel Manipulators[END_REF]. For each operation mode, we find regions where the DKP admits four (in red) or eight (in green) solutions, as it is depicted in Fig. 4. We can also notice in

Non-singular assembly mode changing trajectories

Due to the lack of space and for pedagogical purpose, we only report a slice of the workspace. Letting z = 3, the basic regions are computed by using the cylindrical algebraic decomposition for a given aspect. Figure 5 shows the three basic regions' image of basic components with 8 solutions for the DKP and a single basic region's image of a basic components with 4 solutions for the DKP connects these three previous basic regions. Table 1 presents the roots of the DKP for det(A) > 0 for a joint position in each operation mode.

For each of them, we find out that three roots have their z coordinate close to 3. A non-singular assembly mode changing trajectory can be obtained between three basic regions coming from eight solutions to the DKP. Due to symmetrical properties, there are also three roots of the DKP for det(A) < 0 with z = -3. For OM 1 , we construct a path between P 1 , P 2 , P 3 and for OM 2 between P 5 , P 6 and P 7 . When a straight line between two poses cross a singularity, we add an intermediate point as shown in Fig. 5. The connections between the basic regions depicted in red are the projections of the cusp points in the workspace, i.e. the tangent between the singularity surface and the characteristic surface [START_REF] Chablat | Uniqueness domains and non singular assembly mode changing trajectories[END_REF]. The variation of the det(A) is plotted in the Table 1 Solutions of the DKP for det(A) > 0 OM 1 OM 2 ρ 1 = 3.90, ρ 2 = 3.24, ρ 3 = 3.24 ρ 1 = 3.79, ρ 2 = 3.24, ρ 3 = 3.24 P z q 2 q 3 q 4 P z q 1 q 2 q 3 P 1 Figure 6 and shows the existence of a non-singular assembly mode changing trajectory. The image of this trajectory in the joint space is illustrated in the Figure 7. The projection of the cyclic trajectory defined by (P 1 , P 2 , P 3 , P 1 ) (resp. (P 5 , P 6 , P 7 , P 5 )) onto the joint space encloses three curves of cusps. This behavior is similar to that of the 3-RPR robot described in [START_REF] Zein | Non-Singular Assembly-mode Changing Motions for 3-RPR Parallel Manipulators[END_REF] or the 3-RPS robot in [START_REF] Husty | Non-Singular Assembly Mode Change of 3-RPS Manipulators[END_REF]. The path to connect the fourth solutions is not presented in this paper. The method introduced in [START_REF] Moroz | On the determination of cusp points of 3-RPR parallel manipulators[END_REF] is used to compute the cusp curves.

q 2 q 3 P 1 P 2 P 3 S OM 1 S C (a) q 2 q 3 P 5 P 6 P 7 S OM 2 S C (b)

Conclusions

This article presents a study of the joint space and workspace of the 3-RPS parallel robot and shows the existence of non-singular assembly mode changing trajectories. First, we have shown that each of the two operation modes is divided into two aspects, which is a necessary condition for non-singular assembly mode changing trajectories. Moreover, it turns out that this mechanism has a maximum of 16 real solutions to the direct kinematic problem, eight for each operation mode. Then, by computing the characteristic surfaces, we have shown that we can describe the basic regions for each operation mode. We construct a path going through several basic regions which are images of the same basic component with 8 solutions for the DKP. The analysis of the determinant of Jacobian shows that a non-singular assembly mode change exists for each motion type.

Fig. 3

 3 Fig. 3 Aspects for OM 1 with det(A) < 0 (a) and det(A) > 0 (b) and aspects for OM 2 with det(A) < 0 (c) and det(A) > 0 (d)

  Fig 4 the existence of cusp points.

Fig. 4

 4 Fig. 4 Slice of the joint space for ρ 1 = 3 for OM 1 (a) and OM 2 (b)

Fig. 5 Fig. 6 Fig. 7

 567 Fig. 5 Slice of an aspect for z = 3 and det(A) > 0 for OM 1 (a) and OM 2 (b) with in blue (resp. in red) a basic region coming from a basic component with four DKP (resp. eight)

-24z 2 q 2 q 2 3 q 2 4 -6z 2 q 2 q 4 4 -32q 2 q 2 3 q 4 4 -10q 2 q 6 4 + 2z 3 q 3 4 + 96zq 4 3 q 4 +72zq 2 3 q 3 4 + 23zq 5 4 + 16z 2 q 2 q 2 3 + 10z 2 q 2 q 2 4 + 8q 2 q 4 4 -z 3 q 4 -36zq 2 3 q 4 -21zq 3 4 -4z 2 q 2 + 4zq 4 ) = 0 (14) S OM2 : q 2 1 (6q 7 1 q 3 + 8q 5 1 q 3 3 -2zq 6 1 + 36zq 4 1 q 2 3 + 96zq 2 1 q 4 3 + 64zq 6 3 -18z 2 q 3 1 q 3 -24z 2 q 1 q 3 3 -18q 5 1 q 3 -16q 3 1 q 3 3 + 2z 3 q 2 1 + 3zq 1 4 -72zq 2 1 q 2 3 -96zq 4 3 + 18z 2 q 1 q 3 + 12q 3 1 q 3 -z 3 + 3zq 2 1 + 36zq 2 3 -4z) = 0 [START_REF] Zein | Non-Singular Assembly-mode Changing Motions for 3-RPR Parallel Manipulators[END_REF]