K V Kostas

A G Ginnis

C G Politis

P D Kaklis

A I Ginnis

Shape-optimization of 2D hydrofoils

Keywords: Isogeometric Analysis, NURBS, Potential Flows, Lifting Flows, Shape optimization

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Performance of many ship-hull types is dependent on the hydrodynamic properties of their keels, rudders, hydrofoils and similar appendages operating as lifting bodies. Therefore, the interaction of free-surface flows with lifting bodies constitutes an interesting problem, finding applications in the design of yachts and sailing boats and the performance of stabilizers, hydrofoils and similar devices. As a common feature, the functionality of these lifting appendages depends on their shape and even small variations of it may have significant impact on their performance. Fortunately, the current availability of computing power has allowed the employment of sophisticated Computational Fluid Dynamics solvers and their coupling with advanced Geometric Modeling techniques and Optimization strategies for the design optimization of both ship hull and lifting appendages. In this context, computational tools that represent exactly body shapes, securing that mesh generation is an error-free process provide designers and researchers with a great challenge. Moreover, the idea of bridging the gap between Computer-Aided Design (CAD) and Analysis by the introduction of IsoGeometric Analysis (IGA) by [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF], see also [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF] [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF], which directly uses analysis suitable geometric models from the CAD representation, can be efficiently exploited, especially in problems of shape optimization. Shape optimization in the context of IGA has been presented in various works, as e.g., [START_REF] Wall | Isogeometric structural shape optimization[END_REF] [START_REF] Wall | Isogeometric structural shape optimization[END_REF], [START_REF] Nagy | Isogeometric sizing and shape optimization of beam structures[END_REF] [START_REF] Nagy | Isogeometric sizing and shape optimization of beam structures[END_REF], [START_REF] Nguyen | Isogeometric shape optimization for electromagnetic scattering problems[END_REF] [START_REF] Nguyen | Isogeometric shape optimization for electromagnetic scattering problems[END_REF], Nortoft & Gravesen (2013) [START_REF] Nørtoft | Isogeometric shape optimization in fluid mechanics[END_REF], [START_REF] Qian | Isogeometric Shape Optimization of photonic crystals via Coons patches[END_REF] [START_REF] Qian | Isogeometric Shape Optimization of photonic crystals via Coons patches[END_REF], [START_REF] Cho | Isogeometric shape design optimization:exact geometry and enhanced sensitivity[END_REF] [START_REF] Cho | Isogeometric shape design optimization:exact geometry and enhanced sensitivity[END_REF], [START_REF] Lian | Sensitivity Analysis and Shape Optimisation through a T-spline Isogeometric Boundary Element Method[END_REF] [START_REF] Lian | Sensitivity Analysis and Shape Optimisation through a T-spline Isogeometric Boundary Element Method[END_REF], [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF] [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF], [START_REF] Li | Isogeometric Analysis and Shape Optimization via Boundary Integral[END_REF] [START_REF] Li | Isogeometric Analysis and Shape Optimization via Boundary Integral[END_REF], where NURBS control points and/or weights have been used as design variables to control the boundary shape. On the other hand, in [START_REF] Kostas | Ship-hull shape optimization with a T-spline based BEM-isogeometric solver[END_REF] [START_REF] Kostas | Ship-hull shape optimization with a T-spline based BEM-isogeometric solver[END_REF] and [START_REF] Ginnis | A Multi-Objective Optimization Environment for Ship-Hull Design based on a BEM-Isogeometric Solver[END_REF] [START_REF] Ginnis | A Multi-Objective Optimization Environment for Ship-Hull Design based on a BEM-Isogeometric Solver[END_REF], due to the shape complexity of the ship hull, a geometric parametric model using high-level parameters with a direct design meaning has been developed and used for optimization purposes in the context of IGA approach.

The main purpose of the present work is to develop an optimization procedure for the shape optimization of hydrofoils. This task is accomplished by combining three components: An optimization algorithm, an IGA-BEM solver for the potential flow outside the hydrofoil, presented in §2, and a hydrofoil geometric parametric model, presented in §3. [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF].

The formulation of the exterior potential-flow problem is reduced to a Boundary Integral Equation (BIE) for the associated velocity potential on the body boundary, exploiting the nullpressure jump Kutta condition at the trailing edge, [START_REF] Morino | Boundary Integral Equations in Aerodynamics[END_REF] [START_REF] Morino | Boundary Integral Equations in Aerodynamics[END_REF], [START_REF] Gennaretti | A high order boundary element formulation for potential incompressible aerodynamics[END_REF] [START_REF] Gennaretti | A high order boundary element formulation for potential incompressible aerodynamics[END_REF]. Boundary element methods (BEM) in aero/hydrodynamics were established as main tools for the solution of flow problems around lifting and/or non-lifting bodies of arbitrary shape; see, e.g., Katz & Plotkin (1991) [START_REF] Katz | Low Speed Aerodynamics[END_REF], [START_REF] Paris | Boundary Element Method: Fundamentals and Applications[END_REF] [START_REF] Paris | Boundary Element Method: Fundamentals and Applications[END_REF], [START_REF] Dragos | Mathematical Methods in Aerodynamics[END_REF] [START_REF] Dragos | Mathematical Methods in Aerodynamics[END_REF]. Starting from the pioneering work by [START_REF] Hess | Calculation of non-lifting potential flow about arbitrary 3-D bodies[END_REF] [START_REF] Hess | Calculation of non-lifting potential flow about arbitrary 3-D bodies[END_REF], where the 3D panel method based on source-sink distribution is presented for analyzing the flow around arbitrary non-lifting bodies, this approach has been further extended by [START_REF] Hess | Calculation of potential flow around arbitrary three dimensional lifting bodies[END_REF] [START_REF] Hess | Calculation of potential flow around arbitrary three dimensional lifting bodies[END_REF] for lifting flows and applied by [START_REF] Hess | Calculation of steady flow around propellers by means of a surface panel method[END_REF] [START_REF] Hess | Calculation of steady flow around propellers by means of a surface panel method[END_REF] for the simulation of a steadily translating and rotating propeller. Concentrating on marine propulsors, in the years to follow, various alternative formulations of the boundary element method have been applied for the solution of propeller related steady or unsteady problems. Most boundary element methods for lifting and propeller flows solve directly for the unknown velocity potential, a formulation more stable from the numerical point of view, see, e.g., [START_REF] Morino | Subsonic potential aerodynamics for complex configurations -A general theory[END_REF] [START_REF] Morino | Subsonic potential aerodynamics for complex configurations -A general theory[END_REF], [START_REF] Kinnas | Theory and numerical methods of hydrodynamic analysis of marine propellers[END_REF] [START_REF] Kinnas | Theory and numerical methods of hydrodynamic analysis of marine propellers[END_REF], although several alternatives exist, as for example the methods based on surface vorticity distributions developed by [START_REF] Belibassakis | A Boundary Integral Equation Formulation of the Neumann Problem for a Vector Field in R 3 with Application to Potential Lifting Flows[END_REF], 1998) [START_REF] Belibassakis | A Boundary Integral Equation Formulation of the Neumann Problem for a Vector Field in R 3 with Application to Potential Lifting Flows[END_REF], [START_REF] Belibassakis | A non-linear velocity based Boundary Element Method for the analysis of marine propellers in unsteady flow[END_REF]. In the last period, higher-order BEM, based on B-splines or NURBS representations for the geometry and/or the solution, have been presented for the analysis of flow around marine propellers; see, e.g., [START_REF] Lee | A B-spline higher-order panel method applied to two-dimensional lifting problem[END_REF] [START_REF] Lee | A B-spline higher-order panel method applied to two-dimensional lifting problem[END_REF], [START_REF] Kim | A B-spline based higher-order panel method for the analysis of steady flow around marine propellers[END_REF] [START_REF] Kim | A B-spline based higher-order panel method for the analysis of steady flow around marine propellers[END_REF], [START_REF] Gao | A NURBS based high-order panel method for three-dimensional radiation and diffraction problems with forward speed[END_REF] [START_REF] Gao | A NURBS based high-order panel method for three-dimensional radiation and diffraction problems with forward speed[END_REF]. In the present work, the IsoGeometric Analysis -Boundary Element Method (IGA-BEM) initiated by [START_REF] Politis | An Isogeometric BEM for exterior potential-flow problems in the plane[END_REF] [START_REF] Politis | An Isogeometric BEM for exterior potential-flow problems in the plane[END_REF], for the solution of an exterior potential-flow problem without lift, is extended and applied for the corresponding problem with lift. The IsoGeometric Analysis concept, introduced by Hughes et al (2005) [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] in the context of Finite Element Method, was extended to the Boundary Element Method by various authors, see, e.g., [START_REF] Politis | An Isogeometric BEM for exterior potential-flow problems in the plane[END_REF] [START_REF] Politis | An Isogeometric BEM for exterior potential-flow problems in the plane[END_REF], [START_REF] Simpson | A two-dimensional Isogeometric Boundary Element Method for elastostatic analysis[END_REF] [START_REF] Simpson | A two-dimensional Isogeometric Boundary Element Method for elastostatic analysis[END_REF], [START_REF] Scott | Isogeometric boundary element analysis using unstructured T-splines[END_REF] [START_REF] Scott | Isogeometric boundary element analysis using unstructured T-splines[END_REF], [START_REF] Belibassakis | A BEM-Isogeometric method for the ship wave-resistance problem[END_REF] [START_REF] Belibassakis | A BEM-Isogeometric method for the ship wave-resistance problem[END_REF], Peake et al (2013) [START_REF] Peake | Extended isogeometric boundary element method (XIBEM) for two-dimensional Helmholtz problems[END_REF], [START_REF] Simpson | Acoustic isogeometric boundary element analysis[END_REF] [START_REF] Simpson | Acoustic isogeometric boundary element analysis[END_REF], [START_REF] Ginnis | Isogeometric Boundary-Element Analysis for the Wave-Resistance Problem using T-splines[END_REF] [START_REF] Ginnis | Isogeometric Boundary-Element Analysis for the Wave-Resistance Problem using T-splines[END_REF]. The developed Isogeometric-BEM is based on a parametric NURBS representation of the hydrofoil and employs the very same basis of the geometry for representing the velocity potential. The Boundary Integral Equation is numerically solved by collocating at the Greville abscissas of the knot vector of the hydrofoil's parametric representation. Numerical error analysis of the Isogeometric-BEM, using h-refinement technique (knot insertion), is performed and compared with classical low-order panel methods.

