Louis Guillot

Bruno Cochelin

Christophe Vergez

A generic and efficient Taylor series based continuation method using a quadratic recast of smooth nonlinear systems

Keywords: continuation, Taylor series, quadratic recast, asymptotic numerical method, nonlinear systems, finite element method

This paper is concerned with a Taylor series based continuation algorithm, ie, the so-called Asymptotic Numerical Method (ANM). It describes a generic continuation procedure that apply the ANM principle at best, that is to say, that presents a high level of genericity without paying the price of this genericity by low computing performances. The way to quadratically recast a system of equation is now part of the method itself, and the way to handle elementary transcendental function is detailed with great attention. A sparse tensorial formalism is introduced for the internal representation of the system, which, when combines with a block condensation technique, provides a good computational efficiency of the ANM. Three examples are developed to show the performance and the versatility of the implementation of the continuation tool. Its robustness and its accuracy are explored. Finally, the potentiality of this method for complex non linear finite element analysis is enlightened by treating 2D elasticity problem with geometrical nonlinearities.

Introduction

The so-called Asymptotic Numerical Method (ANM), first described in [START_REF] Cochelin | Asymptotic-numerical methods and pade approximants for non-linear elastic structures[END_REF] and [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF], is a continuation technique based on high order Taylor series expansions of the unknowns with respect to a path parameter. The solution branches are computed step by step as in a classical predictor-corrector algorithm (see [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF], [START_REF] Doedel | Numerical analysis and control of bifurcation problems (i): bifurcation in finite dimensions[END_REF] and [START_REF] Seydel | From equilibrium to chaos: practical bifurcation and stability analysis[END_REF]), but as opposed to first order prediction, the accuracy of a high order Taylor series prediction is so high that no correction is needed in general. Each step provides a local continuous representation of the branch whose length is computed afterward directly from the convergence properties of the series. Having a continuous description of the branch makes the continuation very robust and provides some opportunities in the detection of bifurcations [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF].

ANM has first been applied to compute equilibrium branches of geometrically nonlinear finite elements models in structural mechanics [START_REF] Cochelin | The asymptotic-numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF] [33] 1 . It has been later extended to more complex and less regular non linearities such as the one associated to material non linearities [START_REF] Zahrouni | Asymptotic numerical method for nonlinear constitutive laws[END_REF], [START_REF] Abichou | Asymptotic numerical method for problems coupling several nonlinearities[END_REF], contact conditions and friction [START_REF] Kessab | High order predictor-corrector algorithms for strongly nonlinear problems[END_REF], [START_REF] Aggoune | Asymptotic numerical methods for unilateral contact[END_REF], vibrations of viscoelastic shells [START_REF] Duigou | Iterative algorithms for nonlinear eigenvalue problems. application to vibrations of viscoelastic shells[END_REF] and material instabilities [START_REF] Nezamabadi | A multilevel computational strategy for handling microscopic and macroscopic instabilities[END_REF]. Stationary solution of Navier-Stokes flow [START_REF] Allery | Application of the asymptotic numerical method to the coanda effect study[END_REF], [START_REF] Medale | High performance computations of steadystate bifurcations in 3d incompressible fluid flows by asymptotic numerical method[END_REF] has also been adressed, showing that very large scale problems with more than a million degree of freedom can be treated efficiently with this method. The continuation of periodic solutions of dynamical systems is presented in [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF], [START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF] where the ANM is combined with a Fourier series expansion known as the harmonic balance method. Extension to the continuation of quasi-periodic solutions is adressed in [START_REF] Guillot | Continuation of quasiperiodic solutions with two-frequency harmonic balance method[END_REF]. The continuation solver MANLAB-1.0 [START_REF]Manlab -an interactive path-following and bifurcation alnalysis software[END_REF] has been the first attempt to design a general purpose continuation software working on the ANM principle. It allows a user to enter its own algebraic system and to interactively draw the bifurcation diagram. It works well for systems with a few hundreds of equations but the quadratic recast of elementary transcendental functions was a little abstruse and finite element mechanical models could not be implemented because of the lack of performances.

A common denominator of the works mentioned above is the requirement to recast the governing equations in a quadratic format. The motivation for this recast is that the Taylor series computation becomes easy and efficient when a system is quadratic. However, this requirement has turn out to be the most difficult point of the method. Its misunderstanding has prevented many potential users and developpers from succeeding in the use of the method, especially when the system of equation contains elementary transcendental functions (exp, log, sin, cos, power, ...). Introducing Automatic Differentiation (AD) principle inside the ANM process [START_REF] Charpentier | The diamant approach for an efficient automatic differentiation of the asymptotic numerical method[END_REF], [START_REF] Charpentier | Différentiation automatique de la méthode asymptotique numérique typée: l'approche diamant[END_REF] has been the first elegant way of removing this difficulty. The DiaManlab software [START_REF] Charpentier | Towards a full higher order ad-based continuation and bifurcation framework[END_REF] allows the users to write the governing equations in a natural way, without any need of a quadratic recast. However, once again, the price to pay for this genericity is a reduction of the computing performance that may be again a drawback for large systems resulting for instance from the finite element (FE) discretization of a continuous problem.

In this paper, combinig the experiences of twenty years of work on the ANM and this experience on Automatic Differentiation, we present a generic way to solve an algebraic system with the ANM principle applied at best, ie, providing both genericity and efficiency. The quadratic recast of the original system is now a part of the method itself instead of being a pre-processing task left to the user. The declaration of a list of auxiliary variables that follow the linear declaration rule provides a clear and useful distinction between the main and the auxiliary variables, and between the main and the auxiliary equations. Attention is paid to the treatment of elementary transcendental functions [START_REF] Zucker | Elementary transcendental functions: Logarithmic, exponential, circular and hyperbolic functions[END_REF] and the way their (quadratic) differential form should appear in the process. Additionally, this results in a clear separation between quadratic and functional equations. Another important improvement is the use of a sparse tensor formalism for the internal representation of a quadratic system on a computer. This allows an efficient construction of the linear systems that have to be solved for the series computation, including a block inversion for the condensation of the auxiliary variables. Another key point is that the sparse tensors are automatically constructed from the equations by using polarization formula. This allows to write the quadratic equations in a quite natural way.

The generic method described here has been implemented in Matlab as the fourth version of the Manlab suite. It permits to deal with a very large class of algebraic systems containing many elementary transcendental functions as it will be shown on three examples in section 5. Another very important potentiality of this method is the capacity to implement a FE mechanical model in a quite simple way and, that time, without suffering from poor computing performances. Only a quadratic writing of the governing equation, ie, of the residual vector equation, is required.

