M H Malik

D Borzacchiello
email: domenico.borzacchiello@ec-nantes.fr

F Chinesta

P Diez

Inclusion of Frequency Dependent Parameters in Power Transmission Lines Simulation using Harmonic Analysis and Proper Generalized Decomposition

Keywords: Skin Effect, Frequency Dependent Parameters, PGD, Transmission Lines, Fractional Derivative

This study presents the application of Proper Generalized Decomposition on transmission line models involving frequency dependent parameters. Frequency dependence of parameters can be due to a number of reasons including phenomena such as "Ground Return Effects", "Proximity Effects" and "Skin Effects". In the current study, simplified methods of skin effects based on Bessel functions are used in order to showcase the method. Although we implement our study using a specific model for skin effects to demonstrate the effectiveness of the proposed method to accomodate frequency dependent effects in PGD, the present work is not meant to discuss the merits of any particular skin effects model. The method can easily accommodate other effects which induces frequency dependence in the transmission line parameters. In time domain modeling, the parameters are assumed constant

Introduction

Accurate modeling of transmission lines is important in the simulation and subsequent transient analysis of power grids. In this framework, frequency dependence of the line parameters should be accounted for [START_REF] Sérgio Kurokawa | Behavior of Overhead Transmission Line Parameters on the Presence of Ground Wires[END_REF]. The parameters of transmission lines in practice are highly frequency dependent and it is necessary to accommodate the frequency dependence over the whole frequency range of interest.

AC transmission lines are affected by the frequency of the AC voltage which produces a secondary effect on the line resistance and inductance. This effect is related to the non uniform distribution of the current which tends to concentrate more towards the outer surface of the line away from the center as the frequency increases [START_REF] Admane | Simulation of electric transmission lines with skin effect based on a reduced order differential model[END_REF]. Skin effect results in the reduction of the effective cross-section of the line used for the transmission of current and hence the resistance increases. This concentration of current at the surface has an inverse relation with the internal inductance of the line, as the frequency increases the inductance decreases [START_REF] Coelho | Proposal of an alternative transmission line model for symmetrical and asymmetrical configurations[END_REF][START_REF] Bormann | Reluctance Network Treatment 28 of Skin and Proximity Effects in Multi-Conductor Transmission Lines[END_REF]. There are other phenomena like ground return and proximity which also have an effect on the transmission line parameters resulting in the frequency dependence of these parameters [START_REF] Coelho | Proposal of an alternative transmission line model for symmetrical and asymmetrical configurations[END_REF].

If the phenomena which produce frequency dependence are ignored and the parameters are assumed to be constant, this results in a magnification of the higher harmonics [START_REF] Marti | Accurate Modeling of Frequency-Dependent Transmission Lines in Electromagnetic Transient Simulations[END_REF]. The wave shapes are distorted and at higher frequency the peaks have higher magnitudes than that are observed in practice. There are a wide range of frequencies present in the signals during the transient phase and this can be modeled with the frequency-dependent parameters.

The challenge in time domain modeling is to incorporate the frequency dependence in an efficient manner and, as such, several formulations have been proposed. As it will be discussed later in this paper, many of these result in profound modifications of the original model that are sometimes difficult to interpret in the light of their physical meaning. Accounting for frequency related effects is obviously more natural when the governing equations are formulated in the frequency domain rather than in time. However, this implies that incremental time stepping strategies for transient simulation can no longer be applied and harmonic components have to be solved individually, each requiring the solution of a different problem.

In this work we propose a robust parametric solver that allows to determine at once all the harmonics in a given frequency range. While this step is performed "off-line", the actual time response can be efficiently computed "on-line". The rest of this paper is organized as follows, Section 2 presents the literature review. Section 3 presents the method proposed in the current study to deal with the frequency dependent parameters in the frequency domain. In the Section 4 we present the skin effects model for the transmission lines and in Section 5 we present the distributed parameters model of transmission lines that we adopted in this study. Finally, we present the results of our proposed approach in the Section 8.

Transmission Line Modeling

The frequency dependence of the parameters is the consequence of the assumption made in the derivation of Telegrapher's equations from the Maxwell's equations [START_REF] Nucci | Derivation of telegrapher's equations and field-to-transmission line interaction[END_REF]. It is assumed that the conductor carries a uniform current in a uniform time-harmonic magnetic field. Furthermore, it is assumed that the curvature of the wire is negligible and there is no field variation in the wire direction i.e. x-axis in this case. It is considered that the current flows in x-direction only [START_REF] Igarashi | Semi-Analytical Approach for Finite Element Analysis of Multi-turn Coil Considering Skin and Proximity Effects[END_REF]. These assumptions simplify the 3-D Maxwell's equation to a manageable 1-D equation known as Telegrapher's equation given by:

∂v(x, t) ∂x = -Ri(x, t) -L ∂i(x, t) ∂t (1) ∂i(x, t) ∂x = -Gv(x, t) -C ∂v(x, t) ∂t (2)
x

∈ [0, l] ; t ∈ [0, ∞]
where, v, i are the voltage and current, R, L, G and C are the resistance, inductance, conductance and the capacitance per unit length respectively.

Eliminating i and combining the eqs. 1 and 2 yields

LC ∂ 2 v(x, t) ∂t 2 + (LG + RC) ∂v(x, t) ∂t + RGv(x, t) = ∂ 2 v(x, t) ∂x 2 (3)
In order to retain the correct physical behavior, the model parameters like resistance and inductance are often considered as function of the frequency in order to compensate for the simplifications introduced by the assumptions in the derivation of the 1-D model. The other scenario in which these parameters are kept constant will then require the 3-D partial differential equations to be solved.

Several numerical methods have been proposed for the solution of frequency dependent models [START_REF] V Maló Machado | Magnetic field analysis of three-conductor bundles in flat and triangular configurations with the inclusion of proximity and skin effects[END_REF]. An overview is listed in the Table 1.

There are mainly two families of solution procedures, i.e., a time-domain simulation and a frequency-domain simulation. Both families of methods have their pros and cons and have utility suitable to certain conditions. Timedomain modeling is delicate and some numerical instabilities and accuracy issues have been reported in the literature [START_REF] Marti | Accurate Modeling of Frequency-Dependent Transmission Lines in Electromagnetic Transient Simulations[END_REF][START_REF] Suk | Accurate transient simulation of transmission lines with the skin effect[END_REF]. But, the advantage it provides is the straightforward compatibility with models for components of the power system.

Most of the studies, reported in Table 1,avoid convolutions and developed intelligent models to circumvent the difficulties in time-domain modeling by modifying the formulation of the equations in time. For example, in the study by Semlyen and Dabuleanu [START_REF] Semlyen | Fast and accurate switching transient calculations on transmission lines with ground return using recursive convolutions[END_REF], the authors have used approximations based on the exponential functions. The study by Chu-Sun et al [START_REF] Yen | Time-domain skineffect model for transient analysis of lossy transmission lines[END_REF] uses equivalent circuit with each section is a combination of equivalent resistances and inductances to include the skin effects. Some recent studies like Marques da Costa et al [START_REF] Coelho | Proposal of an alternative transmission line model for symmetrical and asymmetrical configurations[END_REF] used fitting methods to include skin effects in lumped parameters model of transmission line.

On the other hand, research studies in favor of frequency-domain modeling point out the ease of formulation. Instead of time convolutions, we can readily multiply the harmonic response of the system with the desired input function. Then obtaining the solution in time-domain is simply just a mat-ter of inverse transformation. Wilcox and Condon [START_REF] Wilcox | A new transmission-line model for timedomain implementation[END_REF] proposed an inverse transformation method based on ARMA (Auto-Regressive-Moving-Average) fitting method to obtain time-domain model, where ARMA is a statistical model used in several fields to fit model to observed data. While other studies like by Kurokawa et al [START_REF] Sérgio Kurokawa | Behavior of Overhead Transmission Line Parameters on the Presence of Ground Wires[END_REF] used modal transformation matrix. In frequencydomain analysis, as described in the study of Meyer and Dommel [START_REF] Meyer | Numerical modelling of frequencydependent transmission-line parameters in an electromagnetic transients program[END_REF], we have to use linear superposition of all the frequency components. The solution is performed for every individual frequency present in the transients, as the unit impulse response in time, "translates" as the unit heaviside response in frequency-domain.