As regards the geometric parametric model for the hydrofoil, it is constructed within Rhinoceros 3D using its VBscript-based programming language. This parametric model generates an instance of the hydrofoil as a closed cubic B-spline curve using a set of eight design parameters. The construction process is presented in detail in §3.2. The optimization loop is driven by an optimizer which uses a multi-objective optimization method based on a controlled elitist Genetic Algorithm (GA) favoring individuals that increase the diversity of the population; see §3.1. Finally, the optimization environment has been tested by optimizing the shape of a hydrofoil with respect to the following two criteria: maximum lift coefficient and minimum deviation from a reference cross-section area; see §3.3.

Continuous and discrete formulation of the problem

We consider a two-dimensional body whose boundary is ∂Ω B , as shown in Figure 1, moving with constant speed U B in an ideal fluid of infinite extent. In a body-fixed coordinate system Oxy this problem is equivalent to a uniform stream with velocity ∇Φ ∞ = U ∞ = -U B , where Φ ∞ (P) = u ∞ x + v ∞ y is the far-field asymptotic form of the velocity potential Φ(P) of the resulting flow at point P=(x, y). The potential Φ(P) is the solution of the following boundaryvalue problem (BVP):

∇ 2 Φ = 0, P = (x, y) ∈ Ω, (1)
∂Φ ∂n = 0, P ∈ ∂Ω B , (2)
Φ -

(u ∞ x + v ∞ y) → 0, as x 2 + y 2 → ∞, (3)
where Ω is the fluid domain outside ∂Ω B and n denotes the unit normal vector on ∂Ω B directed inwards with respect to the body. The above BVP has a unique solution up to an additive constant and, in order to fix a unique solution, we normally consider, for smooth bodies, zero circulation Γ(C) = C ∇Φ•dc of the velocity field ∇Φ over any circuit C surrounding the body. The difference between potential flows around a smooth body and a hydrofoil is that, in order for the flow around the hydrofoil to have a physical meaning, the circulation has to be nonzero and appropriately adjusted until the flow leaves the trailing edge smoothly. More specifically, on the basis of Kelvin's theorem, Prandtl concluded that if an airfoil, which started its motion from rest in an ideal fluid, is later found to possess non-zero circulation Γ, then the component of the boundary of the fluid which coincided with the airfoil initially, must coincide at a later time with the union of the airfoil surface and a surface, the so-called wake, embedded in the fluid which has circulation -Γ; see, e.g., [START_REF] Meyer | Introduction to Mathematical Fluid Dynamics[END_REF] [START_REF] Meyer | Introduction to Mathematical Fluid Dynamics[END_REF]. In contrast to the 3D case (see, e.g., [START_REF] Bassanini | A Boundary Integral Equation Formulation of the Neumann Problem for a Vector Field in R 3 with Application to Potential Lifting Flows[END_REF] [START_REF] Bassanini | A Boundary Integral Equation Formulation of the Neumann Problem for a Vector Field in R 3 with Application to Potential Lifting Flows[END_REF]), the location and shape of the wake in the 2D case can be taken, without loss of generality, to be a straight line emanating from the trailing edge and extending to infinity. This line is a force-free boundary along which the normal fluid velocity and the pressure should exhibit no jump. More accurately, we can write:

∂Φ + ∂n = ∂Φ - ∂n , P ∈ ∂Ω w : kinematic boundary condition, (4)
p + =p -, P ∈ ∂Ω w : dynamic boundary condition [START_REF] Nguyen | Isogeometric shape optimization for electromagnetic scattering problems[END_REF] where, Φ ± / p ± denote the velocity potential / pressure on the upper, ∂Ω w + , and lower, ∂Ω w -, face of the wake boundary, respectively.

Let Ω be the fluid computational domain bounded by ∂Ω = ∂Ω B ∪ ∂Ω w + ∪ ∂Ω ∞ ∪ ∂Ω w -with the radius of the fictitious circular boundary ∂Ω ∞ tending to infinity and the angle θ between the fictitious boundaries ∂Ω w + and ∂Ω w -tending to zero as r → ∞. Applying in Ω Green's second identity between the potential Φ(P), P∈Ω, and the fundamental solution, G(P, Q) = (1/2π) ln P -Q , of the 2D Laplace equation (1), we can reformulate the BVP (1)-(3) as a 2 nd -kind Fredholm integral equation on the hydrofoil boundary ∂Ω B , taking into account the wake-sheet conditions (4) and ( 5):

Φ(P) 2 + ∂Ω B Φ(Q) ∂G(P, Q) ∂n Q ds Q -µ w ∂Ωw ∂G(P, Q) ∂n Q ds Q = Φ ∞ (P), P ∈ ∂Ω B \P T E , (6)
where subscript Q implies an operator with respect to point Q, P T E = (x e , y e ) denotes the trailing edge of the hydrofoil and

µ w = Φ + (P) -Φ -(P) = constant, P ∈ ∂Ω w . (7)
Equation ( 7) is a form of the Kutta condition, the so-called Morino-Kutta condition (see [START_REF] Morino | Subsonic potential aerodynamics for complex configurations -A general theory[END_REF] [START_REF] Morino | Subsonic potential aerodynamics for complex configurations -A general theory[END_REF] , stating that the fluid velocity should be bounded in Ω and in particular near the trailing edge, enabling to fix the proper value of circulation Γ, thus securing a unique solution of the BVP (1)-( 3). Furthermore, due to the simple shape of the wake sheet, the integral term on the wake ∂Ω w , appearing in equation ( 6), can be easily calculated resulting in:

y x O U ∞ n ∂Φ + ∂n = ∂Φ - ∂n µ w = Φ + -Φ - ∂Ω ∞ ∂Ω B ∂Ω w + ∂Ω w - Ω ∇ 2 Φ = 0 ∂Φ ∂n = 0
∂Ωw ∂G(P, Q) ∂n Q ds Q = ∞ xe ∂G(P, (x Q , y Q )) ∂n Q dx Q = 1 2π arctan y P -y e x P -x e , (8)
revealing that the influence of the wake to the fluid flow is equivalent to a point vortex located at the trailing edge P T E = (x e , y e ). Alternatively, the flow around the hydrofoil can be formulated via the perturbation potential φ(P), defined as:

φ(P) = Φ(P) -Φ ∞ (P), (9)
which leads to the following integral equation:

φ(P) 2 + ∂Ω B φ(Q) ∂G(P, Q) ∂n Q ds Q -µ w ∂Ωw ∂G(P, Q) ∂n Q ds Q = - ∂Ω B ( U ∞ • n(Q))G(P, Q)ds Q , P ∈ ∂Ω B \P T E . ( 10
)
The range of the variable point P in the above integral equation, as well as in Eq. 6, does not include the trailing-point P T E in order to circumvent the issue that, due to the existence of the wake ∂Ω w emanating from P T E , one cannot simply appeal to standard potential theory for boundaries ∂Ω B with corners and replace the term Φ(P T E )/2 in the left-hand side of Equations 6 and 10 by [1 -( β/2π)]Φ(P T E ), β denoting the internal angle at the trailing edge. This choice controls the collocation method, adopted in §2.1, for solving the IGA-BEM discretization of [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF], by restricting the first and last collocation points from collapsing on P T E . There exists, however, works in pertinent literature, e.g., [START_REF] Morino | Singularities in BIEs for the Laplace equation: Joukowski trailingedge conjecture revisited[END_REF] [START_REF] Morino | Singularities in BIEs for the Laplace equation: Joukowski trailingedge conjecture revisited[END_REF], where P T E is permitted to be a collocation point in the context of a an enhanced formulation that involves an additional integral equation resulting from taking the derivative of the potential Φ(P), P ∈ Ω, with respect to the normal of the wake and letting P → P T E .

After solving [START_REF] Nørtoft | Isogeometric shape optimization in fluid mechanics[END_REF] or [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF], we can calculate the hydrofoil's pressure coefficient c p using Bernoulli's equation:

c p := p -p ∞ 1 2 ρ| U ∞ | 2 = 1 - v 2 t | U ∞ | 2 , ( 11
)
where v t is the tangential fluid velocity on the hydrofoil's ∂Ω B given by v t = ∂Φ/∂s with s denoting the arc length on ∂Ω B .

A discrete IGA-BEM formulation

The purpose of this subsection is to present a method that combines Boundary Element Method (BEM) with IsoGeometric Analysis (IGA) for solving numerically the boundary integral equations ( 6) and [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF]. IGA philosophy is equivalent to approximating the field quantities (dependent variables) of the boundary-value problem in question by the very same basis that is being used for representing (accurately) the geometry of the involved body-boundary. In our case the dependent variables are the total potential Φ(P), P ∈ ∂Ω B (see eq. ( 6)) or the perturbation potential φ(P), P ∈ ∂Ω B (see eq. ( 10)). For this purpose, we shall presume that the body boundary ∂Ω B can be (accurately) represented as a closed parametric NURBS curve r(t), t ∈ [0, 1], which is regular, i.e., the derivative vector is well defined and not vanishing, with the exception of the trailing edge: r(0) = r(1), where the derivative vector is not defined. More specifically,

r(t) = (x(t), y(t)) := n i=0 d i M i,k (t), t ∈ I = [t k-1 , t n+1 ] := [0, 1], ( 12
)
where {M i,k (t)} n i=0 is a rational B-spline basis of order k, defined over a knot sequence J = {t 0 , t 1 , . . . , t n+k } and possessing non-negative weights w i , i = 0, . . . , n, while d i are the associated control points; see, e.g., Piegl & Tiller [START_REF] Piegl | The Nurbs Book[END_REF]. Equation [START_REF] Qian | Full analytical sensitivities in NURBS based isogeometric shape optimization[END_REF] can then be written in the following form:

φ(t) 2 + I φ(τ )K(t, τ )dτ - µ w 2π arctan y(t) -y e x(t) -x e = g(t), t ∈ (0, 1), (13)
where, for the sake of notational simplicity, we define φ(t) 13) is a function expressed as:

:= φ(r(t)), G(t, τ )) := G(r(t), r(τ )) and K(t, τ ) = (∂G(t, τ )/∂n τ ) ṙ(τ ) , g(t) = -I ( U ∞ • n(τ ))G(t, τ ) ṙ(τ ) dτ . The kernel K(t, τ ) of the integral equation (
K(t, τ ) = ṙ(τ ) × (r(t) -r(τ )) 2π r(t) -r(τ ) 2 if t = τ, (14)
where here the symbol × denotes the so-called "2D exterior product" defined as: (x 1 , y 1 ) × (x 2 , y 2 ):=x 1 y 2x 2 y 1 . Assuming that r(τ ) ∈ C 2 (I), a careful limiting process as τ → t involving a 2 nd -order Taylor expansion of r(τ ) around τ = t, gives:

K(t, t) = ṙ(t) × r(t) 2π ṙ(t) 2 , if t = τ, (15)
revealing that in the neighborhood of regular points of a smooth (at least C 2 ) ∂Ω B , the kernel K(t, τ ) is a regular function.

Aiming to employ IGA for handling equation ( 13), we project, in a suitably defined manner, the perturbation potential φ(t) on the spline space S k (J ( ) ), S k (J (0) ) := S k (J ), expressed in the form:

φ s (t) := P s (φ(t)) = n+ i=0 φ i M ( ) i,k (t), t ∈ I, M (0) i,k (t) := M i,k (t), (16)
where ∈ N 0 denotes the number of knots inserted in I. Recalling the fundamental property of knot insertion, we can say that {S k (J ( ) ), ∈N 0 } constitutes a sequence of nested finite dimensionalspaces, i.e., S k (J ( ) ) ⊂ S k (J ( +1) ). Equation ( 13) can then be written as

1 2 n+ i=0 φ i M ( ) i,k (t) + I n+ i=0 φ i M ( ) i,k (τ )K(t, τ )dτ - µ w 2π arctan y(t) -y e x(t) -x e = g(t), t ∈ (0, 1). ( 17
)
Several methods are available for defining the projection P s (see eq.16) onto the finite-dimensional space S k (J ( ) ) and discretizing equation ( 13), like Galerkin and collocation. In the present work, a collocation scheme is adopted, which consists in projecting on S k (J ( ) ) through interpolation at a set of collocation points t = t j , j = 0, . . . , n + , which are chosen to be the Greville abscissas associated with the knot vector J ( ) . This leads to the following linear system for the unknown coefficients φ i , i = 0, . . . , n + :

1 2 n+ i=0 φ i M ( ) i,k (t j )+ n+ i=0 φ i q i (t j )- (φ n+ -φ 0 ) 2π arctan y(t j ) -y e x(t j ) -x e = g(t j ), j = 0, . . . n+ , (18)
where q i (t j ) = I M ( ) i,k (τ )K(t j , τ )dτ while, taking into account [START_REF] Qian | Isogeometric Shape Optimization of photonic crystals via Coons patches[END_REF], the unknown µ w is expressed through the difference φ n+φ 0 . Since t ∈ (0, 1), one must consider shifting the values of the first and the last Greville abscissa by a small value > 0. This shifting process gives rise to the need of proper handling of the numerical evaluation of the integral terms q 0 (1 -) and q n+ ( ) appearing in [START_REF] Dragos | Mathematical Methods in Aerodynamics[END_REF]. Indeed, by using the fact that r(0)=r(1) and incorporating 1 st -order Taylor expansion of r( ), ṙ( ) and r(1 -), ṙ(1 -) around 0 and 1, respectively, equation ( 14) leads to the following asymptotic estimates for the kernel:

K(1 -, ) = O( 1 ), K( , 1 -) = O( 1 ), K(1 -, 0) = O( 1 ), K( , 1) = O( 1 ). (19)
Thus, since M 1), we conclude that the integrands of q 0 (1 -) and q n+ ( ) are of order of 1/ in a neighborhood of 0 and 1, respectively. Actually, as one can see in the attached media file 1 , this behavior is restricted in an neighborhood of 0 and 1 and the graph of the integrands have a bell-like shape with a height of O(1/ ) and width of O( ). These integrals can be handled by an adaptive numerical integration scheme provided that precautions are taken in order that does not fall under a threshold value for which the integrals cannot be numerically computed.