2 Continuation using the Asymptotic Numerical Method (ANM)

Definitions, notations and theoretical background

Consider the algebraic system

R(u, λ) = 0 (1)
where R : R Neq+1 → R Neq is a function called the residue, u ∈ R Neq is the vector of unknowns and λ ∈ R is a parameter of interest. The notation U = (u, λ) ∈ R Neq+1 is sometimes used in the following parts. R is a real analytic function and can be written as a sum, product, composition of rational and elementary transcendental functions such as polynomials, exp, cos, log, etc... The aim of continuation is to find the solution set of (1). This is a union of various solution branches that possibly cross at bifurcation points, which is called the bifurcation diagram. The procedure to compute a continuation step is now explained. Let U 0 = (u 0 , λ 0 ) be a regular solution of equation ( 1). This means that the rank of the Jacobian matrix at point U 0 ,

∂R ∂U = ∂R ∂u , ∂R ∂λ (2)
is N eq . Let U 0 be the starting point. Let U 1 = (u 1 , λ 1 ) be a unitary tangent vector at point U 0 and a = (u

-u 0 ) • u 1 + (λ -λ 0 )λ 1 = (U -U 0 ) • U 1 be the pseudo arc-length parameter introduced in [6]
, [START_REF] Doedel | Numerical computation of periodic solution branches and oscillatory dynamics of the stirred tank reactor with a yields b yields c reactions[END_REF]. Note that the definition of a is used as a closing equation for the system (1). The analytic implicit function theorem (which complex form can be found in [START_REF] Fritzsche | From holomorphic functions to complex manifolds[END_REF]) allows to search the solution branch of equation ( 1) around U 0 as an analytic functions of a, i.e as Taylor series expansion with respect to a:

U (a) = U 0 + aU 1 + a 2 U 2 + a 3 U 3 + . . . (3)
Note that in (3) λ is also developed as a series of the path parameter a. In practice, the series are truncated at an order p. They are then introduced in equation ( 1) which gives the Taylor series development of R(U (a)) :

R(U (a)) = R(U (0)) + aR 1 + a 2 R 2 + a 3 R 3 + . . . (4)
where

R 1 = dR da a=0 = ∂R ∂U U 1 R 2 = 1 2 d 2 R da 2 a=0 = ∂R ∂U U 2 -F 2 (U 1 ) R 3 = 1 3! d 3 R da 3 a=0 = ∂R ∂U U 3 -F 3 (U 1 , U 2 ) . . . . . . . . . R p = 1 p! d p R da p a=0 = ∂R ∂U U p -F p (U 1 , . . . , U p-1 ) (5)
with F k being functions that depend only on already computed terms of the series.

The system (1) then looks like : ∀k ≤ p, R k = 0. It has to be noticed that all these equations share the same matrix to inverse that is the Jacobian matrix (2) at point U 0 . When the series (3) is computed up to order p, the domain of utility [0, a max ] is determined. The maximum increase of the residue (1) is set to a small value ε 1 ,

such that |R(U (a)) -R(U (0))| < ε 1 for all a ∈ [0, a max ]. With the approximation R(U (a)) -R(U (0)) = a p+1 R p+1 , it gives a max = ( ε 1 ||R p+1 || ) 1 p+1 (6)
Finally, the ending point of the continuation step U (a max ) is computed and it is used as the starting point of the next continuation step. This approach leads to a continuous representation of the solution branch. There are almost no corrections needed to stay on the branch in practice because the Taylor series development (3) is a very accurate approximation of the true solution branch on its domain of utility. This domain of utility can be computed directly from the knowledge of the series with the formula (6).

Implementation methods

To summarize, the Asymptotic Numerical Method (ANM) is a method to compute a solution branch from the knowledge of a starting point U 0 , of a way to compute the Jacobian matrix (2) of the system and of a procedure to compute the Taylor series (3) i.e. to compute U k from U j , j < k. As explained in the Introduction, there are several ways to deal with these requirements. One of the solution is the use of automatic differentiation. In this case, the knowledge of equation [START_REF]Manlab -an interactive path-following and bifurcation alnalysis software[END_REF] only is enough to compute the bifurcation diagram. It has been done in [START_REF] Charpentier | Diamanlab-an interactive taylor-based continuation tool in matlab[END_REF] in an elegant way. The set up is surprisingly simple for the user, however it becomes computationally expensive when the size of the system increases above several hundreds of unknowns.

Here we follow the traditional idea of the ANM that is described in the book [START_REF] Cochelin | Méthode asymptotique numérique[END_REF] and the papers [START_REF] Cochelin | A path-following technique via an asymptotic-numerical method[END_REF], [START_REF] Cochelin | Asymptotic-numerical methods and pade approximants for non-linear elastic structures[END_REF] and [START_REF] Cochelin | The asymptotic-numerical method: an efficient perturbation technique for nonlinear structural mechanics[END_REF]. The system is written in a quadratic format that allows direct computation of both the jacobian matrix (2) and of the coefficients of the Taylor series (3).

Recast of the system of equations

In this section the equations (1) are rewritten in a quadratic format in order to apply the ANM. The general set up follows an explanatory example.

Explanatory example

Let us consider an academic example containing two elementary transcendental functions, exp and tanh for which the residue function ( 1) is given by R = (r 1 , r 2 ) :

r 1 (u 1 , u 2 , λ) = u 1 + λ exp(u 2 ) 1+u 1 r 2 (u 1 , u 2 , λ) = u 2 + u 1 tanh( -5u 1 1+u 1 u 2 ) (7)

Definition of the auxiliary variables

The equations [START_REF] Charpentier | Towards a full higher order ad-based continuation and bifurcation framework[END_REF] are rewritten hereafter in a quadratic format. To achieve this, auxiliary variables are introduced. These additional variables are defined from the main variables u 1 , u 2 , λ.

To compute r 1 from given input values (u 1 , u 2 , λ), one first decomposes the process into elementary operations and stores intermediate results. For instance, an evaluation of exp(u 2 ) is stored. Then, an evaluation of 1 + u 1 is stored as well. After, the first result is divided by the second one. This number is multiplied by λ, etc... Following this, let v 1 , v 2 , v 3 , v 4 , v 5 be:

v 1 = exp(u 2 ) v 2 = v 1 1+u 1 v 3 = 1 + u 1 u 2 v 4 = -5u 1 v 3 v 5 = tanh(v 4 ) (8)
From these definitions we can now write the system [START_REF] Charpentier | Towards a full higher order ad-based continuation and bifurcation framework[END_REF] with only quadratic nonlinearities :

r 1 = u 1 + λv 2 r 2 = u 2 + u 1 v 5 (9)
In the definition of auxiliary variables [START_REF] Charpentier | Diamanlab-an interactive taylor-based continuation tool in matlab[END_REF] there are different types of equations. Some are already quadratic, here v 3 = 1+u 1 u 2 . Some are easily made quadratic, here

v 2 = v 1 1+u 1 and v 4 = -5u 1 v 3 are respectively rewritten v 2 (1 + u 1 ) = v 1 and v 4 v 3 = -5u 1 which are equivalent to the first definitions if u 1 = -1 and v 3 = 0.
The last type of auxiliary variables, here v 1 = exp(u 2 ) and v 5 = tanh(v 4 ) will be treated in the next section.

The auxiliary variables are defined by following the Linear declaration rule i.e. only from the main variables and from auxiliary variables with a smaller index. To write it in a compact way : ∀i,

v i = f (u, λ, v 1 , . . . , v i-1
). It ensures that r 1 , r 2 can be computed from the input u 1 , u 2 , λ. The choice of the auxiliary variables is not unique. For example, the definitions of v 3 and v 4 could be replace by v 3 = u 1 u 2 and v 4 = -5u 1 1+v 3 respectively. Of course, this does not change the final result. The linear declaration rule ensures that the system of equations ( 7) is equivalent to the two systems ( 8) and (9).

Quadratic recast of the elementary transcendental functions

The case of v 1 = exp(u 2 ) and v 5 = tanh(v 4 ) is treated here. To achieve their quadratic recast, these equations are first differentiated:

dv 1 = exp(u 2 )du 2 = v 1 du 2 dv 5 = (1 + tanh 2 (v 4 ))dv 4 = (1 + v 2 5 )dv 4 (10)
Then the auxiliary variable v 6 = (1+v 2 5 ) is defined to obtain the final quadratic recast of the differentiated equations

dv 1 = v 1 du 2 dv 5 = v 6 dv 4 v 6 = 1 + v 2 5 ( 11
)
This formulation is said to be quadratic because it is quadratic with respect to (u, v, λ, du, dv). A table of recast of the most common transcendental functions is available in Appendix A.

Final recast of the example

The unknowns are now split in two families : U = (u 1 , u 2 , λ) for the original main unknowns and

U aux = (v 1 , v 2 , v 3 , v 4 , v 5 , v 6 )
for the auxiliary unknowns. The full vector of unknowns is called U tot = (U, U aux ). From the previous section, the system (7) has been written

R tot (U tot ) =             u 1 + λv 2 u 2 + u 1 v 5 v 1 -exp(u 2 ) v 1 -v 2 -v 2 u 1 1 + u 1 u 2 -v 3 -5u 1 -v 4 v 3 v 5 -tanh(v 4 ) 1 + v 2 5 -v 6             =             R R aux             (12)
where the two first equations correspond to the initial system of equations ( 7) and the other ones to the auxiliary variables. The horizontal line materializes this separation : R tot = (R, R aux ). There is another useful separation of the residue between already quadratic equations, R quad , and the other ones, R fun where an elementary function appears. R tot is split up in these two parts :

R tot (U tot ) =             u 1 + λv 2 u 2 + u 1 v 5 0 v 1 -v 2 -v 2 u 1 1 + u 1 u 2 -v 3 -5u 1 -v 4 v 3 0 1 + v 2 5 -v 6             +             0 0 v 1 -exp(u 2 ) 0 0 0 v 5 -tanh(v 4 ) 0             = R quad + R fun (13)
To this "total" residue is added the differentiated form of its third and seventh components, that is of R fun , in dR fun vector,

dR fun =             0 0 dv 1 -v 1 du 2 0 0 0 dv 5 -v 6 dv 4 0             (14)
The differentiated form of the already quadratic equations does not need to be specified.

General formula

From a general point of view, the residue (1) R can be written in the following way:

R tot (U tot ) = R quad (U tot ) + R fun (U tot ) = [C + L(U tot ) + Q(U tot , U tot )] + [L d (U tot ) -f (U tot )] (15)
where U tot = (u, λ, v) is the vector of all the unknowns, v are the auxiliary variables. R quad is the already quadratic part of the residue and R fun is the part which is not quadratic (but its differential form is). C is a constant operator, L and L d are linear operators, Q is a quadratic operator and each component of f (U tot ) is an elementary transcendental function of U tot or zero. f is a non-quadratic function whose differential form can be written quadratically. Formally it means that the application (U tot , dU tot ) → d Utot f (dU tot ), where d Utot f is the differential application of f at point U tot , is bi-linear. This bi-linear application is called

Q d . One has Q d (U tot , dU tot ) = d Utot f (dU tot ). To the system (15) is added the differentiated form of R fun (U tot ) = L d (U tot ) -f (U tot ) : d Utot R fun (dU tot ) = L d (dU tot ) -Q d (U tot , dU tot ) (16)
This is used in the next section to compute the Jacobian matrix (2) of the system and the Taylor series (3).

Details on the computation of the series

The formulations ( 15) and ( 16) and only these, are used for the following computations.

Computation of the Jacobian matrix

General computation

From equation ( 15), the Jacobian matrix can be computed automatically :

∂Rtot ∂Utot = L + Q(U tot , •) + Q(•, U tot ) + L d -∂f ∂Utot (•) = L + Q(U tot , •) + Q(•, U tot ) + L d -Q d (U tot , •) (17)
In our example the only non-zero elements of Q d are exp(u 1 )du 1 and (1-tanh(v 4 ) 2 )dv 4 .

They can be written v 1 du 1 and v 6 dv 4 respectively, as in [START_REF] Cochelin | Méthode asymptotique numérique[END_REF].

The Linear declaration rule ensures that the right-bottomed sub-matrix ∂Raux ∂Uaux that contains the partial derivatives of the equations defining the auxiliary variables with respect to the auxiliary variables is triangular (see example below). This information is used for efficient block-solving in the computation of the series.

Explanatory example

The jacobian matrix is the matrix of partial derivatives of the full residue with respect to all the variables. The jacobian matrix of the example is computed with the equations ( 13) and [START_REF] Cochelin | Méthode asymptotique numérique[END_REF].