In the current study, we present a simple yet efficient frequency-domain method based on the Proper Generalized Decomposition for the solution of transmission line including the frequency-dependence of parameters due to skin effects. We propose a method based on PGD's separated representation as a monolithic solver for all the frequencies.

Proposed Approach

Refering to eq. ( 3), the second term is referred as the dissipation term and the last term is the dispersion term. For a lossless transmission line we assume G = 0. By introducing a suitable finite element basis, the function v(x, t) can be approximated as

v(x, t) ≈ n i=1 v i (t)φ i (x). (4)
Then, imposing Galerkin orthogonality of the residual of equation 3 on the basis functions φ i (x) leads to the following system of ordinary differential equations in time

[M ]{v(t)} + [C]{ v(t)} + [K]{v(t)} = {f (t)} (5)
where Time integration of Eq. 5 requires the solution of a n×n algebraic system for every time increment for nonlinear or time dependent matrices. Modal analysis provides an efficient strategy to reduce the numerical complexity of the system by transforming it into a set of n uncoupled equations.

The generalized eigenvalue problem to be solved for the system of eq. ( 5),

with damping [C] removed,

λ 2 [M ] + [K] {V i } = {0} (6)
where, the vectors {V i } are the eigenvectors and can be grouped in a ma-

trix [Λ] = [{V 1 }, {V 2 }, • • • , {V n }].
The matrix [Λ] can be used to diagonalize the matrices [M ] and [K] as:

[Λ] T [M ][Λ] = [m] [Λ] T [K][Λ] = [k] (7)
The system presented in Eq. ( 5) involves damping, there are generally three ways to proceed further in cases of non-zero damping. The first method is to use the direct time integration which we already explained could require solution of n×n algebraic equations. The second method is to setup and solve a different complex eigensystem in which eigenvalue problem diagonalizes two matrices that contains [M ], [C] and [K] as sub-matrices [START_REF] Kurdila | Fundamentals of Structural Dynamics[END_REF]. The final and the easiest method is to continue with the diagonalization presented in eq. ( 7).

For that purpose, we need the damping matrix [C] in a linear combination

of the matrices [M ] and [K] such that [C] = a[M ] + b[K]
. This is usually the case for the transmission line models because the parameters are constant along the transmission line. Therefore it can also be diagonalized using the

same matrix [Λ]. [Λ] T [C][Λ] = a[m] + b[k] (8)
Also, the variable of interest is transformed into the modal basis as:

{U } = [Λ] T {V } (9)
This diagonalization gives us an uncoupled ordinary differential equations that can be solved separately. Writing the equations in indicial notation,

m i üi + (am i + bk i ) ui + k i u i = 0 (10)
The solution of eq. ( 9) is straightforward and once we have the solution in terms of modal basis, we can transform the solution back to the physical basis.

{V } = [Λ]{U } (11)
The case we address here is when the parameters resistance R(ω) and inductance L(ω) becomes frequency dependent because of skin or other effects. The problem where we have frequency-dependent parameters can not be solved needs special attention. Studies such as the one performed by Wilcox and Condon used the modal analysis but modeled the system with frequency-dependent parameters first in frequency-domain and careful manipulation of the system itself is required, otherwise the transformation from frequency-domain model to time-domain model may cause the loss of causality as seen by Crandall [START_REF] Crandall | The role of damping in vibration theory[END_REF] in the context of structural dynamics. In his work, Crandall discusses the effect of frequency-dependent parameters on time integration and introduced the concept of non-equations. The better option is to solve the harmonic analysis problem for every frequency ω. Transforming the problem from time-domain to frequency-domain using Fourier transformation of the unknown variable v i (t) and forcing function f (t) that write

V (ω) = ∞ -∞ v(t)e -ιωt dt F (ω) = ∞ -∞ f (t)e -ιωt dt (12)
Therefore, the equation ( 5) after transformation into frequency-domain is given as:

-ω 2 [M (ω)] + ιω[C(ω)] + [K(ω)] V (ω) = F ∀ ω (13)
Now, the harmonic analysis problem requires solution of n systems of linear equations where n represent the number of discrete frequencies considered. A similar problem involving fractional damping was discussed in [START_REF] Germoso | From standard to fractional structural visco-elastodynamics: Application to seismic site response[END_REF].

The use of proper generalized decomposition (PGD) as a paramteric solver as presented in earlier studies by Chinesta et al. [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF] and [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF] is considered here as well. PGD provides a strong tool as a single monolithic solver as compared to solving the harmonic problem using eq. 13 individually for every frequency. PGD has already been implemented for the parametric power flow problem [START_REF] García-Blanco | A reduced order modeling approach for optimal allocation of distributed generation in power distribution systems[END_REF] and an error estimation strategy has also been developed [START_REF] García-Blanco | Monitoring a pgd solver for parametric power flow problems with goal-oriented error assessment[END_REF].Note that if we solve the linear problem of eq. 13 we have to compute the matrices for every single frequency as they are no longer constant compared to the case where the constant parameters was considered and depend upon the frequency. Therefore, in this study, we present the Proper Generalized Decomposition (PGD) as a parametric solver, a method capable to solve the harmonic analysis.

Skin Effects Model

Skin Effect results in frequency dependence of the transmission line parameters. This phenomenon is the result of unequal distribution of current across the cross-section of the transmission line, this is more pronounced in long transmission lines. The parameters that become functions of frequency are resistance and inductance. Inductance of a conductor is related with the flux linkage, the flux linkage in a conductor happens internally and externally. Internal flux linkage is due to the self current of the conductor while the external flux linkage is due to the presence of adjacent conductors.

Therefore, total inductance of a transmission line is the sum of inductances due to both internal and external flux linkages. The evaluation of external inductance is based on the assumption that the total current flow is on the surface of the transmission line, i.e., external inductance is maximum when the skin effect is maximum. However, the external inductance is the result of the current flow in the neighboring transmission line and therefore remains constant. In the current case, we are dealing with the single transmission line assuming no adjacent transmission lines and in the succeeding discussion wherever inductance is mentioned it is assumed internal inductance unless explicitly mentioned.

Typically transmission lines are bundled aluminum strands twisted together with one kind of transmission line consists of steel strands in the core to increase the strength of the transmission line. However, for the sake of simplicity the discussion is based on a solid cylindrical conductor to effectively demonstrate the model.

Several researchers have presented simple models to include skin effects in the transmission line models [START_REF] Suk | Accurate transient simulation of transmission lines with the skin effect[END_REF], [START_REF] Gatous | Frequency-dependent skin-effect formulation for resistance and internal inductance of a solid cylindrical conductor[END_REF] and [START_REF] Admane | Simulation of electric transmission lines with skin effect based on a reduced order differential model[END_REF]. Since, our focus in the current study is not to model the skin effects but rather use the existing model and present a solution that is able to perform transient simulations with skin effects included. Therefore, we present here some of the simplified models and use them in our study.

For a transmission line with large radius compared to the skin depth [START_REF] Lucht | Transmission Lines and Maxwell's Equations[END_REF],

δ << 4r where skin depth δ is calculated as [START_REF] Kurdila | Fundamentals of Structural Dynamics[END_REF].

δ = 2 ωµσ 1/2 (14)
where, r is the radius of the cross-section of the conductor, ω is the angular frequency and related to frequency in Hz as ω = 2πf , µ is the magnetic permeability and σ is the electrical conductivity of the conductor.

The resistance R(ω) and inductance L(ω) are given in the study of Lucht [START_REF] Lucht | Transmission Lines and Maxwell's Equations[END_REF] based on the expansion of Bessel functions with the assumptions listed here. These simplified functions for resistance and inductance for high frequencies is given as:

R(ω) = 1 σ(2πr)δ , L(ω) = 1 σ(2πr)δω (15)
In the technical report by Phil Lucht [START_REF] Lucht | Transmission Lines and Maxwell's Equations[END_REF], a low frequency limit is defined, above this frequency resistance and inductance become the functions of frequency, below this limit the resistance and inductance are given by