( ) 0,k ( )=O(1), M ( ) n+ ,k (1 -)=O(
The solution of the linear system [START_REF] Dragos | Mathematical Methods in Aerodynamics[END_REF] provides the values φ i , i = 0, . . . , n + , which can be used to easily calculate the tangential velocity on the hydrofoil by using the derivatives of the rational B-spline basis functions M ( ) i,k (t), without resorting to the finite difference scheme, as it is the the practice in low-order panel methods:

v t = U ∞ • t + ∂φ ∂s = U ∞ • t + 1 ṙ(t) n+ i=0 φ i dM ( ) i,k (t) dt , (20)
where t denotes the unit-tangent vector of r(t).

1 See online version of the document.

Numerical results & optimization examples

As a first attempt to test the accuracy and effectiveness of the IGA-BEM method, presented in §3.3, we present numerical results for the flow around the standard NACA-4412 profile. This hydrofoil is represented as a cubic B-spline curve with 47 knots and 43 control points. In Figure 2, the total potential Φ on the hydrofoil is presented for an incident flow at an angle of attack of 5 deg. Numerical results of the present IGA method are compared with corresponding results of low-order panel method [START_REF] Moran | An introduction to theoretical and computational aerodynamics[END_REF] [START_REF] Moran | An introduction to theoretical and computational aerodynamics[END_REF] using a large number of elements (N=500). The agreement of IGA solution, using considerably less degrees of freedom (n=120), is very satisfactory. We recall that in IGA method degrees of freedom (DoF) is the number of the control coefficients used in the NURBS approximation of the potential Φ; see eq. ( 16).

In Figure 3, we compare the pressure coefficient distribution around a NACA-4412 profile at 7 deg. angle of attack, with corresponding results with a low-order panel method which adopts either the direct (Morino) formulation (see [START_REF] Moran | An introduction to theoretical and computational aerodynamics[END_REF]), as in the present work, or the classical indirect formulation of Hess & Smith [START_REF] Hess | Calculation of potential flow around arbitrary three dimensional lifting bodies[END_REF]. Again, considering the low-order panel solution with a large number of elements (N=500), we can see the good performance of IGA method, where it is needed considerably less DoF (n=120) for the solution to converge.

In Figure 4, we compare the overall performance of IGA-BEM compared to low-order BEM. This is accomplished by plotting, for both methods, the error of the developed circulation of the flow µ w around the hydrofoil (we recall that µ w = φ n+ -φ 0 ; see eq. ( 7)) as a function of DoF. Since the exact circulation around the hydrofoil is not known, the error calculated is based on a reference solution, µ wr , which is computed using the IGA-BEM method for a large number of DoF ( 5000). In IGA-BEM refinement is performed by knot insertion (h-refinement). It can be observed that IGA-BEM performs better compared to low-order panel method (e.g. for a relative accuracy of 0.001 the IGA approach requires around 120 DoFs, while for the low-order panel method this rises to 1000 DoFs; see also Figures 2 and3). For this level of accuracy the computational time is 17 and 24 seconds for the IGA and panel method respectively, when executed on a standard desktop PC. The rate of convergence exhibited by the IGA method is superior to that of the panel one but not analogous to the corresponding results concerning potential flows without circulation (see [START_REF] Politis | An Isogeometric BEM for exterior potential-flow problems in the plane[END_REF]) and application of the present IGA-BEM method in elasticity problems (see, e.g., [START_REF] Li | Isogeometric Analysis and Shape Optimization via Boundary Integral[END_REF] and [START_REF] Simpson | A two-dimensional Isogeometric Boundary Element Method for elastostatic analysis[END_REF]). This behavior may be attributed to the trailing edge, where the body boundary ceases to be smooth and the wake sheets emanate.

Optimization environment

Simulation-based optimization is of growing importance in engineering, since it allows to improve a component's or system's performance for moderate cost, in comparison with generally costly physical experiments. Moreover, the optimization is conducted in a rigorous algorithmic framework that enhances the experience and intuition of designers. A major difficulty to apply an automated shape optimization for complex engineering systems is the development of a fully automated design loop. Specifically, for each optimizer-generated set of parameters, a valid geometric model has to be constructed, allowing the generation of the computational domain used by the solver to provide the physical response and the performance analysis which needs to be returned back to the optimizer. All these steps should be fully automated, without user intervention, in order to feed the optimization algorithm and finalize the design loop. In this context, the IGA paradigm offers a significant improvement over the classical panel-based methods, since it relies on a direct relationship between the design parameters and the solver, without any geometrical intermediate meshes.

A second obstacle arises from the simulation process: for complex test-cases, CFD simulations are expensive, in terms of computational time. Moreover, the numerical solutions obtained can be polluted by errors arising from the discretization and iterative methods, yielding noisy performance evaluations. Sometimes, this may lead the optimizer to spurious local optima or even yield the failure of the optimization procedure. Here again, the isogeometric context may be helpful because it allows to avoid geometrical approximations, which reduces the error level, and permits to construct high-order solutions yielding a better computational efficiency.

In our case, the optimization loop comprises three components: the optimization algorithm, the IGA solver presented in §2 and the hydrofoil parametric model which will be presented in the following subsection. Hydrofoil shape optimization examples using the developed framework are presented in §3.3. Both the optimizer and IGA solver are implemented in matlab while the hydrofoil parametric model is implemented in Rhinoceros 3D modeling package.

The selected optimization algorithm belongs to the category of evolutionary ones, as experimentation with gradient and hessian-based algorithms has indicated the existence of multiple local minima that makes their usage problematic. Specifically, our optimizer uses the multi-objective optimization method gamultiobj which employs a controlled elitist genetic algorithm (a variant of NSGA-II [START_REF] Deb | Multi-Objective Optimization Using Evolutionary Algorithms[END_REF]). An elitist Genetic Algorithm (GA) always favors individuals with better fitness value (rank). A controlled elitist GA also favors individuals that can help increase the diversity of the population even if they have a lower fitness value. It is important to maintain the diversity of population for convergence to an optimal Pareto front.

Parametric Model for a Hydrofoil

The parametric model for a general hydrofoil has been materialized within Rhinoceros 3D modeling software package with the aid of its VBscript-based programming language, Rhinoscript. Our model generates a closed cubic B-Spline curve that represents a hydrofoil, using a set of 8 parameters; see Table 1. The geometrical definition of these parameters is depicted in Figure 5. All parameters, with the exception of chord's length (L), are defined using appropriate nondimensional ratios so that their values always lie in [0, 1]. This approach eliminates the need of implementing complex interdependent constraints while guaranteeing the robustness of the procedure which is of significant importance in an optimization procedure. The dimensionalized range of the parameters is appearing in the "actual range" column of Table 1.

The parametric model's construction initiates with the definition of four simple Bézier curves that are employed in the generation of the final cubic B-Spline curve. The upper part of the hydrofoil comprises two parts: a simple cubic Bézier curve and a simple quadratic one for the fore and aft part respectively. The first control point of the cubic curve is positioned at (0, 0) while the last one is at (x z max, z max). The first interior point lie on the line segment defined by (0, 0) and (0, z max) while the second one on the segment defined by (0, z max) and (x z max, z max). The exact position of these points is controlled by the tip form factor parameter. The second part is a quadratic simple Bézier adjacent to the first part and ending on the trailing edge, i.e., at (L, 0). Its interior control point is positioned at ( L-z max tan(a b) , z max) which ensures first order geometric continuity at (x z max, z max). The camber curve is constructed in a similar way using two quadratic simple Bézier curves. The aft part construction is equivalent to the aft part of the hydrofoil's suction side with x z max, z max and a b being replaced by x c max, c max and a b p. For the fore part, an angular parameter is introduced, which corresponds to the tangent of the camber curve at the leading edge and uniquely defines the position of the interior control point. However, in our optimization examples this parameter has been kept fixed and is not included in the list of Table 1.

Finally, the pressure side is constructed by mirroring the upper part about the camber curve. The resulting three2 curve pieces can be exactly represented with a single cubic B-Spline curve, if quadratics are degree elevated and multiple knot insertion is performed at their endpoints. Alternatively, if knot multiplicity is to be avoided, the hydrofoil can be approximated by a single B-Spline curve using Rhinoscript's FitCurve3 function.

The developed parametric model is sufficiently flexible for our shape optimization purposes as it can easily generate specific hydrofoils, as the NACA-4412 depicted in Figure 6, or largely varying general hydrofoils; see shape instances in Figure 7. Figure 7(A) depicts a hydrofoil where the longitudinal positions of camber's and suction side's max-width have moved towards the trailing edge. Figures 7(B,D,F,H) depict mainly variations in the values of z max and c max, while Figures 7(A,E,G) depict mainly variations in the relative positions of x z max and x c max. Finally, Figure 7(C) depicts a symmetric hydrofoil with c max equaling 0.

(A) (C) (E) (G) (B) (D) (F) (H)

Shape optimization examples

The optimization environment has been tested for optimizing the shape of a hydrofoil with respect to the following two criteria: a) maximum lift coefficient and b) minimum deviation from a reference area. The IGA-BEM solver, described in §2.1, is used to produce an average lift coefficient calculated for three angles of attack, namely 1, 3 and 5 degrees. The reference area in the second optimization criterion is set to be the one of the NACA-4412 profile. In the course of the optimization process, the parameter Length (L) of the hydrofoil parametric model is assumed to be fixed and is regularized to the value of one, while all remaining parameters, defined in Table 1, are considered to be free and contribute to the shape evolution of the hydrofoil.

As a first test, we consider that the values of the free parameters range between ±0.05 around those corresponding to the approximation of NACA-4412 by the parametric model (see Figure 6). Specifically, the min-max allowed values of the free parameters are given in Table 2. The multiobjective optimization method gamultiobj, offered by Matlab, is employed to calculate the Pareto front of the above optimization problem, which is depicted in Figure 8. The Pareto front is defined to be the set of design points that are Pareto efficient, i.e., it is impossible to make any one of the objective function values better off without making the other (or more generally, at least one of the others) worse off. Obviously, the optimum shapes lie on the Pareto front; see green line depicted in Figure 8. If, for example, we choose to accept a 0.52% deviation of the target area, the optimum solution would be the hydrofoil instance whose parameter values are depicted in Table 2. Table 2 includes both the non-dimensional and the corresponding dimensional (in parentheses) values of the parameters. The resulting "optimum" lift coefficient value is 1.768, i.e., ∼3% better than that of NACA-4412. Finally, Figure 9, depicts the "optimum" instance of the hydrofoil along with an approximation of NACA-4412. As a second test, we consider that the range of variation of the free parameters is not restricted and is the whole of the interval [0, 1]. The resulting Pareto front is depicted in Figure 10. Finally, Figure 11 depicts six instances of the hydrofoil parametric model lying on the Pareto front. The hydrofoil instances in Figure 11 correspond to the points marked with the circled numbers in Figure 10. This second example leads to shapes that are prone to flow separation even for small Reynold's numbers. Hence, the potential flow theory should be carefully applied in this case and we may need to consider a different approach, e.g., boundary layer corrections or viscous flows, if we are to perform shape optimizations in large feasible spaces.

Conclusions and future work

In the present work, an optimization procedure has been developed and used for the shape optimization of hydrofoils. The optimization environment is based on a geometric parametric modeler for the hydrofoil, an Isogeometric BEM solver for the potential flow and an optimizer employing a controlled elitist genetic algorithm. This procedure has been successfully applied for the multi-objective shape optimization of a hydrofoil with respect to the criteria of maximum lift coefficient and minimum deviation of the hydrofoil area from a reference area. The isogeometric BEM solver developed for the steady lifting flow around hydrofoils is robust and effective since it needs much less degrees of freedom for the same accuracy, compared to classical panel methods. Further work is needed for the improvement of the rate of convergence of the solution by adopting adaptive refinement strategies as well as by treating properly the Boundary Integral Equation near the trailing edge in order to take into account both the corner point of the body boundary and the existence of the wake. As future work, the extension of the present methodology to 3D problems seems to be straightforward. Moreover, the adoption of boundary-layer corrections to the potential flow should enrich the obtained solutions, from the physical point of view. More challenging, seems to be the problem of shape optimization of 3D hydrofoils moving near the free surface since such systems face major difficulties due to the simultaneous presence of the free surface and the free-vortex systems generated by the lifting bodies which can lead to non-linear phenomena, as, e.g., wave breaking.

Figure 1 :

 1 Figure 1: A hydrofoil in a incoming uniform stream with velocity U∞ along with its wake boundary ∂Ωw. Dashed lines denote the assumed boundary orientation.

Figure 2 :

 2 Figure 2: Potential distribution around a NACA-4412 profile at 5 deg. angle of attack and comparison with low-order panel method (L denotes the chord length).

Figure 3 :Figure 4 :

 34 Figure 3: Pressure distribution around a NACA-4412 profile at 7 deg. angle of attack and comparison with low-order panel methods

Figure 5 :

 5 Figure 5: Parametric model of a hydrofoil and its defining parameters

Figure 6 :

 6 Figure 6: Approximation of a NACA-4412 profile using the parametric model. Max. deviation 0.1% of chord length

Figure 7 :

 7 Figure 7: Various instances of the hydrofoil parametric model

Figure 8 :

 8 Figure 8: Pareto front for the inverse of the average lift coefficient, C -1 , and the area-deviation criteria. Average lift coefficient computed for 1, 3 and 5 degrees of angle of attack. Parameters ranging ±0.05 around those corresponding to NACA-4412.