J Utot = ∂R ∂U ∂R ∂Uaux ∂Raux ∂U ∂Raux ∂Uaux =             1 0 v 2 0 λ 0 0 0 0 v 5 1 0 0 0 0 0 u 1 0 0 -v 1 0 1 0 0 0 0 0 -v 2 0 0 1 -1 -u 1 0 0 0 0 u 2 u 1 0 0 0 -1 0 0 0 -5 0 0 0 0 -v 4 -v 3 0 0 0 0 0 0 0 0 -v 6 1 0 0 0 0 0 0 0 0 2v 5 -1             (18)

Computation of the series

In this section, the subscripts "tot" are omitted : U tot is simply written U . The focus is on the general case. The equations ( 15) are written with an already quadratic part R quad and a functional part R fun . For the quadratic part of equation [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF], U is replaced by its series and the expansions are truncated at order N :

0 = C + L N k=0 a k U k + Q N k=0 a k U k , N k=0 a k U k = C + N k=0 a k LU k + i+j=k Q(U i , U j ) + O(a N +1 ) (19)
Equating the coefficients of the series for each power of a to 0,

C + LU 0 + Q(U 0 , U 0 ) = 0 LU 1 + Q(U 0 , U 1 ) + Q(U 1 , U 0 ) = 0 LU 2 + Q(U 0 , U 2 ) + Q(U 2 , U 0 ) = -Q(U 1 , U 1 )
. . .

LU k + Q(U 0 , U k ) + Q(U k , U 0 ) = -k-1 i=1 Q(U i , U k-i ) (20)
The first equation is just the definition of a solution point U 0 and the second one the definition of the tangent vector

U 1 at U 0 . LU k + Q(U 0 , U k ) + Q(U k , U 0 ) is the first part of the jacobian matrix ∂R quad ∂U
at point U 0 applied at point U k as shown in equation [START_REF] Bathe | Finite element procedures[END_REF].

For the non-quadratic part of equation ( 15), the differentiated form is used :

L d (dU ) -Q d (U, dU ) = 0 (21)
Replacing U and dU by their series expansion one gets

0 = L d ( ∂U ∂a ) -Q d U, ∂U ∂a = N -1 k=0 (k + 1)a k L d U k+1 -Q d N k=0 a k U k , N -1 k=0 (k + 1)a k U k+1 = N -1 k=0 a k (k + 1)L d U k+1 -i+j=k Q d (U i , (j + 1)U j+1 ) + O(a N ) (22)
Equating the coefficients of the series to 0,

1L d U 1 -Q d (U 0 , 1U 1 ) = 0 2L d U 2 -Q d (U 0 , 2U 2 ) = Q d (U 1 , U 1 )
. . .

kL d U k -Q d (U 0 , kU k ) = k-1 i=1 Q d (U i , (k -i)U k-i ) (23)
Or dividing by k :

L d U 1 -Q d (U 0 , U 1 ) = 0 L d U 2 -Q d (U 0 , U 2 ) = 1 2 Q d (U 1 , U 1 ) . . . L d U k -Q d (U 0 , U k ) = 1 k k-1 i=1 (k -i)Q d (U i , U k-i ) (24) L d U k -Q d (U 0 , U k ) is the second part of the jacobian matrix ∂R fun
∂U at point U 0 applied at point U k as shown in equation [START_REF] Bathe | Finite element procedures[END_REF].

Calling J U 0 the jacobian matrix of R at point U 0 and putting together the preceding computations, the system to solve at order p ≥ 2 is deduced :

J U 0 (U p ) = p-1 i=1 -Q(U i , U p-i ) + p -i p Q d (U i , U p-i ) := F nl p ( 25
)
These systems are solved efficiently using a sparse tensorial formalism and a fast block-solving described hereafter.

Block solving of the linear systems

In the two previous subsections, the computation of the Jacobian matrix (2) of a system of type [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF] was explained. The computation of the Taylor series (3) of a point solution is also explained. Once U k , k ≤ p -1 are known, the formula ( 25) is used to compute U p :

U p = J -1 U 0 F nl p ( 26
)
The focus is on the inversion of the jacobian matrix as it is the most expensive operation in terms of time computation. The system has been written thanks to auxiliary variables and is written separating the auxiliary part from the main part:

R tot (U tot ) := R(U tot ) R aux (U tot ) (27)
with

U tot = (U, U aux ) t , U = (u, λ) ∈ R Neq+1 and U aux ∈ R Neqaux .
The Jacobian matrix of R tot is computed by block :

J U 0 = ∂R ∂U ∂R ∂Uaux ∂Raux ∂U ∂Raux ∂Uaux := K B C K aux (28)
The Linear declaration rule of the auxiliary variables ensures that K aux is triangular, hence it is very easy to inverse. This information is used to solve the linear systems (26) efficiently:

J U 0 tot (U p,tot ) = F nl p,tot ⇔ K B C K aux U p U p,aux = F nl p F nl p,aux (29)
which is equivalent to the system:

KU p + BU p,aux = F nl p CU p + K aux U p,aux = F nl p,aux (30)
Now the second line is multiplied by the inverse of K aux and U p,aux is replaced in the first line accordingly:

K -BK -1 aux C U p = F nl p -BK -1 aux F nl p,aux U p,aux = K -1 aux F nl p,aux -CU p (31)
Finally, there are two systems to solve : one of size N eq + 1 and a triangular one of size N eqaux , instead of one system of size N eq + N eqaux + 1. It goes without saying that this manipulation improves the time computation a lot. The Jacobian matrix J of the main system, ie, without auxiliary variables, is computed with the condensation formula : J = K -BK -1 aux C . In our example, the jacobian matrix [START_REF] Doedel | Numerical computation of periodic solution branches and oscillatory dynamics of the stirred tank reactor with a yields b yields c reactions[END_REF] is divided in four blocks as in equation ( 28)

K = 1 0 v 2 v 5 1 0 B = 0 λ 0 0 0 0 0 0 0 0 u 1 0 C =         0 -v 1 0 -v 2 0 0 u 2 u 1 0 -5 0 0 0 0 0 0 0 0         K aux =         1 0 0 0 0 0 1 -1 -u 1 0 0 0 0 0 0 -1 0 0 0 0 0 -v 4 -v 3 0 0 0 0 0 -v 6 1 0 0 0 0 0 2v 5 -1         (32)
As expected, K aux is triangular. The condensation formula gives :