R = 1 σπr 2 , L = µ 8π (16)
The low frequency resistance is the same as the DC resistance of a conductor in per unit length given by

R = ρ πr 2 ( 17
)
where ρ is the resistivity of the conducting material.

The article by Monteiro et al [START_REF] Humberto | Simplified skin-effect formulation for power transmission lines[END_REF] presents a simplified formulation for the calculation of skin effects on the transmission line wires. The formulation is based upon the Fourier transformation of the Maxwell's equation [START_REF] Chipman | Theory and Problems of Transmission Lines[END_REF].The methodology presented by Monteiro et al is simple in terms of its application without the need to use Bessel functions and also presented good accuracy as it is directly developed from the Maxwell's wave equations. where

R(ω) = ∞ k=1 (R k /(R 2 k + ω 2 L 2 )) ( ∞ k=1 R k /(R 2 k + ω 2 L 2 )) 2 + ω 2 ( ∞ k=1 L/(R 2 k + ω 2 L 2 )) 2 (18)
L(ω) = ∞ k=1 (L/(R 2 k + ω 2 L 2 )) ( ∞ k=1 R k /(R 2 k + ω 2 L 2 )) 2 + ω 2 ( ∞ k=1 L/(R 2 k + ω 2 L 2 )) 2 (19)
R k = ξ 2 k 4πσr 2 ; L = µ 4π and ξ k = (2k -1)π 2 + π 4 (20)
The units of R k is Ω.m and L is H/m. The units of electrical conductivity σ is S/m and of the magnetic permeability µ is H/m. The term ξ k is the k th Bessel function root and is directly proportional to k. As mentioned previously, the transmission line in this study is considered to be a solid cylindrical copper wire, its radius and the electrical properties of the copper are listed in Table 2.

A comparison of different models was performed and the results were compared with the results from the study of Monteiro et al [START_REF] Humberto | Simplified skin-effect formulation for power transmission lines[END_REF]. The results for the resistance and inductance with respect to the frequency is presented in the Figures 1a and1b.

Distributed Parameters Model

Time Domain Simulation

We start with the time domain distributed parameters (DP) model. The electrical transmission lines provide the path for the voltage and current waves to propagate, the dynamics of the propagation is described by a distributed parameters model. This model depends upon the considerations taking into account like the losses and the frequency dependence of the parameters.

The distributed parameters model for a lossy transmission line with constant parameters can be written in the time domain as:

∂v(x, t) ∂x = -L ∂i(x, t) ∂t -Ri(x, t), ∂i(x, t) ∂x = -C ∂v(x, t) ∂t (21)
for 0 ≤ x ≤ l, where v(x, t) and i(x, t) are the voltage and current in the line respectively, and the resistance, the inductance and the capacitance are represented by R, L, C respectively and these constants are positive and independent of frequency. The parameters R, L, C are given in per unit length and in distributed parameters model are used without multiplication with the length of the line. However, these parameters must be multiplied by the length of the line in case one is using any lumped circuit model for transmission lines e.g. pi-equivalent model. The effect of transverse conductance G of the line is neglected in this model. The interconnection between the transmission line, the load and the generator is performed using the Kirschhoff's laws and the boundary values.

v(x, t)| x=0 = V 0 (t) i(x, t)| x=l = 0 (22)
The term V 0 is the voltage source assumed to represent a generator, the boundary condition given above assumes that the voltage source is ideal, that is, the impedance is zero (Z 0 = 0) and the impedance at the receiving end

(x = l) is infinite (Z l = ∞).
The equation is converted to a single variable second order differential equation in time and space.

∂ 2 v ∂x 2 -LC ∂ 2 v ∂t 2 -RC ∂v ∂t = 0 (23)
The equation ( 23) is then discretized in time and space, using finite difference scheme for the time integration and finite elements for the space.

For a copper conductor wire, we can calculate the resistance and inductance from the equations ( 16), these values are listed in the table 3. Solving a space time problem for these values, with a boundary condition representing a step as an input voltage V 0 .

V 0 =    0 for t 0 1 for t > 0 (24)
Also, a sine function with a frequency of f = 100Hz was used as an input voltage for the same system

V 0 =    0 for t 0 sin(2π100t) for t > 0 ( 25
)
∂ V (x, ω) ∂x = -Î(x, ω).(R + ιLω), ∂ Î(x, ω) ∂x = -V (x, ω).(ιCω) (26)
Combining the two equations given in [START_REF] Davila | Practical implementation of a transmission line model for transient analysis considering corona and skin effects[END_REF], we get a second order ordinary differential equations for voltage.

∂ 2 V (x, ω) ∂x 2 + (LCω 2 -ιRCω) V (x, ω) = 0 ( 27
)
Corresponding boundary conditions in the frequency domain is given as

V (x, ω) x=0 = V 0 ∂ V (x, ω) ∂x x=l = 0 (28)

Proper Generalized Decomposition for Constant Parameters Model

The problem presented in Sec. 5 can be solved using PGD that assures the separated representation,

V (x, ω) = m i=1 X i (x).O i (ω) ( 29
)
For the details about PGD formulation and its advantage, one can refer to several publications for example, [START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF][START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF]. In general, PGD methods comprises of an offline and an online part. In the offline stage of the solution, the problem is formulated in the separated representation and the solution database is created. The online part then is just to compute the solution for the particular problem by multiplying the separated coordinates for each enrichment mode.

Once, a solution using PGD with the boundary conditions given in [START_REF] Kim | Compact equivalent circuit model for the skin effect[END_REF] is obtained, it is straightforward to get the desired time response using the discrete inverse Fourier function available in MATLAB. Using the separated representation of (29), the DP model can be expressed in separated repre-sentation. In order to develop a PGD formulation, first we must have a weak formulation of the problem.

Multiply equation ( 27) with a test function δV * (x) where, δV represents the test function and δV * (x) represents the complex conjugate.

Ωx×Ωω δV * ∂ 2 V ∂x 2 dx dω - Ωx×Ωω δV * f (ω) V dx dω = 0 (30)
where, f (ω) = (LCω 2 -ιRCω).

Here, we briefly describe the PGD algorithm, for details the interested readers can refer to the literature [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF][START_REF] Chinesta | PGD-Based Computational Vademecum for Efficient Design, Optimization and Control[END_REF]. PGD follows the greedy algorithm approach, i.e., one mode is evaluated at a time. The mode is assumed to be converged after sufficient iterations result in no further improvement. This procedure of adding modes is termed as the enrichment process. Assuming that n -1 initial modes are already converged, we define the n th mode to be evaluated as:

V (x, ω) = n-1 i=1 X i (x).O i (ω) + X p n (x)O p-1 n (ω) (31)
We assume the test function δV * (x, ω) in the separated form as

δV * (x, ω) = O * (ω)δX * (x) + X * (x)δO * (ω) (32)
In eq. ( 31), the superscript p indicates the number of fixed point iterations. For individual modes evaluation, we apply the fixed point loop which is to find solution for each coordinate space alternatively.

Since, the problem is now separated in the space X(x) and the frequency O(ω), we have to solve two problems using FEM once for the space domain and then for the frequency domain. The weak formulation for the space only problem is:

Ωx×Ωω dδX * n dx dX p n dx O p-1 * n O p-1 n dx dω + Ωx×Ωω f (ω)δX * n X p n O p-1 * n O p-1 n dx dω = - Ωx×Ωω n-1 i=1 dδX * n dx dX i dx O p-1 * n O i dx dω - Ωx×Ωω f (ω) n-1 i=1 δX * n X i O p-1 * n O i dx dω + Ωω δX * n O p-1 * n Î(ιLω + R) dω (33)
and for the problem to be solved in the frequency domain, the weak formulation will be

Ωx×Ωω dX p * n dx dX p n dx δO * n O p n dx dω + Ωx×Ωω f (ω)X p * n X p n δO * n O p n dx dω = - Ωx×Ωω n-1 i=1 dX p * n dx dX i dx δO * n O i dx dω - Ωx×Ωω f (ω) n-1 i=1 X p * n X i δO * n O i dx dω + Ωω X p * n | x=0 δO * n Î(ιLω + R) dω (34)
The tolerances selected for the stopping of the fixed point algorithm ε p and the tolerance for termination of enrichment process ε n , along with the number of modes n in both the constant parameters and frequency-dependent parameters are listed in Table 4.