Figure 9 :

 9 Figure 9: NACA-4412 along with an optimized instance on the Pareto front of Figure 8.

1 Figure 10 :Figure 11 :

 11011 Figure 10: Pareto front for the inverse of the average lift coefficient and the area-deviation criteria. Average lift coefficient computed for 1, 3 and 5 degrees of angle of attack.

Table 1 :

 1 Parameters' definition

	Nr. Name	description	symbol	actual range
	1	Length	Length of hydrofoil's chord	L	free
	2	Max width	Maximum width of suction	max z
			side wrt chord

Table 2 :

 2 NACA-4412 optimization parameters

	Nr. Parameter value for NACA4412	range of variation	optimized
	1	L	1 (1)	fixed	1
	2	max z	0.487 (0.0984)	[0.437, 0.537] ([0.0882, 0.1083])	0.459
	3	max c	0.450 (0.0403)	[0.400, 0.500] ([0.0358, 0.0448])	0.423
	4	x z max	0.300 (0.3500)	[0.250, 0.350] ([0.3250, 0.3750])	0.259
	5	x c max	0.610 (0.4280)	[0.560, 0.660] ([0.4130, 0.4430])	0.588
	6	a b	0.070 (0.2647)	[0.020, 0.120] ([0.1954, 0.3340])	0.110
	7	a b p	0.445 (0.1275)	[0.395, 0.495] ([0.1151, 0.1399])	0.492
	8	tip	0.530 (0.5240)	[0.480, 0.580] ([0.4840, 0.5640])	0.511

Two for the suction side and one for the pressure side.

FitCurve creates a B-Spline curve of specified degree that matches the input curve(s) within a specified tolerance value.


Olivier Maurice
email: olivier.maurice@esigelec.fr

Adaptation of Kron's Tensorial Analysis of Network for the EMC Design and Analysis of Systems *

Keywords:

Gabriel Kron developed tensorial analysis of networks (TAN) in 1939 to deal with electrical machines. This formalism has been extended to applications in the Electromagnetic Compatibility (EMC) domain. Using such a topological approach to deal with complex systems is particularly relevant in the design phases compared to numerical approaches and allows a problem to be studied in an abstract way. In this paper the robust mathematical approach is synthesized. Then, after explaining how cords are added to the TAN method, we apply the formalism to three typical examples.

I. INTRODUCTION

Nowadays, there is an increasing number of electronic systems that are more and more complex. As systems integrate different functions, it becomes essential to prevent electromagnetic disturbances between the different parts. As a result, the design phase is crucial and must be carried out with precision to ensure of the correct functioning when manufacturing the system. One way to attain this objective is to use numerical simulation. For this, different numerical methods have given rise to numerous electromagnetic commercial codes. We can cite Finite Element, Finite Difference or Moment Methods that are very popular in the EMC domain. But when dealing with the design phase, the whole system has to be simulated to take all the interactions into account and to prevent strong parasitic coupling. As a result, it becomes difficult to detect the origin of disturbances responsible for malfunctions by only trying different configurations. This approach is time consuming, and supposes that the numerical tools are adapted to the problem. One alternative is to use the engineer knowledge to have a more complete vision of the system. In fact, as mentioned above, a system is generally consists of subsystems. Thus it is more efficient to have a macro-model vision of each subsystem and to analyse the different connections between them. Such a topological approach provides a metric giving the value for each connection and allows the comparison between the different propagation path of disturbances. To follow such an approach, Gabriel Kron, developed an engineering approach based on Tensorial Analysis of Networks. This approach first developed for electrical machines, has recently been extended to EMC problems [START_REF] Maurice | La compatibilité électromagnétique des systèmes complexes[END_REF] to take into account not only conducting disturbances but also radiation couplings and even non linear ones. Furthermore, this approach allows EMC problems to be studied in an abstract way. In section II, the Einstein notation used in this paper will be reviewed. In section III the principle of the classic Kron's method will be explained in details. Section IV presents three typical examples that show different kinds of models that can be integrated into Kron's approach. Then we will conclude on the interest of the method in section V.

II. NOTATIONS

In all of this paper we use Einstein's mute index notation. It means that each time an index is repeated, it implies a summation operation on the values of this index, i.e.:

III. KRON'S METHOD

Tensorial Analysis of Networks (TAN) as created by Kron [START_REF] Kron | Tensor analysis of networks[END_REF] starts from a graph that represents the engineer's thinking. From this graph, a system of non linear integrodifferential equations is deduced. This system is then studied theoretically. According to Kron, once equations are obtained, the circuit and its associated graph are no longer useful. The whole analysis should be conducted theoretically. Currents in the nodes pair space or the mesh space constitute the orthogonal axes of a multidimensional space that can be studied using differential geometry tools [START_REF] Hoffmann | Kron's non-Riemannian electrodynamics[END_REF]. The concept of topology in Kron's approach is a mathematical one, and is not directly linked to the physical appearance of the system. It calls for a higher degree of theoretical thinking. However it is not constrained by any good shielding approximation [START_REF] Parmantier | Topology based modeling of very large systems[END_REF] or any other specific approximations, for example [START_REF] Besnier | Electromagnetic topology: An additional interaction sequence diagram for transmission line network analysis[END_REF]. The topological side of Kron's method is in related to mathematical algebraic topology and its operators [7][8].

A. Graphs

In graphs nodes, branches and meshes are easily identifiable. These elements are part of a cellular topology [START_REF] Durand | Networks, Topology and Interaction Principle Implemented in the Kron's Method[END_REF]. Each branch is linked to a real electronic object and not at all to its geometrical aspect. Electronic objects are characterized by equations and intervals. A branch is associated with an equation that gives the relationship between an electromotive force e u applied to the object and the current i u flowing through the object under applied electromotive force e u = z uu (i u ). When engineers connect branches together to make a graph, they construct a real system represented by a drawing. Intervals specify minimum and maximum values of the parameters. A single graph is made of various connected networks (space R, where any node is connected to two branches). Each network has branches (space B) connected together by nodes (space N ), where a connected group of branches makes a mesh (space M ) [START_REF] Branin | The algebraic-topological basis for network analogies and the vector calculus[END_REF]. Euler-Poincaré's topological equation implies (see figure 1): Let us consider the case of a circuit made only with passive elements i.e. resistors, inductances and capacitances. In order to characterize an equivalent circuit associated with a branch, one connects an electromotive force across a branch and measures the current. Kron called this a primitive circuit. Using the Laplace operator, branch characterization reduces to simple expressions: u a = z aa i a (figure 2). On a branch, a difference in potential exists between the two ends of the branch but, making a mesh forces the difference in potential to be null. Connecting two branches to the same mesh (figure 2b.), the two differences in potential are equal (figure 2a.) and therefore, their difference all around the mesh is zero.

M -B + N -R = 0 (2)

FIG. 2. Simple primitive circuit

Lagrange's equation

Under our hypothesis and using the Laplace operator p, circuit laws become: E 0 = Ri, E 0 = Lpi, E 0 = 1/(Cp)i for a resistance, an inductance or a capacitance respectively. We now introduce how to couple these primitive circuits. If M, B, N, R are the topological characteristics of one circuit, we represent by i k , k = 1, 2, . . . , M the mesh currents associated with M meshes chosen in the circuit. e k , k = 1, 2, . . . , M are the electromotive forces of the meshes. We associate mesh loads q k with variables i k such as i k = qk . figure 3 shows two coupled meshes where inductances, elastances and resistances are associated with branches or meshes, mesh currents and loads. A cord (defined in more detail below) represents a mutual inductance between meshes.

FIG. 3. Example of network for Lagrange's equation

In the whole circuit, cinetic energy is given by:

T = 1 2 L km qk qm (3)
For the potential energy we have:

U = 1 2 S km q k q m (4)
where S km are elastance values (inverse of capacitance). And finally for dissipations:

W = 1 2 R km qk qm (5)
Lagrange's equations for the circuit are given by:

d dt ∂T ∂ qk + ∂U ∂q k + ∂W ∂ qk = e k (6)
with (3),( 4),(5) and using Laplace's operator p and transformations on i and e, with i k 0 initial conditions in current and q k 0 initial conditions in loads, Lagrange's equation becomes:

L km p (i m -i m 0 ) + R km i m + S km i m p + q m 0 = e k (7)
Defining: Z km = L km p+ S km p +R km , this is equivalent to:

Z km i m -L km i m 0 + S km q m 0 = e k .
Usually, L kk are the self inductances of meshes, L km , k = m the mutual inductances between meshes. S kk the self elastances of meshes and S km , k = m common elastances between meshes. R kk are the self dissipations of meshes and R km , k = m common dissipations between meshes. Note that for null initial conditions, Lagrange's equation can be written:

L km p + S km p + R km i m = e k (8)

Invariant

The fundamental concept of tensorial analysis is the existence of an invariant. In the case of Kron's proposal, this invariant is:

ω = e k i k = f k h k (9)
for two different referentials e k , i k and f k , h k . Suppose a change of referential to use another mesh base defined by:

i a = C a k h k (10)
h k being new mesh currents. Using [START_REF] Denis-Papin | Cours de calcul tensoriel appliqué[END_REF] in [START_REF] Branin | The algebraic-topological basis for network analogies and the vector calculus[END_REF] gives f k = e a C a k . Now using [START_REF] Denis-Papin | Cours de calcul tensoriel appliqué[END_REF] in (8) leads to:

Z km C m a h a = e k (11)
Multiplying the two terms by C i n gives (with the relation on f k ):