J = K -BK -1 aux C = 1 -λv 2 u 1 +1 λv 1 u 1 +1 v 2 v 5 + 10u 1 v 5 v 6 v 3 + 2u 1 u 2 v 4 v 5 v 6 v 3 2v 4 v 5 v 6 u 2 1 v 3 + 1 0 (33)
Replacing the auxiliary variables by their values [START_REF] Charpentier | Diamanlab-an interactive taylor-based continuation tool in matlab[END_REF], one could check that J is the jacobian matrix of the original system [START_REF] Charpentier | Towards a full higher order ad-based continuation and bifurcation framework[END_REF]. For clarity, this computation is not shown here.

Sparse tensorial formalism

Sparse tensors

The constant, linear and quadratic operators are written using a sparse tensorial formalism. For R(X) = C + L(X) + Q(X, X) ∈ R N with X ∈ R M , we write :

R i = C i + M j=1 L ij X j + M j,k=1 Q ijk X j X k , 1 ≤ i ≤ N. ( 34
)
Constant sparse tensor For C ∈ R N , a constant operator, the tensorial representation is simply : C i , 1 ≤ i ≤ N . Its sparsity leads to a representation with two lists : the list of indexes iC = {i|C i = 0} of the positions of non-zero coefficients and the list of associated values vC = {C i |i ∈ iC}.

Linear sparse tensor For L : R M → R N , a linear operator, the tensorial repre-

sentation is L ij , 1 ≤ i ≤ N, 1 ≤ j ≤ M .
Its sparsity leads to a representation with three lists : the list of indexes iL and the list of associated variables jL that are defined by iL × jL = {(i, j)|L ij = 0} with associated values vL = {L ij |(i, j) ∈ iL × jL}.

Quadratic sparse tensor For

Q : R M × R M → R N , a quadratic operator, the tensorial representation is Q ijk , 1 ≤ i ≤ N, 1 ≤ j, k ≤ M .
Its sparsity leads to a representation with four lists : the list of indexes iQ and the two list of associated variables jQ and kQ that are defined by iQ × jQ × kQ = {(i, j, k)|Q ijk = 0} with associated values vQ = {Q ijk |(i, j, k) ∈ iQ × jQ × kQ}.

Explanatory example

The sparse tensors of the example are now given with their lists of indexes. The vector of unknowns is here

U tot = (u 1 , u 2 , λ, v 1 , v 2 , v 3 , v 4 , v 5 , v 6 ) of size M = 9.
The N = 8 equations are written in [START_REF] Cochelin | Asymptotic-numerical methods and pade approximants for non-linear elastic structures[END_REF]. The lists are

iC = [5 8] vC = [1 1] (35) iL = [1 2 4 4 5 6 8] jL = [1 2 4 5 6 7 9] vL = [1 1 1 -1 -1 -5 -1] (36) iLd = [3 7] jLd = [4 8] vLd = [1 1] (37) iQ = [1 2 4 5 6 8] jQ = [3 1 1 1 6 8] kQ = [5 8 5 2 7 8] vQ = [1 1 -1 1 -1 -1] (38)
The index iQ represents the line of the equation, jQ the number of a variable, kQ the number of a variable also and vQ the associated value. For example here one can read that on the first line appears the product between the third variable (U (3) = λ) and the fifth variable (U (5) = v 2 ) with associated value 1. Indeed, one can read the term 1 × λ × v 2 on the first line of the total residue vector [START_REF] Cochelin | Asymptotic-numerical methods and pade approximants for non-linear elastic structures[END_REF].

iQd = [ 3 7] jQd = [ 4 9] kQd = [ 2 7] vQd = [-1 -1] (39)
The index jQd corresponds to the variable which is not differentiated and kQd to the variable which is differentiated. Here one can read : On the third equation, there is the product of the fourth variable (U (4) = v 1 ) with the differential of the second variable (dU (2) = du 2 ) with associated value -1. Indeed, the third equation of ( 14) is the differentiated form of the definition of the first auxiliary variable that looks 0 = dv 1 -v 1 du 2 .

Automatic generation of the tensor lists using polarization formula

The construction of the tensor lists from the equations can be done by hand for this simple example, but it is evident that an automated process has to be used to deal with large systems. The key point in this list generation is the separation between the constant, the linear and the quadratic terms. This can be achieved by using the polarization formula as follows : C is computed by the evaluation of R(0) = C + L(0) + Q(0, 0) = C. The vector L(X) is computed by the evaluation of

R(X) -R(-X) = C + L(X) + Q(X, X) -C -L(-X) -Q(-X, -X) = 2L(X).
(40)

The vector Q(X, Y ) is computed by the evaluation of

R(X + Y ) -R(X -Y ) = C + L(X + Y ) + Q(X + Y, X + Y ) -(C + L(X -Y ) + Q(X -Y, X -Y )) = 2L(Y ) + 2Q(X, Y ) + 2Q(Y, X) = 2L(Y ) + 4Q(X, Y ), if Q is symmetric. (41) With the formula R(Y ) -R(-Y ) = 2L(Y ), one gets : C = R(0) L(X) = 1 2 (R(X) -R(-X)) Q(X, Y ) = 1 4 (R(X + Y ) -R(X -Y ) -(R(Y ) -R(-Y ))) (42)
The lists are then constructed by applying the polarization formula (42) to all the vector of the canonical base. This process is automatic from the knowledge of the residue function [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF].

Examples

The objective of this section is twofold. Firstly, to illustrate the quadratic recast on three more examples, some with elementary transcendental functions. Secondly, to illustrate the robustness of the continuation and the simple bifurcation detection. The input files that a user has to provide for running these examples in Manlab-4.0 are not given for brevity. The input files of the explanatory example, corresponding to ( 13) and ( 14), can be found on the web site Manlab at http://manlab.lma.cnrs-mrs.fr/.

Logistic map

The first example is the well-known logistic map. The sequence (x n ) n∈N is defined by the following recurrence relation :

x n+1 = µx n (1 -x n ) := f µ (x n ) (43) 0 ≤ µ ≤ 4
is the bifurcation parameter. Depending on its value, the sequence (x n ) n∈N converges towards a fixed point, a periodic orbit or its behaviour is chaotic. The logistic map undergoes an infinite number of period doubling bifurcation when µ is increased up to a critical value µ crit ≃ 3.56994567 above which its behaviour is chaotic [START_REF] Li | Period three implies chaos[END_REF]. These bifurcations are detected using the method described in [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF].