The boundary condition in eq [START_REF] Kim | Compact equivalent circuit model for the skin effect[END_REF] is an equivalent of having an impulse function at the source which gives the impulse response or the transfer function Ĥ of the transmission line. The voltage impulse response at the receiving end is given as:

Ĥ(x = l, ω) = n i=1 X i (x = l)O i (ω) (35)
Once, this is done, the time response to any arbitrary condition at the source node is then straightforward multiplication of the voltage function in frequency domain with the particular solution of the PGD at the receiving node. Then the frequency response to an arbitrary function V0 is given as

V (x = l, ω) = Ĥ(x = l, ω) • V0 (36)
The input voltage V0 is usually a step function or a sinusoidal function in time. In practice, these functions are transformed into frequency-domain using the fast Fourier transform algorithm to compute the discrete Fourier transform. Similarly, the response V is inverse transformed to time domain using the same fast Fourier algorithm using discrete samples in the frequency interval.

Verification of PGD results

In this section, we first present the time response using inverse fourier transform of the PGD solution particularized at the receiving end.

The PGD solution in terms of voltage on a space-frequency domain is presented in Figures 3a to 3c.

The results from Figures 5a and5b shows good correlation between the solution from time integration and the PGD solution. Therefore, it is estab- order backward difference scheme. The blue curve of time simulation tends to approach the red curve of PGD solution as finer time discretization is used but even then spikes observed in the PGD solution are not captured by the time simulation. In comparison, PGD based approach uses frequency domain which by default is a higher order approximation and hence able to capture a wider range of dynamics.

PGD for frequency dependent parameters

The benefit of solving DP model in the frequency domain in the separated representation of PGD is the easy, fast and accurate solution with the option of a parametric solution. PGD methods provide the best option for parametric problems both in terms of computational efficiency and ease of solving the problem. In the case of skin effects, the equation ( 27) now has parameters resistance R and inductance L as a function of frequency ω. As explained earlier time integration of such problems is almost impossible to achieve and hence PGD provides real accurate option to solve such problems.

The problem is similar to the case of fractional derivatives which is similar to the fractional RLC circuit [START_REF] Francisco | Fractional RC and LC electrical circuits[END_REF].

Eq [START_REF] Gustavsen | Calculation of electromagnetic transients in transmission cables and lines taking frequency dependent effects accurately into account[END_REF] with R and L as functions of frequency reads as:

∂ 2 V (x, ω) ∂x 2 + (L(ω).C.ω 2 -ιR(ω).C.ω) V (x, ω) = 0 (37)
The increase in the resistance with the frequency and a simultaneous decrease in the inductance results in a system more damped than the system with constant resistance and inductance. This effect can be noticed in the Figures 6a to 6c showing the voltage in the space-frequency domain, there is only one peak visible with a much less amplitude.

The effect of increased damping due to the increase in resistance can be compared in the receiving end voltage for the frequency dependent and constant paramters solution given in the Figures 7. In As an initial step, we introduce length of the transmission line as a new dimension for the parametric problem. The solution now is given as:

V (x, ω, l) = n i=1 X i (x)O i (ω)D i (l) (38)
where D i represents the modes of the dimension introduced to represent length.

The mathematical formulation of the problem remains similar and there is no significant additional computational cost to the problem. For the case, we have reported in this study, the new dimension of transmission line length is considered for a range (200, 300) km. The discretization for the space and frequency have been kept the same as in previous case. The discretization and the other information valuable to gauge the PGD solution is listed in Table 5.

The advantage of parametric solution is presented in Figure 10 where the variation of voltages along the transmission line with different lengths for a frequency of 60 Hz is plotted. Another advantage of such a solution is in the sensitivity or uncertainty analysis, in such an analysis a fixed point is selected and a variation of usually 5% in the independent variable and its effect on the variable of interest is studied. With a PGD solution, we already have a complete solution as a curve and the derivative can be easily evaluated for the sensitivity analysis.

Conclusions

The work presented in this study, shows the effectiveness of frequencydomain methods for problems involving frequency-dependent parameters. The main contribution of the study presented in this paper is the development of a reduced order model able to include phenomenon like skin effects that cause the parameters to become functions of frequency. The method is equally applicable for the simulation of transmission line models including the ground return and proximity effects which has a similar effect on the parameters, i.e., the frequency dependence. In particular, harmonic analysis using PGD is fast, accurate and simple method for this kind of parametric problem. In this paper, we have demonstrated the method of PGD, constructing the database of the parametric solution offline and particularize the solution online for a specific set of parameters online in real time. PGD is an a-priori model order reduction method and is particularly effective to obtain solutions online from the pre-computed database of solutions. With the ability of including parameters such as the length of the transmission line or material properties of the conducting material as one of the coordinates of the separated representation, it truly provides the capability of having real time solutions using not so power platforms. In contrast to the method presented in the current study, there are solution methods available that include the frequency-dependence of parameters and simulate the model in time-domain, but these methods are complicated to implement and since they use time integration it is computationally costly and time consuming.

PGD is already established as a fast and accurate method for reduced order modeling. Another potential application of this method could be fractional RLC circuits.

Figure 1 :

 1 Figure 1: Frequency dependent resistance and inductance of a conductor

Figure 2 :

 2 Figure 2: Time simulation of constant parameters model in space-time domain

  .

6 .

 6 Frequency-Domain Model Skin effects phenomena results in frequency dependence of DP model, direct time integration of which is difficult to perform. But, the problem in the frequency domain becomes a linear second order differential equation, where the parameters resistance and inductance are now a function of frequency. Fourier transformation of the time domain distributed parameters model gives the equation in the frequency domain as:

Figure 3 :

 3 Figure 3: Voltage in space-frequency domain with constant parameters

Figure 4 :Figure 5 :

 45 Figure 4: Time response by inverse Fourier of the PGD solution

Figure 8 ,

 8 the two responses from the constant parameters and frequency dependent parameters are plotted on top of each other for a direct comparison. The increased damping removes almost all the oscillations observed in the constant parameters solution of Figures 4a and 4b. The difference between the constant parameters model and frequency dependent parameters solution are very different. Therefore, just to elaborate the effect of different values of the parameters we used the constant parameters model with parameters values evaluated at 60 Hz and 120 Hz and plotted in the Figure 9.

8. 1 .

 1 Parametric SolutionThe advantage of PGD based harmonic analysis is that any parameter can be introduced as a model extra-coordinate into the separated representation.

Figure 6 :

 6 Figure 6: Voltage in space-frequency domain with skin effects included

Figure 7 :

 7 Figure 7: Time response from PGD solutions for frequency dependent parameters

Figure 8 :

 8 Figure 8: Comparison of response between constant and frequency dependent parameters

Figure 9 :

 9 Figure 9: Comparison of response for constant parameters with parameter values at 60 Hz and 120 Hz

Figure 10 :

 10 Figure 10: Variation of voltage along the transmission line for different lengths at 60 Hz

Authors

Phenomena

Methods Time Domain Modeling Semlyn and Dabuleanu [START_REF] Semlyen | Fast and accurate switching transient calculations on transmission lines with ground return using recursive convolutions[END_REF] Ground return effect Approximation of time variable characteristic admittance using exponential functions J. R. Marti [START_REF] Marti | Accurate Modeling of Frequency-Dependent Transmission Lines in Electromagnetic Transient Simulations[END_REF] Ground return effect Weighting Functions, Characteristic Impedance

Chu-Sun et al [START_REF] Yen | Time-domain skineffect model for transient analysis of lossy transmission lines[END_REF] Skin effect Equivalent Circuit using parallel equivalent resistors and inductors Suk Oh [START_REF] Suk | Accurate transient simulation of transmission lines with the skin effect[END_REF] Skin


Rubén Ibáñez

Adrien Scheuer

Emmanuelle Abisset-Chavanne
email: emmanuelle.abisset-chavanne@ec-nantes.fr

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Antonio Huerta
email: antonio.huerta@upc.edu

Roland Keunings
email: roland.keunings@uclouvain.be

Pimm

A simple microstructural viscoelastic model for flowing foams