C k n Z km C m a h a = C k n e k = f n (12)
As C k n are constant, in general C k n Z km C m a = Z na , which transforms [START_REF] Maurice | On mathematical definition of cords between networks[END_REF] into:

Z na h a = f n ( 13
)
The equation is unchanged under any choice of mesh base. It means that the mesh space, under our hypothesis, can be considered as a Riemannian space with a Euclidian metric as all the Z km are constant and Z km is its fundamental tensor [10][11].

3. Application of Kron's method to EMC: Generalized interactions through cords Each time a current i o creates a field, this field induces an electromotive force e k at the considered part of the network. This field is coupled after a transport process associated to any adequately defined Green's function. As for branches, these functions can be anything as long as their mathematical expression is defined. Starting from this proposal, the generic name of cord [START_REF] Maurice | On mathematical definition of cords between networks[END_REF][13] was submitted to represent the coupling functions between branches or meshes. Green's functions can be those of free space, inside cavities between wires and field modes, and many other interactions. The coupling can be associated with the source covector e k or represented using an impedance coupling function obtained from the ratio e k /i m = z km . From this general approach of coupling functions, several models were developed, such as [START_REF] Leman | Simulation d'une cavité EM par réduction de circuits électriques couplés formant une super matrice impédance[END_REF]: i Branin's model for transmission lines ii field to wire interactions iii coupling to objects within cavities These physical processes are the major ones for EMC work. But also unusual EMC problems can be addressed, such as the interaction between particles and power supply on electronics [START_REF] Hubert | Combined effects induced by electromagnetic pulses and single event effect on electronics[END_REF] or the emissions of printed circuit board tracks [START_REF] Rogard | Comparison of radiation modeling techniques up to 10 GHz Application on a microstrip PCB trace[END_REF]. Any near field or far field interaction can be easily incorporated in Kron's formalism using cords. Many cases were tested with good results providing a very elegant formulation of the problem under the mesh space. One particular example to illustrate this is the near field interaction between antennas achieved in [START_REF] Alaeldine | EMC-oriented analysis of electric near-field in high frequency[END_REF].

A good illustration for general use of cords can be the coupling between two matched antennas in free space. The first mesh consists of the source generator, source impedance and radiation resistance. The second mesh is similar to the first but without a generator and a matched load. The cord is defined by Friis's equation:

P 2 P 1 = G 1 G 2 λ 4πR 12 2 e -R 12 c p (14)
where P 2 is the received electromagnetic power on the second antenna, P 1 the transmitted power to the first antenna, G 1 the first antenna gain, G 2 the second one, R 12 the distance between the two antennas, λ the wavelength, c the speed of light and p the Laplace operator. To define the cord z 21 we define e 2 as a received electromotive force and i 1 as a transmitted current. Since P 1 = Z 0 (i 1 ) 2 and P 2 = (e 2 ) 2 /Z 0 , Z 0 being the radiation resistance and source impedance, we easily obtain: For antennas, the coupling is symmetric. The fundamental tensor linked with this system is:

e 2 i 1 = z 21 = Z 0 λ 4πR 12 G 1 G 2 e -R 12 c p (15
Z = Z 0 z 21 z 21 Z 0 (16)
C. Intrinsic mesh properties

The advantage of having an additional space dimension compared to nodal techniques is to have a direct correspondence with surface flux through Hodge's operator [START_REF] Tarhasaari | Some realizations of a discrete Hodge operator: a reinterpretation of finite element techniques [for EM field analysis[END_REF]. There is no need for partial inductances: one inductance of a loop can be directly associated with this loop in the mesh space. A proposal to show the inductance in the center of the mesh in the graph was made [START_REF] Maurice | La compatibilité électromagnétique des systèmes complexes[END_REF]. Mathematically, once the impedance matrix Z is defined in the branch space, its transformation by the matrix C to the mesh space is carried out using the bilinear product ζ = [C] T ZC. Now if some properties belong to the mesh space, a µ matrix can be defined where self inductance of meshes, mutual inductances, etc., are computed and in the mesh space, the equation becomes:

e a = (ζ ab + µ ab ) i b (17)
D. Building a mesh set, additional branches and virtual meshes

On any graph G including nodes, branches, meshes and networks, one spanning tree ST can be defined for each network. It is a sequence of branches that go one and one time only through all the nodes of the network. If we impose to choose of a spanning tree made of successive branches, then if b u is an unused branch:

∃ (n 1 , n 2 ) ∈ {∂ST } / (n 1 , n 2 ) ∈ b u (∂ being the bound- ary operator). Then if A is the incidence matrix so that n 1 = A k 1 b k , n 2 = A q 2 b q and (b k , b q ) ∈ ST 2 . If {b k , .
. . , b q } is the set of branches included between b k and b q on ST , then {b k , . . . , b q , b u } makes a mesh M u using what is called a closing branch b u . This technique provides a guide to choose meshes in any network. It leads to more adequate matrix conditioning except when two meshes are coupled through a filter structure. figure 5 shows the two cases of a shared branch or a spanning tree process to define mesh-to-branch connectivity.

FIG. 5. Mesh definitions

In the case of the filter structure, it is clear that the boundary made by the filter must be the region where currents are shared and distributed depending on frequency and filter characteristics. We consider two branches of branch currents i 1 and i 2 , connected through a common mesh (they can be the boundaries of two larger networks). We can always add a branch of current i a connected to the same two nodes of branches i 1 and i 2 to make two meshes: the first one of currents i 1 , i a and the second one of currents -i a , i 2 . The impedance of this added branch appears as an extra-diagonal component of the impedance matrix of the two mesh networks. In many cases, this added branch has a capacitor as its impedance, which is the classical first order filter in EMC (figure 6).

FIG. 6. Shared branch detail

Another kind of mesh can be attached to the branch of ST : virtual meshes. These are meshes of known mesh currents J v [START_REF] Angot | Compléments de mathématiques à l'usage des ingénieurs de l'électrotechnique et des télécommunications[END_REF]. The general equation becomes: e a = z ab i b + z av J v .

E. Kron's approach

In Kron's approach, each component of a system is seen as a primitive element involved in the system construction [START_REF] Kron | Diakoptics: the piecewise solution of large-scale systems[END_REF]. The set of all primitive networks constitutes the diagonal part of Kron's fundamental tensor. All the electromagnetic interactions between the primitive elements remain to be added. Kron's fundamental tensor is a mathematical image of the EMC interaction matrix [START_REF] Morgan | A handbook for EMC testing and measurement[END_REF]. Kron's method can be broken down into six steps:

1. build up the set of primitive circuits z mn involved in the system as defined in the branch space (figure 7);

2. carry out the direct summation of all these components of impedance tensor z mm : z complete = ⊕ m z mm (figure 7);

3. add coupling between branches A mn (figure 8);

4. define connectivity C, choose a spanning tree to identify meshes then transform to the mesh space 9). Capital letters are now used to identify objects that belong to the mesh space; 5. add self mesh properties S µν , interactions between meshes and current sources: (Z µν + S µν ) = G µν , to obtain