2 N -periodic orbits are sought for with N ∈ N ; in other words, fixed points of f 2 N µ :

x 1 = f 2 N µ (x 1 ) (44)
As this equation ( 44) is not quadratic, 2 N -1 auxiliary variables that are the iterates of x 1 by f µ along the orbit are introduced : µ is defined so that the equations are recast quadratically as follows :

∀n ∈ {1, . . . , 2 N -1}, x n+1 = f µ (x n ), (45)
λx 1 = x 2 N (1 -x 2 N ) λµ = 1 λx 2 = x 1 (1 -x 1 ) λx 3 = x 2 (1 -x 2 ) . . . . . . λx 2 N = x 2 N -1 (1 -x 2 N -1 ) (46)
The horizontal lines materializes the separation between the main equation and the auxiliary equations. Feigenbaum [START_REF] Feigenbaum | The universal metric properties of nonlinear transformations[END_REF] showed that the ratio of consecutive intervals between two period doubling converges toward a constant δ now called Feigenbaum's constant. This convergence is shown on figure 2. The last two bifurcations are closer than the machine precision as can be seen on figure 1. It is probable that the last bifurcation position was not computed with enough accuracy. This explains the position of the last point on figure 2.

This example is interesting since when N is set to 12 as it is done here, there is only one main equation (44) and 2 N = 4096 auxiliary variables (46). The inversion of a sparse triangular matrix of size 4096 is immediate and then the system to solve becomes of size 1.

Layne-Watson

The Layne-Watson system of equations has been first introduced in [START_REF] Watson | A globally convergent algorithm for computing fixed points of c2 maps[END_REF]. The problem is to find the fixed point of a function g with a homotopy technique. This example has been chosen as it is very tough for continuation algorithms. There are no branch crossing but the bifurcation diagram is very intricate as can be seen on the last plot of figure 3. For N ∈ N and x = (x 1 , . . . , x N ) ∈ R N , g is defined by ∀i ∈ {1, . . . , N } , g i (x) := exp cos i

N k=1 x k (47)
To find the fixed point, the system

R(x, λ) = x -λg(x) (48)
is introduced. Obviously R(0, 0) = 0. This point is used as a starting point of the continuation. If λ = 1 is reached throughout the continuation then a fixed point of g is known. The recast of the equations is, for i ∈ {1, . . . , N },

x i -λg i = 0 c i -cos i N k=1 x k = 0 s i -sin i N k=1 x k = 0 g i -exp(c i ) = 0 (49)
And the differentiated form of the last three vectorial equations, for i ∈ {1, . . . , N },

dc i + s i d i N k=1 x k = 0 ds i -c i d i N k=1 x k = 0 dg i -g i dc i = 0 (50)
For N = 10, 11 fixed points were find. For N = 50, 4000 continuation steps were needed to find 4 fixed points. The computation time was about 45 seconds on a laptop computer (two processors 2,5 GHz Intel Core i7 with 16 Go RAM). In [START_REF] Nedialkov | Solving differential-algebraic equations by taylor series (iii): the daets code[END_REF], the last value tried was N = 130 and some issues appeared while dropping the tolerance under 10 -11 . This issue seems to persist here but above N = 150 approximately. For N = 200, a fixed point of g was found at the tolerance 10 -11 , after more than 75000 continuation steps i.e. 35 minutes of computation, but it was not possible to find a fixed point at a smaller tolerance.

Stick-slip motion of a one d.o.f. oscillator

Let consider a one d.o.f. damped spring-mass oscillator which is driven by a belt with a constant velocity v b . Its position u(t) is governed by the (dimensionless) Newton law

ü + 2η u + u = µ(v) (51)
where η is the damping ratio, µ(v) the friction force and v = u -v b the sliding velocity. The friction force is assumed to follow a Coulomb law with constant static and dynamic friction ratio denoted by µ S and µ d . Precisely, -µ s < µ(v) < µ s when v = 0 (stick), µ(v) = -µ d when v > 0 and µ(v) = µ d when v < 0 (slip). In the following, we look for a periodic solution made of a single stick phase and a single slip phase. Without loss of generality, the origin of time is taken at the transition between stick and slip, where u(0) = µ s -2ηv b and u(0) = v b . The slip phase occurs for t ∈ [0, t 1 ] with µ(v) = µ d . The solution to (51) with above mentioned initial conditions is

u g (t) = µ d + e -ηt A cos( 1 -η 2 t) + B sin( 1 -η 2 t) (52)
with

A = µ s -µ d -2 η v b and B = v b +ηA √ 1-η 2 .
The stick phase occurs for t ∈ [-t 2 , 0] with u ′′ (t) = 0, and the solution is

u s (t) = (µ s -2ηv b ) + v b t (53)
The main equation are the following periodicity conditions involving the two unknowns t 1 , t 2 and the two free parameters η and v b (µ s and µ d are constant).

u g (t 1 ) = u s (t 2 ) ug (t 1 ) = us (t 2 ) (54)
We introduce the following auxiliary variables

ω D = 1 -η 2 θ = ω D t 1 C = cos θ S = sin θ A = µ s -µ d -2 η v -b B = v b +η A ω D D = A * ω D + η B X = A C + B S Y = v b C -D S α = -η t 1 E = exp α (55)
and recast the main equations as

r 1 (t 1 , t 2 , η, vb, ω D , . . . , E) = E X + 2η v b + v b t 2 + µ d -µ s r 2 (t 1 , t 2 , η, vb, ω D , . . . , E) = E Y -v b (56)
Continuation with respect to η or v b is performed by adding a third equation to stick one of these two parameters to a constant value. These diagrams are not shown here.

Application to finite element analysis in mechanics

As mentioned in the introduction, the Asymptotic Numerical Method has first been design for nonlinear finite element analysis in solid and structural mechanics. Several specialized finite element programs have been developed by the researcher working on the ANM each time a new formulation or a new class of nonlinearity was under study.

The interest of the generic Taylor series based continuation method presented here is to provide a unified framework for implementing a large class a finite element analysis whatever the complexity of the formulation or of the constitutive law. Firstly, recalling that the input of this generic method is only the quadratic recast of the equation ( 15) (( 13) and [START_REF] Cochelin | Méthode asymptotique numérique[END_REF] for the explanatory example), implementing a finite element analysis only requires a quadratic writing of the residual equation vector with suitable auxiliary variables. The consistent tangent stiffness matrix which is always a key issue in a finite element implementation is here automatically constructed from the sparse tensor lists and the block solving procedure. Secondly, the internal representation of the system with the sparse tensor lists allows to keep a very good performance of the computations even for complex mechanical model.