Keywords: Flowing foams, Viscoelasticity, Conformation, Microstructural description

The numerical modelling of forming processes involving the flow of foams requires taking into account the different problem scales. Thus, in industrial applications a macroscopic approach is suitable, whereas the macroscopic flow parameters depend on the cellular structure: cell size, shape, orientation, etc. Moreover, the shape and orientation of the cells are induced by the flow. A fully microscopic description remains useful to understand the foam behaviour and the topological changes induced by the cell elongation or distortion, however, from an industrial point of view, microscopic simulations remain challenging to address practical applications involving flows in complex 3D geometries. In this paper, we propose a viscoelastic flow model where the foam microstructure is represented from suitable microstructure descriptors whose evolution is governed by the macroscopic flow kinematics.

Introduction

Aqueous foams are concentrated dispersions of gas bubbles in a surfactant solution. Their structures are organized over a large range of length scales and complex flows take place at different scales [START_REF] Cohen-Addad | Flow in foams and flowing foams[END_REF].

The proposal of macroscopic constitutive equations allows for the efficient modelling and simulation of industrial processes involving the flow of foams [START_REF] Cheddadi | Numerical modelling of foam Couette flows[END_REF][START_REF] Benito | An elasto-visco-plastic model for immortal foams or emulsions[END_REF]. Usually, such descriptions remain however too phenomenological, and even though they predict accurately the flow kinematics, microstructure information remains often unaccessible. On the opposite side, fully microscopic simulations allow for very detailed descriptions of the foam microstructural evolution [START_REF] Bikard | Numerical prediction of the foam structure of polymeric materials by direct 3D simulation of their expansion by chemical reaction based on a multidomain method[END_REF]. However such approaches fail to address scenarios of industrial interest that usually involve the flow of foams in very large and complex 3D geometries.

The macroscopic flow model is expected to depend on the cellular structure: cell size, shape and orientation, as well as on the fluid rheology and the surface tension. Moreover, cell shape and orientation are induced by the flow. This microscopic information could be introduced into a macroscopic flow model by using standard upscaling and homogenization techniques. Thus inspired by [START_REF] Feyel | A multilevel finite element method (FE2) to describe the response of highly nonlinear structures using generalized continua[END_REF], at some locations in the domain in which an effective homogeneous fluid flows, we could attach a representative volume containing several cells, whose size, shape and orientation depend on the considered location. Now, a detailed microscopic calculation could be carried out in order to determine the effective fluid rheology. However, such a route, widely and successfully considered in a variety of fields, remains expensive from a computational point of view despite some attempts at combining it with advanced model reduction techniques [START_REF] Lamari | Routes for efficient computational homogenization of nonlinear materials using the Proper Generalized Decomposition[END_REF].

The most appealing description consists of a macroscopic flow model making use of some conformational variables describing the main microstructural features, as widely considered in the field of multiscale polymer modelling [START_REF] Keunings | Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory[END_REF][START_REF] Binetruy | Flows in polymers, reinforced polymers and composites. A multi-scale approach[END_REF]. Thus, in [START_REF] Karimi | PUFoam: A novel open-source CFD solver for the simulation of polyurethane foams[END_REF] the authors study flows during foaming, considering cell evolutions but without addressing the shape orientation and then, the induced anisotropy. Richer microstructure descriptions can be obtained by using a set of configurational coordinates, from which a conformation tensor can be derived and a macroscopic constitutive equation established. This route was successfully considered in [START_REF] Tlili | Mechanical formalisms for tissue dynamics[END_REF]. In the present work, we propose an alternative simpler microstructure description and its coupling with the macroscopic flow.

Remark. In the sequel, we consider the following tensor products, where Einstein's summation convention is assumed:

-if

Cell conformation

A very simple description of a 3D cell consists of a deformable ellipsoid (with constant volume) ranging from the spherical shape to the infinite aspect ratio ellipsoid (rod). Such an ellipsoid could be represented by means of three orthogonal extensible springs, with reference length 2L 0 and stiffness K.I nt h e sequel, we restrict our analysis to 2D scenarios (by considering ellipses instead of ellipsoids) but the derived models and their numerical solution procedures can be straightforwardly extended to 3D. First, we consider the kinematics of a single linear elastic dumbbell as starting point for elaborating the cell conformation.

Kinematics of an extensible rod

The extensible rod, of reference length 2L 0 and assumed aligned in direction p (p having a unit norm), is represented by an elastic spring of length 2L (in the deformed state) and stiffness K equipped with two beads at its extremities where hydrodynamic forces act. In the sequel, the word hydrodynamic refers to the viscous drag force and not to the one considered in other works to describe the effects of a bead kinematics on the others from the use of the Oseen tensor. These forces scale with the fluid -bead relative velocity, the former given by v 0 + ∇v • pL and the latter by v G + ṗL + p L,w h e r ev 0 is the unperturbed fluid velocity at the rod center of gravity and v G the velocity of the rod centre of gravity. A sketch of the rod and the forces acting on it is depicted in Fig. 1.

The system is assumed inertialess, that implies the equilibrium of forces and torques. The first implies F H (pL)+F H (-pL)=0, leading to v 0 = v G , that is, the rod centre of gravity moves with the fluid. Now, to prevent a resultant torque, force F H (pL) must align with p, i.e.

F H (pL)=λp, λ ∈ R.T h u s ,w eh a v e F H (pL)=ξ(∇v • pL -ṗL -p L)=λp, (1)
that multiplying by p and taking into account that p • p = 1 and consequently p • ṗ = 0, yields expression that introduced into Eq. ( 1) reads

ξ(∇v :(p ⊗ p)L -L)=λ, (2)
ξ(∇v • pL -ṗL -p L)=ξ(∇v :(p ⊗ p)L -L)p, (3)
that leads to the rotary velocity ṗ

ṗ = ∇v • p -∇v :(p ⊗ p)p, (4)
that is nothing else than the standard Jeffery expression for ellipsoids of infinite aspect ratio (rods) [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF]. Now, by equating the force acting on the beads λ with the one within the spring, we have

2K(L -L 0 )=ξ(∇v :(p ⊗ p)L -L), (5) or
L = - 2K ξ (L -L 0 )+∇v :(p ⊗ p)L). (6)
Thus, the kinematics of an elastic dumbbell of reference length 2L 0 with conformation at time t given by its orientation p and length 2L,r e a d

ṗ = ∇v • p -∇v :(p ⊗ p)p L = -2K ξ (L -L 0 )+∇v :(p ⊗ p)L) . ( 7
)
2.2 From rigid ellipses to orthogonal elastic bi-dumbells

In [START_REF] Chinesta | From single-scale to two-scales kinetic theory descriptions of rods suspensions[END_REF] it was proven that in order to represent a rigid ellipse whose kinematics are given by the Jeffery equation it suffices to consider a rigid system composed of two mutually orthogonal rods whose lengths correspond with the length of the ellipse axes.

In this case, if p refers to the direction of the ellipse largest axis, and

F = r 2 -1 r 2 +1
,w i t hr the ellipse aspect ratio, we have ṗ

= Ω • p + F (D • p -D :(p ⊗ p)p) , (8)
where Ω and D are respectively the vorticity and the rate of strain tensors, 2Ω = ∇v -(∇v) T and 2D = ∇v +(∇v) T . Now, we address the more general case in which both rigid rods are replaced by two extensible and mutually perpendicular springs of reference lengths L 0 1 and L 0 2 . In the sequel, the same reference lengths are assumed for both dumbbells, i.e. 2L 0 1 =2L 0 2 =2L 0 . In the present configuration, and considering that as proven in our former works the centre of gravity moves with the fluid, the hydrodynamic forces applying at beads L 1 p 1 and L 2 p 2 , F H 1 and F H 2 read respectively

F H 1 = ξ(∇v • p 1 L 1 -ṗ1 L 1 -p 1 L1 ), (9)
and

F H 2 = ξ(∇v • p 2 L 2 -ṗ2 L 2 -p 2 L2 ), (10)
with p 1 ⊥ p 2 , and with their orientation rates of change expressed from

ṗ1 = ω × p 1 ṗ2 = ω × p 2 . (11)
The angular momentum balance implies now