Z µν = C m µ C n ν z mm | m=n + A mn | m =n (figure
E σ = G σβ i β + G µα J α (18)
(figure 10). ;

6. solve the system. The result of using this method is this system of equations. Kron's method is a rigorous mathematical basis to theoretically study the behavior of electronic systems in EMC. Once the system of equations is available, it can be solved using numerical techniques. In Kron's fundamental tensor G σβ , components coming from various techniques may be included: Z matrix for moment method, macro-models derived from measurements, and so on. Due to its mathematical origin, Kron's method may be applied to many different problems [START_REF] Leman | Use of the circuit approach to solve large EMC problems[END_REF] [START_REF] Kondo | Memoirs of the Unifying Study of the Basic Problems in Engineering Sciences by Means of Geometry: Part One[END_REF].

IV. THREE EXAMPLES OF KRON'S METHOD USE FOR EMC

In order to illustrate the use of Kron's method for EMC, we have chosen three examples. The first one focuses on a shielded cable problem. Internal pigtails still pose a difficult problem that cannot be solved using 3D codes or SPICE-like software. Kron's method is an efficient way to approach the problem, to write its mathematical expression rigorously and to compute the order of magnitudes. The second example focuses on cavities. Numerical software doesn't always take all losses into account inside cavities. Even if techniques are now available, it remains difficult to analyse the result without the order of magnitudes. Once again Kron's method gives the means to make this analysis. And finally the third case concerns abstract study. It shows how Kron's formalism leads to the rigorous analytical equations of

A. A multi-braid cable with pigtail

Working in the mesh space is an efficient technique to obtain equations in the case of cable transfer impedances. The technique was applied in a difficult case, where the internal braid was connected to the external one through a pigtail [START_REF] Casagrande | Shielded cables modelling with internal pigtails[END_REF]. As long as the wavelength is greater than the cable length, an external domain including the external braid, ground plane and shield short circuits can be reduced to one mesh. The same representation can be used for any other intermediate domain. figure 11 shows the cable configuration and figure 12 how the case was managed with two braids linked through pigtails.

The problem concerns three domains: the external domain with the external braid and the ground plane, the intermediate domain with the internal braid and pigtails and the internal domain with the internal braid and the wire. Each domain includes the braid resistance and self inductance and for the internal wire, the same elements but given for a wire inside a braid. Couplings are transfer impedances, computed using equations given in [START_REF] Demoulin | Blindages électromagnétiques[END_REF] and refer to equations [START_REF] Angot | Compléments de mathématiques à l'usage des ingénieurs de l'électrotechnique et des télécommunications[END_REF] to [START_REF] Leman | Use of the circuit approach to solve large EMC problems[END_REF]. For the direct coupling between the external braid and the wire through the internal pigtails, a special equation was developed in [START_REF] Casagrande | Shielded cables modelling with internal pigtails[END_REF] and it shows how the method easily accepts any kind of newly developed equation. The transfer impedance is given by:

Z pigtail = p µL p 2π Log a q (19)
L p being the pigtail length, a the supposed constant distance between the straight pigtail and the height reference q, located at the surface of the pigtail (see figure 13). The internal braid has a transfer impedance given by:

Z t2 = (R α + ξp) x ( 20
)
where x is the cable length, R α the DC resistance of the braid and ξ its transfer inductance value. However, since the internal braid is terminated by the internal pigtail, the total transfer impedance between this braid and the inner wire is given by Z t2 = Z t2 + Z pigtail ; so:

Z t2 = R α (x -L p ) + p ξ (x -L p ) + µL p 2π Log a q (21)
The external braid is coupled with the internal braid over its whole length and is coupled with the inner wire only along the pigtail length. Transfer impedance of braids is given, for braid angles around 45 degrees(ψ) and a diffusion process, by:

Z t = 4 πd 2 N Cσcos(ψ) (1 + j) d δ sh (1 + j) d δ ( 22
)
where σ is the braid conductivity, ψ the braid angle, C the number of conductor groups in the braid, d the conductor diameter and δ the skin depth given by:

δ = 2 ωµ 0 µ r σ 1 2 (23)
To the diffusion process must be added the diffraction process, represented by a transfer inductance. It is given, for angles around 45 degrees by:

L t = ν α m µ 0 π 2 D 2 ( 24
)
where D is the braid diameter, ν the number of holes per meter and α m the magnetic polarisability of each hole. α m is approximated by:

α m ≈ L 3 g 6 ( 25
)
L g being the hole length. A correction term ct is given by:

   Kp = 1,16 N Cd atan N 3 sin π 2 -2ψ µ0 σ ct = Kpe -j π 4 ω 0,7 (26)
and Z t1 is given by:

Z t1 = (Z t + ct + L t p) (x -L p ) ( 27
)
while Z t3 is given by Z t3 = (Z t + ct + L t p) L p . These equations provides the means to compute both Z t1 , Z t2 and Z t3 values. Z t3 is applied on the pigtail length L p .

The global impedance tensor of the problem is:

Z =   R 1 + R 2 + L 1 p 0 0 -Z t1 (x -L p ) L 2 p + R 3 + L 3 p 0 -Z t3 L p -Z t2 R 4 + R 5 + L 4 p   (28
) R k with k = 1, 2, 3, 4, 5 are the loads. For the internal braid, the only significant resistance is the pigtail connection to the end of the cable. L 3 is the pigtail inductance (approximated by µπr/2, with r radius of the equivalent circle made by the pigtail and its distance to the inner wire).

The measurement was made in a triaxial test bench for transfer impedances. Inductances L k for k = 1, 2, 4 were computed using the fundamental coaxial laws. The first one gives the cable characteristic impedance:

Z c = 60ln r e r i ( 29
)
r e and r i are the internal and external radii of the braid and the internal wire. Then, we can obtain the inductance using:

L = Z c v x ( 30
)
where v is the propagation speed in the cable. The measurement was made on a 10 cm diameter and 30 cm long triaxial test bench, which can make measurements up to 100 MHz [START_REF] Casagrande | Shielded cables modelling with internal pigtails[END_REF]. figure 14 shows the measurements and computation curves. The average difference between measurement m i and computation Note that if we had used an approach by topological decomposition of domains, it would have been difficult to take into account the Z t3 interaction. This interaction is of major significance. By not taking it into account the computation differs enormously from the measurement as shown in figure 15.

c i (|m i -c i | /m i ) is 33%

B. Coupling inside cavities

Various models were developed to compute problems in cavities (without including those that use 3D codes) [START_REF] Kasmi | Stochastic Kron's model inspired from the Random Coupling Model[END_REF]. To model confined fields many techniques may be incorporated in Kron's formalism such as chain matrices, Branin's model for guides, RLC circuit networks or transfer functions that behave like random couplings. One of these techniques may be more efficient than the others depending on the problem being investigated. Modelling the cavity itself is not difficult since an EMC engineer naturally looks mainly for order of magnitudes. Estimating the coupling between wires and field inside cavities is much more difficult. Using Kron's method, once the coupling impedances between the field modes and wires are obtained, the problem is entirely solved. For simple geometries, it is quite easy to estimate the analytical formulation in order to compute the coupling functions, and various works show this capability [START_REF] Breant | MKME: Simulation dun système complexe, de la cavité à l'électronique[END_REF]. For more complex cavities, one technique consists in following Collins [START_REF] Collin | Field theory of guided waves[END_REF] to approach cavity geometry using several guided waves. With more than one guided wave and various interfaces, many geometries can be approached sufficiently well to give good results compared to measurements [START_REF] Leman | Contribution à la résolution de problèmes de compatibilité électromagnétique par le formalisme des circuits électriques de Kron[END_REF]. We describe here the "RLC" approach (resonators) which is very efficient, giving accuracy of the same order as when using 3D codes. It requires a great effort to determine the circuit components [START_REF] Demaison | Démonstration de la non équivalence formelle entre les mesures d'émissions rayonnées en CRBM et en chambre anéchoïques, remarques générales sur l'emploi des CRBM[END_REF]. But it is possible to benefit from all the works done previously in high frequency modelling involving cavities (see for example [34] [35]). figure 16 shows our problem made of two small antennas inside a cavity. Each field mode is represented by an "RLC" oscillator. Various elements inside the cavity interact with these modes, and if coupling exists between modes, these oscillators can be coupled with each other.

Each time a coupling occurs, a cord represents it on the graph making an electromagnetic link between the coupled elements. The emitting antenna is a simple loop. This kind of primitive circuit is directly given in the mesh space with an impedance equal to:

z e = R e + L e p + 1 C e p (31)
p being the Laplace operator. The same equation is used for the receiving antenna with index "r". Under the hypothesis of weak losses, the resonance frequencies for a rectangular cavity of length L, height H and width w are given by:

f mnq = c 2 m L 2 + n H 2 + q w 2 (32)
Under the hypothesis that only one polarisation exists in one dimension H, we can obtain the capacitor value for one resonator at one field mode making the equivalence between the electric field energy stored in the volume of the cavity and 1/2CV 2 formulation, V being the field work along the dimension H and C the capacitor we want to calculate. It means that:

0 w 0 L O H 0 dxdydz sin mπx w sin qπy L 2 . . . . . . = C H O dz 2 (33)
After integration we obtain for an mnq mode: 34) From C mnq we can compute L mnq using:

C mnq = 0 2H wL -wL 2mqπ 2 (cos (2qπ) cos (2mπ) . . . . . . -cos (qπ) -cos (mπ) + 2)} (
L mnq ≈ 1 4π 2 f 2 mnq C mnq (35)
Losses are not taken into account in this equation, always under the hypothesis that they are sufficiently weak. These losses have to be calculated, to obtain the resistance values for the resonators. We start from the composite quality coefficient Q. It involves different loss mechanisms such as losses in the walls (Q W ) and losses associated with antennas (Q a ). Q W accounts for the contribution of overlapping orders of modes given by:

Q W = 3v 2δSµ r 1 + 3λ 16 1 L + 1 w + 1 H -1 (36)
where v is the cavity volume, δ the skin depth, λ the wavelength and S the wall surface. For the receiving and emitting antennas, the quality coefficient is:

Q a = L a ω R a (37)
The whole quality coefficient for our problem is obtained through:

1 Q r = 1 Q W + 2 Q a (38)
This gives the resistance for each resonator. Now, the resistance associated with each mode is given by:

R mnq = L mnq 2πf mnq Q r (39)
Once the cavity model is determined (we suppose that models for the antennas are available and given by z a ), the coupling coefficient between all the components involved in the problem remains to be obtained. At this level, the impedance tensor results from the direct summation of all the primitive elements involved in the problem. It means that:

z = ⊕ k z k =     z a 0 . . . 0 0 R mnq + L mnq p + 1 Cmnqp . . . 0 . . . . . . . . . . . . 0 0 . . . z a     (40)

Coupling function calculation

In general, if we consider a coupled system, it is identified by a fundamental impedance tensor given by:

g = a β β b (41)
The β function represents the coupling phenomenon. In the mesh space of mesh currents J k , the system of equations for one source on the first mesh is given by:

e 1 = aj 1 + βj 2 0 = βj 1 + bj 2 (42)
To identify the condition that leads to the maximum energy on the receiving part, we look for:

∂j 2 ∂β = 0 ( 43
)
This gives β = √ ab. For magnetic energy, the coefficient is in general -α √ L 1 L 2 p. The coefficient α accounts for losses in the coupling process. So in general we have β = α √ ab. Determining the coupling function is reduced to the identification of this coefficient α. In the case of a closed circuit of normal n, we consider the coupling with the magnetic field. The coupling coefficient is in this case given by:

α = drdθ n • B 1 (x, y, z) (44)
For example, locating a piece of wire as the emitting antenna at the center of the cavity with its normal following y on x 0 point for a length χ leads to the next coefficient α = α e :

α e = 1 χ χ 0 dxcos mπx w(t) sin π(x -x 0 ) R b (46)
The same calculation has to be made for the receiving antenna. Adding these coupling functions to the impedance tensor gives the complete matrix of this tensor and solves the complete coupled system. The fundamental tensor looks like: α xy are the coupling function between antennas and field modes, or between modes, etc. Note that we accept here that a direct coupling occurs between the antennas through α aa . In [START_REF] Demaison | Démonstration de la non équivalence formelle entre les mesures d'émissions rayonnées en CRBM et en chambre anéchoïques, remarques générales sur l'emploi des CRBM[END_REF], complete exercise was done and gives very good results in comparison with measurements, we describe now this experiment. The experiment consists in locating a small rectangular open end wire antenna on one disk of a cylindrical copper cavity. The cavity is 10 cm in diameter and 30 cm long (figure 17), and we can see on the figure the half-wave antenna centred on the bottom disk of the cavity. The half-wave antenna is 3 cm long, parallel to the metal plate of the disk, at 5 mm above it. Its base is centred and connected to an SMA connector, in order to be power-supplied using a 50Ω high frequency amplifier. The cavity is excited at T EM 010 first mode (cylindrical one). The coupling coefficient between the first antenna and the cavity was computed using the scalar product between the magnetomotive force of the antenna and the mode reluctance:

z =     z a α e . . .
α = l 2 π 1 a+l cos π 2 (a+l) a -1 -. . . 1 a-l cos π 2 (a-l) a -1 (48)
where a is the cavity diameter and l the antenna length. The antenna is defined using:

             Z c = 60log 4h d R a = 1 lσ Cu l π.d.e Z a = Zc p.tan( 2πf l c ) ( 49
)
where h is the antenna height on the plate, σ Cu the copper conductivity, d the antenna diameter and e the skin depth. c is the light speed and f the frequency.

The equivalent circuit for the mode is given by its com-FIG. 18. Measurement versus calculation in cavity ponents:

           L 1 = 1 2πf1 C 1 = L 1 R 1 = L 1 2π f1 Qr ( 50
)
These relations lead to the coupling coefficient β = α √ L a L 1 . So the fundamental tensor for this problem directly written in the mesh space is:

z m = R g + R a + 1 Cap + L a p -β -β R 1 + 1 C1p + L 1 p (51)
figure 18 shows the results. A difference, graphically measured between calculation and measurement for the S 11 parameter obtained at the first antenna input is 10% in frequency and 2% in amplitude.

C. EMC analysis in a dynamic scenario

A set of available analytical equations allows us to make in-depth theoretical studies of problems. In various cases it is possible to demonstrate general results without making any numerical computation or optimisation in the design phase [START_REF] Breant | To improve the variability of one complex system with the MKME[END_REF][37] [38][39]. It would be difficult to write the equation of a system using nodal, direct or other methods. From this point of view, Kron's method is quite unique. In [START_REF] Maurice | La compatibilité électromagnétique des systèmes complexes[END_REF] for example, a simple exercise on general behaviour of a system submitted to external sources was proposed. One example can be given, inspired from Kron's work on electrical machines [START_REF] Kron | Equivalent circuits for the hunting of electrical machinery[END_REF]. We imagine a first electronic device A emitting a magnetic field that can disturb a second one B seen as a receiver. In the first step we look at the coupling relation between the two electronic devices. A turns around B, so that the magnetic interaction between A and B is not only FIG. 19. Interaction between two loops in movement reduced to the mutual inductance M 12 between them, but also includes a cut flux coupling G 12 Ω where G 12 is the magnetic coefficient of the cut flux and Ω the angular speed. We suppose in this simple illustration that all parameters are constant (for example Ω(t) = Ω). figure 19 shows the system analysed.

We write the equations of the system when each electronic device is reduced to a simple mesh and in its original configuration, like a magnetic loop:

   e A = r AA i A + L AA di A dt + M AB di B dt + G AB Ωi B e B = M AB di A dt + r BB i B + L BB di B dt (52)
This system constitutes the primitive system for Kron. Now we make a connection between currents of the same system but for a different angle at any other time:

i A = i A cos (θ(t)) i B = i B (53)
This leads to the connectivity:

C = cos (θ(t)) 0 0 1 (54)
System (50) can be generalized by writing:

e x = r xx i x + L xy di y dt + G xy Ωi y (55)
where L tensor includes M components. Using [START_REF] Casagrande | High frequency bundles modeling[END_REF] we have i x = C x β i β . With replacement in [START_REF] Durand | The Kron formalism of tensor analysis, applied to graphs, networks and antennas[END_REF] this leads to: dt . The ratio r between the electromotive force received in the moving loop compared to the electromotive force received in a static one in the harmonic domain becomes:

e x =
r = r αβ + L αβ p + L αy pC y β Ω 2 + G αβ Ω r αβ + L αβ p (61)
If r αβ is low, G ≈ L, this leads for θ = 0 to:

r = Ω 2 + Ω p + 1 (62)
If Ω << ω, r ≈ 1 + Ω 2 , it shows that the influence increases when the angular speed is squared. The phenomenon may be significant in magnitude. In this case, its contribution should be taken into account for risk prediction. This kind of theoretical analysis can be done thanks to Kron's formalism and linked with tensorial algebra, giving access to complex situations, even relativistic ones [START_REF] Maurice | Introduction dune théorie des jeux dans des topologies dynamiques[END_REF][42] Far beyond EMC, Kron's work has influenced many communities through his approach for systems [START_REF] Happ | Gabriel Kron and Systems Theory[END_REF]. His formalism makes it possible to synthesize general interactions in an easier way than any other, even in multiphysics, since his fundamental tensor leads directly to Lagrange's equations of dissipative systems.

V. CONCLUSION

It is little known that the nodal method usually employed for circuits was replaced by the "modified nodal analysis" in order to take into account driven sources. But in [START_REF] Peikari | Fundamentals of network analysis and synthesis[END_REF], it was demonstrated that this problem can be taken into account in classical nodal analysis, through the mesh space! In fact, Kron's mesh space is able to solve any circuit today. Furthermore the method makes it possible to carry out in-depth analyses of a problem using the geometrical approach in the real mathematical meaning of topology. We have shown through three examples that Kron's method applied to EMC gives the technique to study complex systems using all the theoretical knowledge developed in EMC. It can be the first task done by EMC engineers before using numerical tools to validate their designs. No doubt many other techniques will be found based on Kron's works usable partly for EMC, such as the recent "xTAN" method to take into account the human factor in the studying of complex systems from the EMC point of view [START_REF] Maurice | Introduction dune théorie des jeux dans des topologies dynamiques[END_REF][51][52] [START_REF] Casagrande | High frequency bundles modeling[END_REF].

α a α u α → a α u α ( 1 )

 1 Vectors projected in a base b k are identified by their coordinates: u = u k b k → u k , with the index up. The k component of a vector appears as a superscript. Covectors are noted α u • b k = ω k . The k component of a covector appears as a subscript.

FIG. 1 .

 1 FIG. 1. Euler-Poincaré's relation

) figure 4

 4 figure 4 illustrates this example.

FIG. 4

 4 FIG. 4. cord illustration

FIG. 7 .FIG. 9 .

 79 FIG. 7.A set of primitive circuits in the mesh space (1 st step of Kron's model)

FIG. 11 .

 11 FIG. 11. Shielded cable with internal pigtails

FIG. 13 .

 13 FIG. 13. Pigtail dimensions

  for a standard deviation of 26% (ψ = 45.2, σ = 10 8 , N = 16, C = 120, d = 112.5µm, D = 15mm).

FIG. 14 .

 14 FIG. 14. Comparison computation -measurements

FIG. 16 .

 16 FIG. 16. Interaction inside a cavity: graph

B 1

 1 is a normalized field (here in T •m -3 ) in order to have α value without dimensions reaching 1 maximum. The magnetic field modes b are defined by: b = cos qπy L sin mπx w(t) u x + . . . . . . + sin qπy L cos mπx w(t) u y

  α aa α e R mnq + L mnq p + 1

				FIG. 17. Cavity photograph
				
	. . .	. . .	Cmnqp . . . α ca . . . . . .	   (47)
	α aa	α ca	. . . z a

  r xx C x β i β + L xyNow we can multiply each member on the left by C x α to obtain:e α = r αβ i β + L αy dC y β dθ Ωi β + L αβ di β dt + G αβ Ωi β (59)We note H αβ,θ = L αy dC y β dθ to obtain:e α = r αβ i β + H αβ,θ Ωi β + L αβ di β dt + G αβ Ωi β (60)Now we can wonder what the added term H αβ,θ Ωi β represents when compared to a classic EMC evaluation of coupling between two loops in static configuration. For a static system, the equation is: e a = r ab i b + L ab di b

			dC y β i β dt	+ G xy ΩC y β i β	(56)
	But:
	dC x β i β dt	=	dC x β dθ	dθ dt	i β + C x β	di β dt	(57)
	Knowing that θ = Ω we obtain:
	e dC y β dθ	Ωi β + C y β	di β dt	+ G xy ΩC y β i β
								(58)

x = r xx C x β i β + L xy

ACKNOWLEDGMENTS

The author would like to specially thank Philippe Besnier who greatly helped him to present his work on Kron's method.

Appendix A: Kron's diakoptics

Diakoptics is perhaps the best known technique developed by Kron[54] [START_REF] Happ | Diakoptics and networks[END_REF]. There are various ways to apply diakoptics. We consider here a general principle among the simplest ones.

From a first set of solutions, we obtain mesh current values depending on sources:

-1 E m . Then we add loads on the boundaries of the previous network. This can be written:

because loads do not have self sources and i f are the boundary currents. Which leads to the system

The new currents with added loads can be obtained from the previous circuit computing:

Yuri Sohor has developed very efficient techniques to use diakoptics for large problems [START_REF] Sohor | Application of Tensor Decomposition in Studying Linear Induction Motors[END_REF].