As an illustration of the performance of this method for finite element nonlinear analysis, the problem of geometrically nonlinear 2D elasticity in the framework of large displacements, small strains and a total Lagrangian formulation is addressed. A linear elastic constitute law is assumed for the simplicity of the presentation.

Consider an elastic body Ω submitted to a load growing proportionally to a load parameter λ. Let u be the displacement field, X = ∇u the displacement gradient tensor, F = I d + X the transformation gradient tensor, E = 1 2 (X + X t + X t .X) the Green-Lagrange strain tensor, D the elastic constant tensor, S and P the second and first Piola-Kirchhoff stress tensor.

Taking the displacement field u as the main variable and the tensor field X, S and P as the auxiliary variables, a possible formulation of this continuous problem is : Find u and X, S, P satisfying Auxiliary equations

   X = ∇u S = D : ( 1 2 (X + X t + X t .X)) P = (I d + X).S (57) main equation Ω δF : P dv -λ < F ext , δu >= 0 ∀δu (58)
The auxiliary equations are quadratic and follows the linear declaration rule. They allows to compute P from a given u. The main equation corresponds to the virtual work theorem and expresses the equilibrium inside Ω. Notice that this equation is linear with respect to P .

We now introduce a classical displacement based Finite Element method (see for instance [START_REF] Bathe | Finite element procedures[END_REF], [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF]) with classical isoparametric elements. After discretization, the main variable u becomes the vector of active degree of freedom u that collects displacements at the free nodes, and the auxiliary variable becomes the value of the tensor X, S and P at the Gauss point of the elements. Let introduce the notation u e for the vector of d.o.f of an element, P = P 11 P 12 P 21 P 22

T for the vector of components of the tensor P , (same notation for X, and S = [S 11 , S 22 , S 12 ] T ), G for the gradient-displacement matrix on an element.

The discrete form of the virtual work equation becomes the residual equilibrium at the free node

∪ Ωe G T P dv -λ F ext (59)
where ∪ stands for the classical assembly of the elementary internal forces vector. A classical Gauss integration is performed to evaluate these integrals at element level. Finally, the discrete problem is : Find U := u , λ and the auxiliary variables U aux = [ X , S , P , . . . ] at each Gauss point verifying Auxiliary equations at Gauss point

                       X = G u e S = D   X 11 + (X 2 11 + X 2 21 )/2 X 22 + (X 2 12 + X 2 22 )/2 (X 12 + X 21 + X 11 X 12 + X 21 X 22 )/2   P =     (1 + X 11 ) * S 11 + X 12 * S 12 (1 + X 11 ) * S 12 + X 12 * S 22 X 21 * S 11 + (1 + X 22 ) * S 12 X 21 * S 12 + (1 + X 22 ) * S 22     (60) main equation ∪ Gauss G T P det(J e )w e -λ F ext (61)
From these quadratic equations, the sparse tensor list are generated thanks to the use of the polarization formulas (42). Since the main variable is the vector of active d.o.f of the model, the jacobian matrix J = K -BK -1 aux C which is introduced for the block solving process is exactly the consistent tangent stiffness matrix of the finite element model.