L 2 1 p 1 × (∇v • p 1 -ṗ1 )+L 2 2 p 2 × (∇v • p 2 -ṗ2 )=0, (12)
which coincides with the expression obtained in the case of rigid rods [START_REF] Chinesta | From single-scale to two-scales kinetic theory descriptions of rods suspensions[END_REF], proving the validity of the Jeffery equation in the case of orthogonal elastic bi-dumbells.

Introducing the Jeffery equation (8

) with r = L 1 L 2 that implies F = L 2 1 -L 2 2 L 2 1 +L 2 2 , i.e. ṗ1 = Ω • p 1 + F D • p 1 -p T 1 • D • p 1 p 1 , (13)
into the expression of the hydrodynamic force acting on bead p 1 L 1 ,w eh a v e

F H 1 = ξ(∇v • p 1 L 1 -ṗ1 L 1 -p 1 L1 ) (14) = ξ ∇v • p 1 L 1 -Ω • p 1 L 1 -F D • p 1 L 1 -p T 1 • D • p 1 p 1 L 1 -p 1 L1 (15) = ξL 1 (1 -F)D • p 1 + F p T 1 • D • p 1 p 1 -ξp 1 L1 . ( 16
)

Modelling incompressible ellipses from orthogonal elastic bi-dumbells

The projection of force F H 1 in the direction p 1 is the one that causes the spring extension, i.e.

2K(L 1 -L 0 )+F I 1 = ξL 1 p T 1 • D • p 1 -ξ L1 , (17)
where the force F I 1 ensures the incompressibility constraint. Similar calculations lead to

2K(L 2 -L 0 )+F I 2 = ξL 2 p T 2 • D • p 2 -ξ L2 , (18)
where it was assumed that

K 1 = K 2 = K.
Thus, the extension velocities read

L1 = - 2K ξ (L 1 -L 0 )+L 1 p T 1 • D • p 1 - 1 ξ F I 1 , (19) and L2
= - 2K ξ (L 2 -L 0 )+L 2 p T 2 • D • p 2 - 1 ξ F I 2 . ( 20
) T a k i n gi n t oa c c o u n t p T 1 • D • p 1 = D :(p 1 ⊗ p 1 ) p T 2 • D • p 2 = D :(p 2 ⊗ p 2 ) , (21)
and the fact that, since p 1 and p 2 are mutually orthogonal, (p 1 ⊗ p 1 )+ (p 2 ⊗ p 2 )=I (with I the identity tensor), the incompressibility constraint

d dt (L 1 L 2 )=0reads L1 L 2 +L 1 L2 = - 2K ξ (L 1 -L 0 )L 2 - 2K ξ (L 2 -L 0 )L 1 -F I 1 L 2 ξ -F I 2 L 1 ξ =0,( 22
)
which simplifies to

L1 L 2 + L 1 L2 = - 4K ξ L 1 L 2 + 2K ξ L 0 (L 1 + L 2 ) -F I 1 L 2 ξ -F I 2 L 1 ξ =0. ( 23
)
Forces related to the incompressibility constraint are expected to contribute isotropically to the resulting macroscopic stress, and then taken in the flow incompressibility constraint. Thus, if we compute the contribution of F I 1 and F I 2 to the stress by using the Kramers rule, we have

σ I = F I 1 p 1 ⊗ L 1 p 1 + F I 2 p 2 ⊗ L 2 p 2 = F I 1 L 1 (p 1 ⊗ p 1 )+F I 2 L 2 (p 2 ⊗ p 2 ),(24) which suggests considering F I 1 L 1 = F I 2 L 2 ,s i n c ep 1 ⊗ p 1 + p 2 ⊗ p 2 = I. Thus, considering F I 2 = F I 1 L 1
L 2 in Eq. ( 23), we have

F I 1 L 2 1 + L 2 2 ξL 2 = 4K ξ L 1 L 2 - 2K ξ L 0 (L 1 + L 2 ), (25)
or

F I 1 =4K L 1 L 2 2 L 2 1 + L 2 2 -2K L 0 (L 1 + L 2 )L 2 L 2 1 + L 2 2 , (26)
that vanishes for the relaxed case

L 1 = L 2 = L 0 , i.e. F I 1 (L 1 = L 2 = L 0 )=0 and consequently F I 2 (L 1 = L 2 = L 0 )=0.
Thus, finally the governing equations for the orthogonal elastic bi-dumbbell representing an incompressible ellipse read:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ṗ1 = Ω • p 1 + L 2 1 -L 2 2 L 2 1 +L 2 2 D • p 1 - L 2 1 -L 2 2 L 2 1 +L 2 2 p T 1 • D • p 1 p 1 L1 = -2K ξ (L 1 -L 0 ) -1 ξ F I 1 + L 1 p T 1 • D • p 1 L2 = -L 2 L 1 L1 F I 1 =4K L 1 L 2 2 L 2 1 +L 2 2 -2K L 0 (L 1 +L 2 )L 2 L 2 1 +L 2 2 F I 2 = F I 1 L 1 L 2 . ( 27
)

Conformation descriptor

When considering control volumes (small enough with respect to the flow but large enough with respect to the cell size) in a flowing foam, it can be observed that the cells in each control volume have similar shapes and orientations. In that case, the cell population in each volume element can be described using p 1 and L 1 (when L 0 is assumed known).

The contribution of a cell to the stress using the Kramers rule can be obtained from the elastic forces. The contribution due to the incompressibility constraint being isotropic, it can be aggregated to the pressure term. Thus, the conformation contribution σ c results

σ c =2K(L 1 -L 0 )p 1 ⊗ L 1 p 1 +2K(L 2 -L 0 )p 2 ⊗ L 2 p 2 (28) =2K(L 1 -L 0 )L 1 (p 1 ⊗ p 1 )+2K(L 2 -L 0 )L 2 (p 2 ⊗ p 2 ) (29) =2K (ΔL 1 L 1 (p 1 ⊗ p 1 )+ΔL 2 L 2 (p 2 ⊗ p 2 )) , (30)
that vanishes in the relaxed configuration L 1 = L 2 = L 0 ,w i t hΔL 1 = 0 and ΔL 2 =0. In these circumstances, the simplest choice for the conformation tensor consists of the second-order symmetric tensor c defined from

c = ΔL 1 L 1 (p 1 ⊗ p 1 )+ΔL 2 L 2 (p 2 ⊗ p 2 ). (31)
It is important to note that the conformation does not involve an averaging process affecting the different cells inside the control volume because we assumed that all of them share almost the same conformation. However, as it is usual in the modelling of suspensions, this hypothesis could easily be relaxed and we would thus consider the average of different cell conformation, as addressed in the micro-macro simulations by Keunings and coworkers using the Lagrangian Particle Method -LPM - [START_REF] Halin | The Lagrangian particle method for macroscopic and micro-macro viscoelastic flow computations[END_REF][START_REF] Keunings | Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory[END_REF].

Even though many rheological behaviours could be associated to flowing foams, this paper focuses on the consideration of induced anisotropy due to the cell deformation. For this reason, and without loss of generality, we consider in what follows the simplest rheological behaviour consisting of a Newtonian behaviour complemented with an elastic contribution related to the cell deformation. In any case, more complex bulk rheologies could be considered, as for example the one related to a viscoplastic behaviour.

Thus, the macroscopic viscoelastic constitutive equation reads

σ = -pI +2ηD + μc, ( 32
)
where p is the pressure, that can be viewed as the Lagrange multiplier associated with the macroscopic flow incompressibility constraints, η the effective homogenized fluid viscosity, and μ the rheological parameter affecting the contribution of the microscopic conformation that scales with K and the volume concentration of cells. The origin of the elasticity introduced into the conformation field evolution is related to the surface tension that resists the cell deformation from its spherical reference conformation. Topological changes can operate at the elemental cell level. However, on average their net effect is to avoid too large cellular distorsions. Thus, the coefficient μ should describe both the surface tension as well as the microscopic topological changes. In the proposed model, it remains purely phenomenological and should be identified from appropriate rheological tests.

Even though the present model does not address relaxation mechanisms, they could easily be incorporated by including a viscous component at the spring-beads level.

Rheological behaviour

In order to study the rheological response of the proposed model, we carried out a numerical study to obtain the loss and storage moduli of a foam. To proceed, consider a periodic shear strain

γ(t)=γ 0 sin(ωt) (33)
whose time derivative provides the evolution of the shear rate

γ(t)=γ 0 ω cos(ωt). (34)
In simple shear flow, the velocity gradient thus reads ∇v = 0γ(t) 00 ,a n dw e could use the model equations ( 27) to obtain the time evolution of the shape and orientation of a cell subject to such an oscillatory flow. We can thus obtain the time evolution c(t) of the conformation tensor introduced in the previous section.

We can now decompose the off-diagonal component of tensor c(t), that is c 12 (t) in its in-phase and out-of-phase contributions according to

c 12 (t)=c ′ 12 sin(ωt)+c ′′ 12 cos(ωt), (35)
where c ′ 12 and c ′′ 12 are analogous to the so-called storage (G ′ )a n dl o s s( G ′′ ) modulus respectively (they actually correspond to the storage and loss moduli associated with the stress contribution σ c ).

Figure 2 shows the value of c ′ 12 and c ′′ 12 for an oscillation angular frequency ranging from 0.1 to 100 rad/s. The parameters are ξ =0.1, K =1,L 0 =0.005 (metric system of units). We can recognize a Maxwell-like model, with the classical slopes of 1 (G ′′ )a n d2( G ′ ) at low frequency.

To obtain properly the storage (G ′ )a n dl o s s( G ′′ ) moduli, we should also consider the viscous component of the stress tensor. Thus the moduli are given by Finally, we also conducted another rheological test to observe the behaviour of the conformation tensor when a constant shear flow is brutally stopped. As expected, we obtain an exponential relaxation of the tensor components towards zero.

G ′ = σ ′ 12 γ 0 = μ c ′ 12 γ 0 (36)
G ′′ = σ ′′ γ 0 =2ηω + μ c ′′ 12 γ 0 . (37)

Macroscopic flow problem

Neglecting inertia, the macroscopic flow problem is defined as follows:

⎧ ⎨ ⎩ ∇•σ = 0 σ = -pI +2ηD + μc ∇•v =0 , ( 38
)
where σ is the Cauchy's stress tensor, I the unit tensor and v the macroscopic fluid velocity field. It can be noticed that the problem reduces to the standard Stokes flow model as soon as μ = 0. Incompressibility is justified by the fact that drainage is neglected.

The flow model must be complemented with the conformation expression

c = ΔL 1 L 1 (p 1 ⊗ p 1 )+ΔL 2 L 2 (p 2 ⊗ p 2 ), (39)
whose evolution is governed by the microstructural model

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ṗ1 = Ω • p 1 + L 2 1 -L 2 2 L 2 1 +L 2 2 D • p 1 - L 2 1 -L 2 2 L 2 1 +L 2 2 p T 1 • D • p 1 p 1 L1 = -2K ξ (L 1 -L 0 ) -1 ξ F I + L 1 p T 1 • D • p 1 L2 = -L 2 L 1 L1 F I 1 =4K L 1 L 2 2 L 2 1 +L 2 2 -2K L 0 (L 1 +L 2 )L 2 L 2 1 +L 2 2 F I 2 = F I 1 L 1 L 2 . ( 40
)
In order to solve the resulting flow model defined in Ω ⊂ R 2 , appropriate boundary conditions must be enforced at the domain boundary

Γ ≡ ∂Ω v = v g in Γ D t = σ • n = t g in Γ N , (41)
with

Γ D ∪ Γ N = Γ and Γ D ∩ Γ N = ∅.

Numerical solution strategy

Cells are assumed represented at the initial time by a spherical conformation tensor c = 0, since the cells are initially in the relaxed state,

L 1 = L 2 = L 0 .
1. The associated Stokes problem is solved at time t n , from the microstructural term given at the previous time step μc(t n-1 ). The problem is solved using a standard mixed velocity-pressure formulation using any finite element satisfying the stability conditions, the so-called LBB conditions. In the example reported later, we considered Q9/Q4 finite elements; 2. Then, from the just computed velocity field v(t n ) the orientation and extension fields, p 1 (x,t n )a n dL 1 (x,t n ) respectively, are updated by using a first-order discontinuous Galerkin formulation; 3. Finally, the conformation tensor is updated c(t n ).

These three steps are repeated until reaching the maximum simulation time or the steady state.

Updating the conformation field

In the numerical experiments described below, we consider 2D flows defined in 2D geometries Ω ⊂ R 2 . In the 2D case, the unit vector p 1 can be expressed by p 1 =( c o sθ, sin θ) T . By taking the time derivative, we get ṗ1 = θ(sin θ, cos θ) T .U s i n gθ as orientation descriptor, the Jeffery equation can be expressed as θ(sin θ, cos θ)

T = G(θ, ∇v,L 1 ). ( 42
)
Multiplying the previous expression by (-sin θ, cos θ), we obtain the scalar equation θ = G(θ, ∇v,L 1 ),

where the material derivative can be expressed in an Eulerian framework by introducing the orientation field θ(x,t) whose evolution is governed by

∂θ ∂t + v •∇θ = G(θ, ∇v,L 1 ). ( 44
)
The scalar equation governing the evolution of the spring length can also be written as

∂L 1 ∂t + v •∇L 1 = H(θ, ∇v,L 1 ). ( 45
)
Both equations ( 44) and ( 45) are purely advective, and thus appropriate discretization taking into account their hyperbolic character must be used. We make here the simplest choice, a first-order discontinuous Galerkin scheme, that considers the generic variable P (θ or L 1 ) constant in each element Ω e belonging to the mesh M of Ω. The test function is assumed also constant in each element, vanishing outside. The balance in element Ω e reads, taking into account the flow incompressibility (∇•v =0):

Ω e ∂P ∂t dx + Ω e ∇•(vP)dx = Ω e J dx, (46)
where the source term J represents G or H depending on the considered equation.

Using the divergence theorem, the second term of the left-hand side can be written from the boundary flux, i.e.

Ω e ∂P ∂t dx + ∂Ω e Pv • n dx = Ω e J dx, ( 47
)
where n is the unit outward vector normal to the element boundary.

As P is not defined on ∂Ω e , we consider the element boundary decomposition ∂Ω e = ∂ + Ω e ∪ ∂ -Ω e ,w h e r e∂ -Ω e and ∂ + Ω e represents the inflow and outflow element boundaries, both defined from v • n < 0a n dv • n > 0r espectively. Then we assume that the property on the inflow boundary is given by its value at the upstream element, i.e. P(x ∈ ∂ -Ω e )=P e -, whereas on the outflow element boundary, it is given by the property at element Ω e , i.e. P(x ∈ ∂ + Ω e )=P e .

Discussion

-In the flow model just proposed there is no size effects, which implies that the characteristic length of cells must be small in relation to the one characterizing the spatial variation of the macroscopic velocity field. -The mesh considered for integrating the velocity field must be small enough to capture all the macroscopic velocity field details. -The mesh considered for calculating the conformation field must be small enough for assuming cells in each element of the mesh described by the conformation field inside the element, and it must be large enough for assuming that it represents a population of cells. However, as soon as the model is described in a continuous way, there is not contradiction with the fact of considering elements smaller that the characteristic size of the cell (in the case of polymer flows one can consider elements smaller that the size of a molecule). 5 Numerical results

Uncoupled microstructure-flow calculations

In this section, we address some simple flows in order to evaluate the response of the cells. The model is uncoupled in the sense that the flow induces cell deformation, but the cell conformation does not affect the flow kinematics. Thus, even though some of the flows addressed here could exhibit rich kinematics in practice, in what follows the uncoupled solution does not allow to capture such a rich kinematics, as for example shear banding. Moreover, the only interest of this section being the evolution of the cell deformation and its macroscopic description, rheology is not considered in the analyses carried out.

In the numerical examples considered in this section a unit square is considered as depicted in Fig. 3. The initial (relaxed) conformation is given by a zero conformation tensor c(x,t =0)=0.