In figures 5, the buckling behaviour of a sandwich-like beam with X stiffeners is considered as an illustrative example. The beam is clamped on the left and submitted to a uniform compressive force on the right. The bifurcation diagram is composed of a symmetric fundamental branch from which secondary branches bifurcate. The fundamental branch has been computed with 23 steps. Six bifurcation points have been automatically detected and located with the method explained in [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations[END_REF]. The first three bifurcated branches have been computed with 20 steps of continuation each. The deformed structure consists of a global buckling mode with local wrinkling of skins in the compressed parts. The three first deformed modes are shown in figure 5.

For this example, the mesh is composed of 282 quadrilateral 8 nodes elements with a 2 × 2 reduced integration scheme, the number of active d.o.f is 2616 and the number of auxiliary variables is 12408. The number of nonzero coefficients of the quadratic sparse tensor Q is here 31584. The computing time for these 83 steps of continuation is 25 seconds on a laptop computer (two processors 2,5 GHz Intel Core i7 with 16 Go RAM) with the Matlab software Manlab-4.0.

As discussed so far, the formulation of a mechanical model is not unique and it is also the case for its quadratic recast. A natural question is then "what is the best recast" and also "is it worth to have a minimal number of auxiliary variables". As an answer to these questions, the problem considered above has also been written by taking only the stress tensor S for the auxiliary variables. The consequence of not storing the gradient X and the first Piola-Kirchhoff P is a tremendous increase of the size of quadratic sparse tensor Q which pass from 32 terms per Gauss point for the previous formulation to more than a thousand. This results in a severe increase of the computing time of the right hand side F nl p of the linear systems to be solved [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF]. The corresponding gain of the computing time in the block solving procedure that results from a smaller system is not significant. Hence, it follows that the best formulation is not the one with a minimal number of auxiliary variables but the one with the minimal size for the quadratic sparse tensor here.

To finish these comments on the use of the proposed method for finite element analysis, we come back to the elementary transcendental function. For the sake of simplicity, the mechanical model used here for the illustration does not contain any elementary transcendental function. This is however not the case for most non linear finite element analysis. for instance, the functions ln, exp and power are frequently used for the modelization of the constitutive law of nonlinear material. For the structural models such as beams, plates, and shells, rotational degrees of freedom are frequently introduced to described the rotation of the section, and as a consequence the functions sines and cosines also frequently appear in the formulation. From the preceding presentation on how to manage these functions by using the adequate auxiliary variables and adequate companion functions (see the table in appendix A), it is clear that these more complex mechanical models can also be treated with this generic continuation method. It should be noticed that, whatever the complexity of the model and whatever the use of transcendental functions, the consistent tangent stiffness matrix is automatically constructed as explained above.

Conclusion

An efficient implementation of a Taylor series based continuation algorithm is presented. Its need of a quadratic formulation of the equations is explained and extended to the case of elementary transcendental functions. The decades of experience using the ANM continuation resulted in fourth version of the continuation software Manlab that is used here to treat successfully different systems. The accuracy of the computations is shown following the branches of the classical Logistic map and the twelve first flip bifurcations. The robustness of the method is demonstrated when applied to the Layne-Watson problem. The algorithm is then applied on a stick-slip problem, showing the possibility to make the continuation of a smooth by part system with this method. The last section presents a generic application to mechanical problems issued from FEM.

The quadratic recast of the system of equations is the key point of this continuation tool that allows to compute the Jacobian matrix and the Taylor series of the ANM automatically. The separation of the variables into main or auxiliary variables, the condensation of the auxiliary variables using a block solving technique and the sparse tensorial formalism ensure a very significant improvement of the computation time as compared to the previous versions of generic implementation of the ANM in the Manlab suite2 .

Figure 1 :

 1 Figure 1: Bifurcation diagram of the logistic map obtained with Manlab 4.0. The red dots are period doubling bifurcations. The x-axis is log(1 -µ µ crit ) where µ crit is the critical value above which the system becomes chaotic. The 12 first bifurcations were computed.

Figure 2 :

 2 Figure 2: Convergence of the ratio of the length of two consecutive intervals between two period doubling bifurcations to Feigenbaum constant δ. δ n = log 10 | µ n+2 -µ n+1µ n+1 -µn -δ| is plotted to show the exponential rate of convergence unambiguously.

Figure 3 :

 3 Figure 3: Bifurcation diagram of the Layne-Watson homotopy system (48). On the top with N = 10 and on the bottom with N = 50. On the plots, x 1 is in red and x N is in blue.

Figure 4 :

 4 Figure 4: Scheme of the system and a phase portrait for η = 0.1171 and v b = 0.1 taken from the full bifurcation diagram (not shown here).

Figure 5 :

 5 Figure 5: Bifurcation diagram of a sandwich-like beam with X stiffeners showing the three first deformed modes. Only an excerpt of the bifurcation diagram is shown here. The first, the second and the third mode are shown here and corresponds respectively to the red dot A, B and C.

The term "numerical" stands in the name of the method because of this finite element discretization.

Manlab 4.0 is available online on the dedicated website https://manlab.lma.cnrs-mrs.fr/ .

Acknowledgment

The authors want to thank Stéphane Bourgeois for fruitful discussions on the finite element analysis part. This work has been carried out in the framework of the Labex MEC (ANR-10-LABX-0092) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the Investissements d'Avenir French Government program managed by the French National Research Agency (ANR).

A Table of quadratic recasts of the most common elementary transcendental functions

auxiliary variable companion variable(s) full quadratic recast differentiated form


The aim of the present work is to determinate an accurate formulation for the seaspray source function

Introduction

Material and Methods

Results and Discussion

The simulations presented in Fig. 3 were made for a local wind speed of 9,8 m s -1 and for a particle radius of 5 µm using various S3F and compared to our experimental data. For each fetch, 20 experimental data at a minimum were used for our calculations:

The experimental results reported in Fig. 3 show a peak in the evolution of the sea-spray source function with respect to the fetch. At very short fetches, the seaspray generation processes tend to increase with fetch whereas at large fetches, they slightly decrease with increasing fetch. This can be compared to the evolution of the whitecap fraction with the wave age as already noted by Lafon et al. (2004;2007). The wave age is defined as the ratio between the wave phase velocity Cp and alternatively the wind speed U or the wind friction velocity 𝑈 * . We have therefore plotted in Fig. 4 the variation of the sea spray production versus the wave age which confirm the results previously plotted in Fig. 3.

Conclusion References

Field site and instrumentation

The data set used in this present work is issued from the measurements of aerosol particle size distributions which took place on the island of Porquerolles (Fig. 1) and on board of the ship Atalante during the Miramer campaign (Fig. 2).

The aerosol data were acquired in the 0.1-45 µm size ranges using two particle measuring systems (PMS): one active scattering spectrometer probe (ASASP) and one classical scattering spectrometer probe (CSASP).

The MACMod model

The general governing equation for the number concentration of aerosols N at a given radius r, can be written as the budget equation:

𝜕𝑁 𝜕𝑡 + 𝑈. 𝛻 𝑁 = -𝛻. 𝜙 + 𝑆
where 𝑼 is the air velocity, 𝝓 is the aerosol flux and 𝑺 is the usual volume source-sink term accounting for condensation, evaporation, coagulation and nucleation.

According to Fairall and Davidson (1986), the vertical flux 𝝓 can be expressed as:

𝜙 = 𝑤′𝑛′ -𝑉 𝑔 𝑋 -𝐷 𝑚 𝜕𝑁 𝜕𝑧
where 𝒘′ and 𝒏′ denote turbulent fluctuations, 𝑽 𝒈 is the gravitational settling velocity, 𝑫 𝒎 is the molecular diffusion coefficient and 𝒛 is the height above the surface. With 𝐴 1 = 4,5 ; 𝐴 2 = 0,408 ; 𝐴 3 = 0,931

Sea-spray source function methods

And 𝐵 1 = -1,53-log (𝑟 80 ) 0,55

; 𝐵 2 = -0,51-log (𝑟 80 ) 0,57

; 𝐵 3 = -0,57-log (𝑟 80 ) 0,52 is the size dependent SSA production flux, 𝒊 is the mode number, 𝑭 𝒊 (𝑹𝒆 𝑯 𝑾 ) is the flux for mode 𝒊 , 𝑫 is the dry particle diameter, 𝜎 is the geometrical standard deviation and 𝑪𝑴𝑫 is the count median diameter.

Mediterranean

  Institute of Oceanography (MIO-UM 110) -University of Toulon S.Laussac, J.Piazzola, G.Tedeschi Development of a sea salt emission parameterization for the North-western Mediterranean A sea-state dependent model for the sea-spray generation function has been implemented in the MACMod model and compared to experimental data acquired in the Mediterranean area. The simulations provide a relevant variation of the aerosol concentrations with the fetch evolution.

• The Ovadnevaite et al. (2014) flux formulation

  The S3F formulation byOvadnevaite et al. (2014) is based on the revisited dimensionless breaking wave parameter proposed by Zhao and Toba (2001) 𝑹𝒆 𝑯 𝑾 . The S3F formulation is expressed as a combination of five lognormally distributed modes for different droplet sizes:

		𝑑𝐹 𝑑𝑙𝑜𝑔𝐷	=	5 𝑖=1	𝐹 𝑖 (𝑅𝑒 𝐻 𝑊 ) 2𝜋 𝑙𝑛𝜎 𝑖	exp -	1 2	𝑙𝑛	𝐷 𝐶𝑀𝐷 𝑖 𝑙𝑛𝜎 𝑖	2
	with	𝒅𝑭 𝒅𝒍𝒐𝒈𝑫