Four different flows are considered: (i) a simple shear; (ii) a contraction flow; (iii) the driven cavity flow problem and finally (iv) the flow around a square obstacle. In all cases, the microstructure was computed on the basis of the associated Stokes kinematics, i.e. considering the uncoupled flow associated with μ =0.

In order to quantify the way the flow behaviour (shear, rigid motion or elongation) affects the microstructure evolution we consider a flow criterion. For that purpose, first, we introduce the relative rate of rotation W from

W = ω -w, (48)
where ω is related to the flow vorticity Ω according to where ǫ is the third-order permutation tensor (also known as the Levi-Civita tensor). The vector ω c a na l s ob ew r i t t e ni nt e r m so ft h ec u r lo ft h ev e l o c i t y ,

Ω = ǫ • ω, (49)
ω = - 1 2 ∇×v. (50)
The vector w in Eq. ( 48) represents the angular velocity of the eigenvectors of the rate of strain tensor D.

A simple local descriptor of the type of flow can then be constructed from the second invariant of D,γ = √ 2D : D, and from the norm of W, W , according to

χ = 4 W γ +2 W . (51)
We have 0 ≤ χ ≤ 2, and more specifically χ = ⎧ ⎨ ⎩ 0 in planar extension, 1 in pure shear, 2 in rigid motion.

(52)

Figure 4 shows the evolution of the microstructure conformation for different K in a pure shear with ξ =0 .1a n dL 0 =0 .005 (metric system of units). The qualitative analysis performed here does not require a precise determination of these coefficients. As it can be noticed, the more stiffness the spring has, the less deformation the conformation presents. In order to prove that the conformation is not only accommodating a deformation but also that it is rotating, we included a small cercle on the ellipsoid surface to appreciate the way in which that point is evolving in time.

Simple shear flow

We consider V =2m• s -1 on the N-boundary, zero velocity on the S-boundary and a linear velocity evolution on the E-boundary and W-boundaries. The microstructure effects were introduced by using K =1.

Figure 5 depicts the velocity field, that evolves almost linearly through the domain thickness, as well as the flow criterion, that as expected corresponds to a perfect shear behaviour. The pressure field is constant in the whole flow domain.

Extrusion-like flow

In the present case, we consider again K = 1. A Poiseuille velocity profile (parabolic) is enforced on the W-boundary (-2y(y -1), 0) T , whereas at the Eboundary the fluid leaves the square domain throughout an exit where tensionfree boundary conditions are enforced.

Figure 7 depicts both components of the velocity field whereas Fig. 8 shows the steady state conformation and the flow criterion. As it can be seen, the largest conformation axis remains aligned with the streamlines. A pure extension is noticed along the symmetry axis in agreement with the flow criterion.

Driven cavity flow

The only difference with respect to the previously analyzed flow is that now a unit horizontal velocity is applied on the top wall (N-boundary) that induces the fluid flow within the cavity.

Figure 9 depicts both components of the velocity field. Figure 10 shows the conformation at a given time as well te steady state flow crieria. The conformation evolves periodically and no steady state is reached.

Flow around an obstacle

In the present case, a unit horizontal velocity is enforced on the W-boundary, with non-slip conditions on the top and bottom walls and a free boundary condition is specified on the E-boundary. The domain contains in its center an obstacle where non-slip boundary conditions are enforced. Figure 12 depicts the conformation and flow criteria. As it can be noticed, elongation is specially present upstream and downstream, the shear being located in the neighbourhood of the upper and lower sides of the obstacle, with the expected effects on the microstructure conformation.

Coupled simulations

This section adresses a coupled simulation in which flow kinematics induces microstructure evolution and the latter affects at its turn the flow kinematics.

Figure 13 shows two snapshots taken from films of flowing foams where the flow and microstructures were assumed almost at steady-state. The microstructure of these snapshots will be used first to identify the model parameters and then to test the agreement between the computed and observed microstructures. Of course, from the information of kinematics only, we cannot conclude on rheological aspects but at least this serves to evaluate the proposed model in terms of conformation evolution.

In order to extract from these images the conformation tensor, both images were segmented (every pixel is transformed into either black or white) and then cell boundaries were easily identified. From that, the centre of gravity of each cell can easily be obtained as well as the conformation tensor. If X e i , i =1 ,...,Q e are the points defining the wall of cell C e , the centre of gravity and the inertia tensor J e are obtained from

X e G = 1 Q e Q e i=1 X e i ( 53
)
and

J e = 1 Q e Q e i=1
(X e i -X e G ) ⊗ (X e i -X e G ) (54) respectively. The numerical inertia tensor can be calculated as soon as p 1 and L 1 are known (p 2 and L 2 derive from p 1 and L 1 ).

J = L 2 1 (p 1 ⊗ p 1 )+L 2 2 (p 2 ⊗ p 2 ). (55)
In Fig. 14, the experimental ellipses related to J e are superimposed to the cells, and also to the solution predicted by using the proposed model. The model parameters K, ξ, η and μ are, as previously indicated, chosen in order to obtain cell shapes close to those observed experimentally: (metric system of units) K =1 ,ξ =0 .1, η = 1 and μ =0 .1. L 0 =0 .005 was identified from the average cell surface observed experimentally. In absence of rheological data, we noticed that many choices of those parameters lead to similar microstructures and that the impact of those on the kinematics was almost negligible. Thus, a proper rheological characterization seems compulsory to obtain an adequate flow-microstructure coupling.

Conclusions

This work proposes a simple model for flowing foams, where the microstructure is introduced from a conformation tensor that describes the structural elasticity. The numerical predictions have been compared with some experimental results proving the ability of the model to describe the effective kinematics as well as the flow induced microstructure evolution.

More applicative analyses, in particular in the study of industrial processes, require an appropriate rheological characterization and very probably the proposal of appropriate rheometric devices. All these aspects will be addressed in future works.

Compliance with Ethical Standards. The authors declare that they have no conflict of interest.

  a and b are first-order tensors, the single contraction • reads (a•b)=a j b j (); -if a and b are first-order tensors, the dyadic product ⊗ reads (a ⊗ b) jk = a j b k ; -if a and b are respectively second and first-order tensors, the single contraction • reads (a • b) j = a jm b m ; -if a and b are second-order tensors, the single contraction • reads (a•b) jk = a jm b mk ; -if a and b are second-order tensors, the double contraction : reads (a : b)=a jk b kj .

Fig. 1

 1 Fig. 1 Extensible rod immersed in a flow.

5 Fig. 2

 52 Fig. 2 Oscillatory shearc ′ 12 and c ′′ 12 components

Fig. 3

 3 Fig. 3 2D domain.

Fig. 4

 4 Fig. 4 Conformation evolution for different spring stiffness in a simple shear flow.

Fig. 5

 5 Fig. 5 Velocity field (left) and flow criterion (right) in a simple shear flow.

Fig. 6

 6 Fig. 6 Microstructure conformation in a simple shear flow.

Figure 6

 6 Figure6shows the steady state conformation. Because of the shear, the conformation is expected to rotate clockwise. The microstructure shows a significant variation along the domain as a consequence of the different velocity and constant shear rate.

Fig. 7

 7 Fig. 7 Extrusion-like flow problem. Velocity field: (left) x-component and (right) ycomponent.

Fig. 8

 8 Fig. 8 Conformation distribution (left) and flow criterion (right) in the extrusion-like flow.

Fig. 9

 9 Fig. 9 Driven cavity flow problem. Velocity field: (left) x-components and (right) ycomponent.

Fig. 10

 10 Fig. 10 Conformation (left) and flow criterion (right) in the driven cavity flow problem.

Fig. 11

 11 Fig. 11 Flow around an obstacle. Velocity field: (left) x-component and (right) ycomponents.

Figure 11

 11 Figure 11 depicts both components of the velocity field. The velocity field presents two stagnation points located at the intersection between y =0 .5 and the obstacle. These stagnation points correspond with the maximum and minimum pressures.

Fig. 12

 12 Fig. 12 Conformation (left) and flow criterion (right) in the flow around an obstacle.

Fig. 13

 13 Fig. 13 Experimental snapshots (courtesy of F. Graner).

Fig. 14

 14 Fig. 14 Experimental microstructure with the associated conformation superimposed (left) and the associated numerical predictions (right).

Aknowledgements. This project has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 675919.