Nicolas Dupin

Rémi Parize

El-Ghazali Talbi

Matheuristics to stabilize column generation: application to a technician routing problem

This paper considers a simplified Technician Routing and Scheduling Problem with skill constraints, where a Column Generation (CG) heuristic was shown effective. This work proposes new CG schemes to stabilize the CG process and thus to accelerate the CG heuristics. Solving CG subproblems with matheuristics allows to diversify the column generation. Furthermore, a tabu search matheuristic allows to generate aggressively columns at each iteration. Both techniques imply a better stability of the CG scheme: the diversification is interesting for the first iterations whereas tabu intensification is especially useful for the last iterations. It implies a significant acceleration of the CG convergence. A perspective of these CG strategies proposed is an extension to other problems where CG induces independent and heterogeneous subproblems.

Introduction

The technician routing and scheduling (TRS) problem was defined for the Challenge ROADEF [START_REF] Dutot | Technicians and interventions scheduling for telecommunications[END_REF]. The optimization assigns technicians to vehicles to serve customers' requests, with constraints for skills and number of technicians required for the interventions as well as capacity constraints to use tools. Matheuristics were used to tackle this problem, using Mixed Integer Programming (MIP) for subproblems [START_REF] Fırat | An improved MIP-based approach for a multi-skill workforce scheduling problem[END_REF][START_REF] Kovacs | Adaptive large neighborhood search for service technician routing and scheduling problems[END_REF][START_REF] Pokutta | France telecom workforce scheduling problem: a challenge[END_REF]. This paper considers a relaxation, a multiple depot vehicle routing problems with time windows (VRPTW) where skill constraints define the jobs that vehicles can process. This relaxation was already studied in [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF], analyzing the impact of the skill constraints in highly constrained instances and providing tailored matheuristic approaches for such instances.

Column Generation (CG, see [START_REF] Lübbecke | Selected topics in column generation[END_REF]) approaches are known to be efficient to solve VRPTW to optimality [START_REF] Feillet | A tutorial on column generation and branch-and-price for vehicle routing problems[END_REF]. Facing large instances, a CG scheme can give rise to primal heuristics [START_REF] Archetti | A survey on matheuristics for routing problems[END_REF][START_REF] Taillard | A heuristic column generation method for the heterogeneous fleet VRP[END_REF]. A bottleneck of CG algorithms is the difficulty to have a stable convergence, with "bang-bang" effects and erratic dual variables between consecutive iterations. To stabilize the convergence of dual variables and to reduce the computational effort, stabilization techniques are essential. In that goal, mathematical properties can be exploited as in [START_REF] Amor | On the choice of explicit stabilizing terms in column generation[END_REF][START_REF] Du Merle | Stabilized column generation[END_REF][START_REF] Rousseau | Interior point stabilization for column generation[END_REF], but also heuristics as in [START_REF] Hansen | Variable neighborhood search: Principles and applications[END_REF].

The contribution of this paper is to investigate matheuristic techniques for the stabilization of the CG algorithm in [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF], to accelerate both CG computations of dual bounds and CG heuristics. Similarly to [START_REF] Dupin | Dual Heuristics and New Lower Bounds for the Challenge EURO/ROADEF 2010[END_REF], matheuristics are not only techniques using exact methods within a meta-heuristic, but also techniques increasing the size of instances where dual bounds for MIPs are computable.

This paper is organized as follows. In section 2, we introduce the applicative technician routing and scheduling problem and its compact MIP formulation. In section 3, we present the extended MIP formulations of the problem and its standard CG algorithm. In section 4, an alternative CG scheme is proposed to enhance the diversification among generated columns. In section 5, a matheuristic algorithm is proposed to generate intensively columns for CG subproblems. In section 6, the computational results are discussed. In section 7, our contributions are summarized, discussing also future directions of research.

Problem statement

In the following simplification of [START_REF] Dutot | Technicians and interventions scheduling for telecommunications[END_REF], we assume that the jobs are realized by a single technician. In this case, the skill constraints are binary: either a vehicle can process a job, or it cannot. Furthermore, we relax scheduling and capacity constraints from [START_REF] Dutot | Technicians and interventions scheduling for telecommunications[END_REF], focusing on the Technician Routing problem.

The routing optimization concerns a set of technicians, who have to realize jobs at customers' places with time window constraints on the arrival time. The jobs require specific skills, which limits the possible routes of the vehicles. The optimization horizon is one day. The objective function to minimize is the total distance of the routes, adding high penalties for non allocated jobs. We note that a similar relaxation of [START_REF] Dutot | Technicians and interventions scheduling for telecommunications[END_REF] with additional constraints was also studied with a matheuristic approach in [START_REF] Pillac | A parallel matheuristic for the technician routing and scheduling problem[END_REF]. Notations are gathered in the Table 1. Distance between the location of jobs j and j .

D(i, j)

Distance between the depot di to the location of job j. T (j, j )

Transportation time between the location of jobs j and j . T (i, j)

Transportation time between the depot di to the location of job j. A first MIP formulation can be written similarly to [START_REF] Agra | The robust vehicle routing problem with time windows[END_REF]. Binaries u i j,j ∈ {0, 1} are defined for all technician i ∈ I and all jobs (j, j ) ∈ J 2 i related to i , where u i j,j = 1 indicates that technician i realizes job j and then job j . By convention u i j,j = 0. We extend these notations with binaries for all technicians i ∈ I and all jobs j ∈ J i with u i di,j (and u i j,di respectively) indicating if job j is the first (respectively the last) job of technician i, and u i di,j = 0 if i = i . Furthermore, continuous variables t j represent the start time for job j, bounded with the time-windows constraints: for all j ∈ J , tmin

j t j tmax j . min u i j,j
,t j j,j ∈J ∪I D(j, j )u i j,j

+ j∈J i Pj 1 -i,j u i j ,j (1) ∀j
, i∈I j ∈J i ∪{d i } u i j j 1 ( 2
)
∀i ∈ I, j,

j ∈J i ∪{d i } u i j j = j ∈J i ∪{d i } u i jj (3) ∀i ∈ I, j, j ∈J i ∪{d i } u i j j 1 ( 4
)
∀i ∈ I,

j ∈J i u i j d i = j ∈J u i d i j 1 (5) ∀(j, j ) ∈ J 2 , tj + Dj + T (j, j ) t j + 1 -i u i j,j • M 1 j,j (6)
∀i ∈ I, j ∈ Ji,

tstart i + T (i, j) tj + 1 -u i d i j • M 2 i,j (7)
∀i ∈ I, j ∈ Ji, tj + Dj + T (j, i)

tend i + 1 -u i jd i • M 3 i,j (8)
∀i ∈ I, j, j u i j,j ∈ {0, 1}, tj ∈ [ tmin

j , tmax j ] (9)
The objective (1) to minimize is the distance of the total routing, with penalties for non attributed jobs. Constraints (2) are elementary constraints. Constraints (3) are flow constraints to model tours. Constraints (4) are capacity flow constraints, expressing that one job is done at most one time for an optimal solution. Indeed, having d(j, j ) + d(j , j ) d(j, j ), removing jobs planned twice decrease the objective function, there is an optimal solution allocating exactly one technician per job. Constraints (5) express that the depots d i are the flow source and arrival for all technician i. Constraints ( 6), ( 7) and ( 8) express time windows constraints, respectively between two consecutive jobs, and the respect of the technicians' starting and finishing time. Coefficients M 1 j,j , M 2 i,j , M 3 i,j can be chosen as following:

M 1 j,j = tmax j + D j + T (j, j ) -tmin j , M 2 i,j = tstart i + T (i, j) -tmin j and M 3 i,j = tmax j + D j + T (j, i) -tend i .
To increase the solving limits of this weak MIP formulation, matheuristics are efficient to tackle larger instances, even for highly-constrained instances, as shown in [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF]. MIP neighborhoods restrict the number of technician and the number of jobs in the MIP computations, for greedy algorithms constructing solutions from scratch, or for Variable Neighborhood Descent (VND) strategies.

Extended formulation and CG scheme

Similarly to [START_REF] Feillet | A tutorial on column generation and branch-and-price for vehicle routing problems[END_REF], we have an extended MIP formulation by enumerating all possible (and non empty) routes P i for each technician i. The cost of a route k ∈ P i is denoted C i,k . The variables of the MIP formulation are z i,k ∈ {0, 1}, such that z i,k = 1 if route k ∈ P i is chosen, and y j ∈ {0, 1} with y j = 1 if job j is not planned. It leads to the following MIP:

min y,z i∈I,k∈Pi C i,k z i,k + j∈J P j y j (10)
s.t : i,k:j∈k z i,k + y j 1 ∀j ∈ J , (11)
k∈Pi z i,k 1 ∀i ∈ I, (12)
z i,k , y j ∈ {0, 1} (13)
The difficulty is that P i cannot be enumerated. The Linear Programming (LP) relaxation is solved by the CG algorithm, adding iteratively new routes. Having a subset of routes C i ⊂ P i , the Restricted Master Problem (RMP) denotes the previous LP relaxation applied to a subset of variables:

min z i,k ,yj 0 i∈I,k∈Ci C i,k z i,k + j∈J P j y j s.t : i,k:j∈k z i,k + y j 1 (π) ∀j ∈ J , k∈Ci z i,k 1 (σ) ∀i ∈ I, (14)
π, σ denotes respectively the dual variables associated to the constraints ( 11), [START_REF] Hansen | Variable neighborhood search: Principles and applications[END_REF].

Variables must be added in the RMP if their reduced cost is strictly negative. If the remaining variables have a positive reduced cost, the RMP gives the LP relaxation of [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems[END_REF][START_REF] Fırat | An improved MIP-based approach for a multi-skill workforce scheduling problem[END_REF][START_REF] Hansen | Variable neighborhood search: Principles and applications[END_REF]. This negativity question among a non enumerable set is formulated as an optimization problem, minimizing the reduced costs of the columns that are not in ∪ i C i . These negativity question are decomposed for all technicians, RC * = min i RC * i , where RC * i is the minimization problem over the possible routes for technician i, RC * 0 ensures that the RMP gives the LP relaxation of [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems[END_REF][START_REF] Fırat | An improved MIP-based approach for a multi-skill workforce scheduling problem[END_REF][START_REF] Hansen | Variable neighborhood search: Principles and applications[END_REF]. Computations of RC * i are independent for the technicians, and can be formulated similarly to the MIP formulation of Section 2:

RC * i = minu σi + j,j (D(j, j ) -πj)u j,j s.t : ∀j ∈ J , j ∈J ∪{d i } u j j = j ∈J ∪{d i } u jj j ∈I u j j d i = j ∈I u j d i j 1 ∀(j, j ) ∈ J , tj + Dj + T (j, j ) t j + (1 -u j,j ) • M 3 i,j j ∈ J , tstart i + T (i, j) tj + (1 -u d i j ) • M 2 i,j j ∈ J , tj + Dj + T (j, i) tend i + (1 -u jd i ) • M 1 j,j ∀j, j u j j ∈ {0, 1}, tj ∈ [ tmin j , tmax j ] (15)
The LP relaxation of (10-12) is calculated by the Algorithm 1, with an initial set of columns as input. The initial set of columns can be empty, in this case the first RMP is feasible thanks to the penalty costs. Matheuristics can initialize the CG algorithm adding the columns corresponding to primal solutions.

To solve CG subproblems [START_REF] Lübbecke | Selected topics in column generation[END_REF] to optimality, a dynamic programming algorithm is valid formulating [START_REF] Lübbecke | Selected topics in column generation[END_REF] as a shortest path problem with resource constraints (ESPPRC), we refer to [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems[END_REF]. However, the computations to optimality are required only to prove the optimality of the RMP at the last iteration. For the first CG iterations, we just need to generate negative reduced cost solutions. The dynamic programming computations are time consuming in the Algorithm C set of initial columns. do :

solve RMP ( 14) with columns defined in C store dual variables σ and π and optimal cost from ( 14) for each technician i ∈ I : solve (15) to optimality with last (σ, π) values if CR * i < 0 then add the optimal column to C end for while : columns are added in C return the last cost of the RMP [START_REF] Lazic | Variable Neighbourhood Decomposition Search for 0-1 Mixed Integer Programs[END_REF] 1, it can be replaced by heuristics to generate quickly negative solutions. For instance, the matheuristics [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF] apply.

A general difficulty to implement CG algorithms is that the dual variables converge erratically. This is closely related to the properties of linear optimization as the optimums are extreme points. A consequence for the CG algorithm is that many iterations are necessary to have good dual variables to generate the most appropriate columns. Stabilization techniques enforce to smooth the convergence of dual variables to reduce the number of iterations to converge for the CG algorithm, and thus to reduce the total computation time for the CG convergence. To stabilize the convergence of dual variables and to reduce the computational effort, stabilization techniques are essential. Mathematical properties can be exploited as in [START_REF] Amor | On the choice of explicit stabilizing terms in column generation[END_REF][START_REF] Du Merle | Stabilized column generation[END_REF][START_REF] Rousseau | Interior point stabilization for column generation[END_REF], but also heuristics as in [START_REF] Hansen | Variable neighborhood search: Principles and applications[END_REF]. In the following, tailored stabilization schemes are proposed using matheuristics.

Algorithm 1 computes the LP relaxation of [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems[END_REF][START_REF] Fırat | An improved MIP-based approach for a multi-skill workforce scheduling problem[END_REF][START_REF] Hansen | Variable neighborhood search: Principles and applications[END_REF]. Such iterations, even truncated before the stopping criterion, can also be used to design heuristics. Integer computations of the RMP with the generated columns give primal solutions, which is efficient in the computational experiments of [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF]. Similarly to Lagrangian heuristics, the partial LP relaxations can be repaired to build primal solutions. The effort to accelerate the convergence of the Algorithm 1 has also an impact to accelerate CG heuristics.

POPMUSIC column generation stabilization

This section is motivated by the following fact: dual variables induce to generate columns with jobs with the highest values π j in [START_REF] Lübbecke | Selected topics in column generation[END_REF]. Hence, independent computations of (15) are likely to use the same jobs, with the highest values of π j . These generated columns are likely to be redundant in the recombination induced by the next RMP computation. It suggest to incorporate diversifications in the generation of columns. To achieve that goal, we investigate combined subproblems for a subset of technicians I 0 ⊂ I (typically two or three technicians), imposing a total diversification in these subproblems, with following MIP subproblems:

RC * I 0 = maxu i∈I 0 Ci (16.1) s.t : ∀j ∈ J i∈I 0 ,j ∈J (u i jj + u i jd i ) 1 (16.2) ∀i ∈ I0, -Ci +σi + j,j (D(j, j ) -πj)u i j,j (16.3) ∀i ∈ I0, ∀j ∈ J , j ∈J ∪{d i } u i j j = j ∈J ∪{d i } u i jj (16.4) ∀i ∈ I0, j ∈I u i j d i = j ∈I u i d i j 1 (16.5) ∀i ∈ I0, ∀(j, j ) ∈ J , tj + Dj + T (j, j ) t j + 1 -u i j,j • M 3 i,j (16.6) ∀i ∈ I0, j ∈ J , tstart i + T (i, j) tj + 1 -u i d i j • M 2 i,j (16.7) ∀i ∈ I0, j ∈ J , tj + Dj + T (j, i) tend i + 1 -u i jd i • M 1 j,j (16.8) ∀i ∈ I0, Ci 0 (16.9) ∀j, j u i j j ∈ {0, 1}, tj ∈ [ tmin j , tmax j ] (16)
This subproblem can be seen as a concatenation of the technician subproblems in I 0 ⊂ I, with (16.2) to enforce that each job can be affected to at least one technician in I 0 . Furthermore constraints (16.3) and objective function (16.1) induce an objective function where a technician without negative reduced cost counts for 0, and the global optimization minimizes the sum of the negative reduced cost among the technicians. It gives rise to the Algorithm 2, where the partitions of technicians are varying for each iterations, similarly to POPMUSIC (Partial optimization metaheuristic under special intensification conditions [START_REF] Taillard | POPMUSIC-Partial optimization metaheuristic under special intensification conditions[END_REF]).

Algorithm 2: POPMUSIC column generation algorithm Input: C set of initial columns. do :

solve RMP ( 14) with columns defined in C store dual variables σ and π and optimal cost from ( 14) compute PI a partition of I in small subsets for each subset I0 ∈ PI : solve ( 16) with a matheuristic with last (σ, π) values for each column c with a negative reduced cost add the column to C end for end for while : columns are added in C return the last cost of the RMP [START_REF] Lazic | Variable Neighbourhood Decomposition Search for 0-1 Mixed Integer Programs[END_REF] The subproblems ( 16) can be solved to search columns with a negative reduced cost using matheuristics, the number of technician is restricted and matheuristics restrict the size of the MIP computations removing jobs similarly to [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF]. With several technicians in the subproblems [START_REF] Pillac | A parallel matheuristic for the technician routing and scheduling problem[END_REF], the dynamic programming algorithm [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems[END_REF] does not apply anymore. The POPMUSIC CG scheme applies only to generate columns with a negative reduced cost. To prove the optimality of the RMP, Algorithm 1 applies for the remaining iterations. -I0 a subset of technician.

-s a cyclic permutation of I0 with order(s) = |I0|.

-σ, π the dual variables of the last RMP computation. Initialization: C = ∅, the columns to add in the RMP for each technician i ∈ I0 :

Let i = i0, J0 = J for k = 1 to |I0| : solve ( 15) for technician i with (σ, π) values and the remaining jobs in J0 if the solution induces a column with a negative reduced cost add the column in C remove the jobs of the column in J0 i = s(i0) end for : end for return C Algorithm 3 decomposes the search of columns with a negative reduced cost using only computations with one technician. These decomposed subproblems can also be solved as ESPPRC with [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems[END_REF]. This scheme allows to generate the best individual columns as in Algorithm 1, and also good complementary columns for these best individual columns. The solutions of [START_REF] Lübbecke | Selected topics in column generation[END_REF] with Algorithm 3 are uneven in the reduced costs following the order of priority of the technicians, whereas a straightforward matheuristic search applied to [START_REF] Lübbecke | Selected topics in column generation[END_REF] would have a smoother repartition of the reduced costs. A stake of the computational experiments is to determine whether these two generation schemes are complementary or if one dominates the other. To generate both types of columns, Algorithm 3 can be applied to generate initial columns before a smoothing local search phase to optimize the cost (16.1) with a MIP-VND similar to [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF].

Tabu Search Matheuristic intensification

This section is motivated by the following fact: many good columns can be obtained with slight modifications from a good column. At one iteration, high values of π indicate the jobs that are highly encouraged to be generated in a new column. A kernel of jobs with high π j values is likely to be combined in several interesting columns regarding the reduced costs. Using such property, it is interesting to generate aggressively columns with local search moves from interesting columns. This section introduces a Tabu Search (TS) matheuristic to generate quickly a pool of good solutions.

We denote the binaries v i,j = j u i jj indicating if technician i ∈ I 0 realizes job j. Having N feasible solutions previously calculated, we denote with ṽN i,j the value of these variables. To forbid already generated columns, we add the following "no-good-cuts" constraints similarly to [START_REF] Lazic | Variable Neighbourhood Decomposition Search for 0-1 Mixed Integer Programs[END_REF]:

∀n ∈ [[1, N ]], i,j:ṽ n i,j =1 (1 -v n i,j ) + i,j:ṽ n i,j =0 v n i,j 1 (17)
To search around the last solution ṽN i,j , allowing k modifications from the N -th solution, it can also be written as a linear constraints: Add constraint [START_REF] Rousseau | Interior point stabilization for column generation[END_REF] in MIP with columns of p Add constraint [START_REF] Pokutta | France telecom workforce scheduling problem: a challenge[END_REF] in MIP with columns of p solve MIP remove constraint [START_REF] Rousseau | Interior point stabilization for column generation[END_REF] in MIP with columns of p update p, the optimal columns in the last MIP update l adding the columns of p with a negative reduced cost in l while n < N and ReducedCost(p) < 0 return l // the list of column to add in the next RMP These constraints allow to generate a TS procedure in Algorithm 4 to generate aggressively columns for the classic CG or the last POPMUSIC CG scheme. In the case of the classic CG algorithm, the input for I 0 is a singleton in the Algorithm 4. Adding such constraints in the CG subproblems can be solved with Branch&Bound for small MIP computations. For larger MIP computations, these constraints respect the structures of the MIP matheuristics from [START_REF] Dupin | Matheuristics for the Discrete Unit Commitment Problem with Min-Stop Ramping Constraints[END_REF] to find quickly good solutions.

i,j:ṽ N i,j =1 (1 -v N i,j ) + i,j:ṽ N i,j =0 v N i,j k (18)

General algorithm

Algorithm 5 gathers POPMUSIC and TS stabilization of the CG algorithm, both techniques can be processed together. With parameter N = -1, it deactivates TS stabilization. To deactivate POPMUSIC stabilization, the partitioning procedure always split I into singletons. the ones from [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF], with a graduation of highly-constrained instances to instances more similar to VRPTW instances with no skill constraints.

Algorithm 5 provides a very significant acceleration of the CG iterations for the CG heuristics of [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF]. Solving subproblems with matheuristics provides firstly a significant acceleration. We note that the quality of MIP-based matheuristics to solve the subproblems [START_REF] Lübbecke | Selected topics in column generation[END_REF] allows to have fewer iterations with the exact computations at the end of Algorithm 5. To run repairing heuristics based on the RMP, two iterations were enough for the last phase of Algorithm 5.

First iterations of all CG schemes

Last CG iterations of all CG schemes

Comparison between GC2 and GC5 Comparison between GC2 and GC4 The acceleration of the CG algorithm concerned also the reduction of the number of iterations to converge. The aggressive generation of several columns for each subproblems with several matheuristics is a first way to stabilize the CG algorithm. To understand the specific contributions of the different CG schemes presented in this paper, we provide here comparative analyses based on several schemes solving exactly subproblems for small and mid-size instances, and a proven termination of the CG algorithm. Following CG schemes are implemented and compared:

-CG1: it is the classic CG scheme, as written in the Algorithm 1.

-CG2: The classical CG scheme is deployed with a TS intensification to solve single technician subproblems. It corresponds to the Algorithm 5 with only partitions into singletons, with parameters k = 3 and N = 5. -CG3: it is the POPMUSIC CG scheme without TS intensification, where the subproblems ( 16) are solved only using the Algorithm 3 and exact computations. -CG4: it is the POPMUSIC CG scheme without TS intensification, where the subproblems ( 16) are solved twice: Algorithm 3 gives first solutions and a second phase optimize the summed reduced cost (and thus balancing the reduced costs) with a VND similar to [START_REF] Dupin | Matheuristics for a VRPTW with competence constraints[END_REF], generating all the columns with a negative reduced cost got from these two phases. -CG5: it corresponds to the strategy CG4 with TS intensification activated with k = 3 and N = 5.

The dominant computational results are illustrated in Fig. 1. The known hierarchies among CG schemes are significant: CG1 is significantly dominated by the other schemes, POPMUSIC and TS intensification induce a very significant acceleration of the CG convergence. Comparing CG3 and CG4 allows to conclude that the optimization of the summed reduced cost in subproblems ( 16) induces a better CG convergence than the CG with the hierarchical decomposition of Algorithm 3. Combining POPMUSIC and TS intensification dominates the single stabilization strategies: POPMUSIC column generation is mostly useful for the first iterations to recombine columns in the next RMP, where it can be difficult to find a solution allocating all the jobs (and thus paying no prohibitive penalty). TS intensification is crucial when dual variables have a relative stability, which is helpful for the latest phases of the CG convergence. To compute the extended LP relaxation to optimality, TS strategies are prominent to reduce the number of CG iterations.

Conclusion and perspectives

This work outlined how matheuristics can be used to solve different subproblems from the classical CG schemes to reduce the time required for the convergence of a CG scheme. It induces a better stability of the dual variables and reduces the number of iterations required for the CG convergence. More specifically, POPMUSIC CG imposes a diversification of columns which is useful in the first stages, whereas an aggressive generation of solution with a Tabu Search matheuristic is efficient when dual are stabilized. Combining both techniques showed promising results.

The perspectives of these results are to decrease the time necessary for a CG convergence. It increases the size of instances where dual bounds are computable. It accelerated also CG matheuristics having earlier good dual variables to generate columns for primal heuristics. As a perspective, the CG strategies proposed in this work can be applied to other problems where CG induce independent and heterogeneous subproblems. The comparison/combination of the matheuristic stabilization with exact stabilization methods will also be investigated.

]]

  Working time window for technician i, to leave and come back to the depot. Time windows to begin the job j ∈ J .

Algorithm 1 :

 1 Standard column generation algorithmInput:

Algorithm 3 :

 3 Diversification of subproblem solutions Input:

Algorithm 4 :

 4 Tabu search intensification Input:-I0 a subset of technicians the current value in the RMP of dual variables (σ, π) -a set of initial columns c ∈ i∈I 0 Pi with a negative reduced cost in (15) -an integer N ∈ N, a maximal number of TS iterations an integer k ∈ N, a number of maximal modifications TSintensification(I0, k, N ) //Initialization: MIP, a MIP formulation for (15) related to technician i p = c initial columns Taboo list of columns l = {c} an integer n = 0 to denote iterations //Loop to generate columns with negative reduced costs do:

Fig. 1 .

 1 Fig. 1. Comparison of the CG convergence for an instance with 40 jobs and 15 technicians with 10 different skills

Table 1 .

 1 Table of notations, set and indexes

	I	Set of technicians.

J Set of customers, the jobs for the technicians. Ji Subset of jobs that technician i can complete with skill constraints. Dj Duration of job j. Pj Cost penalization (or outsourcing cost) if job j is not planned. di Starting and finishing depot for technician i. D(j, j )

solve RMP ( 14) with columns defined in C store dual variables σ and π and optimal cost from ( 14) Compute PI a partition of I in small subsets.

for each subset I0 ∈ PI : solve ( 16) with a matheuristic with last (σ, π) values if N > 0 then apply TSintensification(I0, k, N ) for each column c with a negative reduced cost add the column to C end for end for while : columns are added in C // Final column generation phase do :

solve RMP [START_REF] Lazic | Variable Neighbourhood Decomposition Search for 0-1 Mixed Integer Programs[END_REF] with columns defined in C store dual variables σ and π and optimal cost from ( 14) for each technician i ∈ I : solve (15) to optimality with last (σ, π) values if N > 0 then apply TSintensification({i}, k, N ) add in C the columns with a negative reduced cost end for while : columns are added in C or the specific stopping criterion is reached return the current RMP solution and the last cost of the RMP (14)

Computational results

The computational experiments were computed with a laptop Intel Core2 Duo processor, 2.80GHz, running Linux Ubuntu. The code was implemented in Python using PuLP to call MIP solvers like Cplex or open source tools with the same programming interface. The instances for the computational experiments are


Ilyes Mnassri
email: ilyes.mnassri@ensam.eu

Adil El Baroudi
email: adil.elbaroudi@ensam.eu

Elasto-Acoustic Coupling Between Two Circular Cylinders and Dense Fluid

Keywords: Elasto-acoustic interaction, three-dimensional coupled vibration, thick and thin structures, frequency analysis, finite element method

This paper describes a theoretical method for free vibration analysis of two elastic and isotropic cylinders filled with a dense fluid. The free vibration of two cylinders is studied on the basis of the linear three-dimensional elasticity theory. The compressible fluid is assumed to be nonviscous and isotropic which satisfy the acoustic wave equation. In this paper, the coupled dispersion equations of longitudinal, flexural and lobar modes are deduced and analytically solved. The finite element results computed by the Comsol Multiphysics software are compared with the present method for validation and an acceptable match between them was obtained. It is shown that the results from the proposed method are in good agreement with the numerical solutions. With this method, the effects of the cylinder parameters, such as the circumferential wave, the axial wavenumber, the thickness-to-radius and the length-to-radius on the coupled frequencies are investigated.

Introduction

The dynamic behavior of multilayered cylindrical structures in contact with dense fluid have been the focus of recent research, since they are extensively used as one of the basic components in many different engineering disciplines like the marine, mechanical, nuclear, aerospace, bioengineering and civil industries. Many of these studies are based on classical or thin-shell theories [START_REF] Leissa | Vibration of Shells[END_REF][START_REF] Zhang | Frequency analysis of submerged cylindrical shells with the wave propagation approach[END_REF]. For dynamic analysis of complex structures, thin or shell theories are most appropriate, but they lose their accuracy for thick structures, making the use of three-dimensional theory of elasticity inevitable.

The vibration analyses of multilayered cylindrical structures have been gaining an increasing attention only in the recent years. Much of the interest could be attributed to the important significance of the multilayered cylindrical structures in various practical problems contributing to noise and structural vibration. Missiles, rockets, pipes and power transmission shafts are typical cylindrical structures. When the frequency of an acting force on a cylinder coincides with it's natural frequency, the resonance phenomenon occurs and may cause failure in the hole system. Therefore, the prediction of natural frequencies is crucial and has a direct impact on the design of multilayered cylindrical structures.

The shell theory could be used for thin cylindrical structures, while thicker cylinders could only be accurately studied by using the three-dimensional elasticity theory. It should be noted that the theory of elasticity could be applied to thick or thin cylindrical structures. Gazis [1959a,b] studied the free harmonic waves in a hollow cylinder of circular cross-section, within the framework of the linear elasticity theory. He derived a characteristic equation using the Helmtholtz potentials. [START_REF] Armenakas | Propagation of harmonic waves in composite circular cylindrical shells. Part I: Theoretical investigation[END_REF] derived the frequency equation for harmonic waves, with an arbitrary number of circumferential nodes, traveling in traction-free composite, circular cylindrical structures. [START_REF] Armenakas | Free Vibrations of Circular Cylindrical Shells[END_REF] presented the transmission of elastic energy by means of elastic waves, and formulated the eigenvalue problem for stress-free cylinders for different circumferential wavenumbers. Their analysis was based on the linear three-dimensional elasticity theory. [START_REF] Xuebin | Study on free vibration analysis of circular cylindrical shells using wave propagation[END_REF] presented a wave propagation approach for free vibration analysis of cylindrical shell based on Flugge classical thin shell theory. A domain decomposition technique for solving vibration problems of uniform and stepped cylindrical shells with arbitrary boundary conditions was presented in [START_REF] Qu | Free and forced vibration analysis of uniform and stepped circular cylindrical shells using a domain decomposition method[END_REF]. [START_REF] Brischetto | An exact 3D solution for free vibrations of multilayered cross-ply composite and sandwich plates and shells[END_REF] proposed a three-dimensional free vibration analysis of multilayered structures and developed an exact solution for the differential equations of equilibrium written in general orthogonal curvilinear coordinates.

Other approaches for the vibration analysis of cylindrical structures have been proposed. [START_REF] Mcniven | Axially symmetric waves in finite, elastic rods[END_REF][START_REF] Mcniven | Axially symmetric waves in hollow, elastic rods: Part I[END_REF] used a three-mode theory to obtain a good agreement with the exact solution based on the three-dimensional theory of elasticity. El [START_REF] El Baroudi | Fluid structure interaction within three-dimensional models of an idealized arterial wall[END_REF] investigated the propagation of torsional, longitudinal, flexural and lobar waves in arterial wall layers numerically by using finite element method. [START_REF] Hamidzadeh | Vibrations of Thick Cylindrical Structures[END_REF] studied the vibrations of thick cylindrical structures. They employed the three-dimensional elastic theory and developed the governing equations of motion in terms of Helmholtz potentials. [START_REF] Gladwell | Natural frequencies of free finite length circular cylinders[END_REF] used a finite element approach to the three-dimensional vibration of finite length circular cylinders. [START_REF] Singal | A theoretical and experimental study of vibrations of thick circular cylindrical shells and rings[END_REF] developed approximate solutions for natural frequencies of thick cylinders and rings, and validated their analysis by experimental investigations. [START_REF] Ye | Three-dimensional free vibration analysis of thick cylindrical shells with general end conditions and resting on elastic foundations[END_REF] presented a unified method based on the three-dimensional theory of elasticity developed for the free vibration analysis of thick cylindrical shells with general end conditions and resting on elastic foundations. In the paper [START_REF] Alibeigloo | Three-dimensional static and free vibration analysis of carbon nano tube reinforced composite cylindrical shell using differential quadrature method[END_REF], a free vibration analysis of carbon nanotubes reinforced composite cylindrical shell was carried out using the three-dimensional theory of elasticity.

In most of the studies mentioned above, the frequency equation of the vibration problem was obtained without taking into account the fluid's effect. The originality of this paper is to investigate the effect of a dense fluid on the modal analysis of two elastic circular cylinders (Fig. 1). This paper also discusses the influence of some model parameters such as the axial wavenumber k z , the circumferential mode n, the thickness-to-radius parameter h/a and the length-to-radius parameter L/a on the natural frequencies. Indeed, the frequency equation was obtained for the free vibration model using the wave propagation method. The comparison of the results by the present method with those obtained by numerical finite element method has been carried out. It is shown that the present approach is simple, robust and gives accurate natural frequencies.

Formulation

Governing equations

As shown in Fig. 1, an isotropic circular cylinder of finite length with two layers filled with a compressible fluid in the cylindrical coordinate system (r, θ, z) is considered in the present work. The thickness of every layer is determined by radius a, a 1 and a 2 where a and a 2 are the radius of the inner and outer surfaces of the cylinder, respectively. It is assumed that the cylinder is free on the outer surface. The layers are assumed to be perfectly bonded. The system displacements and stresses are defined by the cylindrical coordinates r, θ and z. The length of the cylinder is L (0 ≤ z ≤ L). The values related to the inner and outer cylinders will be denoted by the upper indices (1) and (2), respectively. The r, θ and z axes are taken in the radial, circumferential and axial directions of the cylinder, respectively. The corresponding displacement components at any point of the cylinder in the radial, circumferential and axial directions are separately denoted by u

(i) r , u (i) θ and u (i)
z . According to the generalized linear Hooke's law, the elastic material under consideration is assumed to be linear, homogeneous, and isotropic for which the constitutive equation during deformation may be written as

σ (i) = λ i (∇ • u (i) )I + µ i ∇u (i) + ∇u (i) T (1)
where ρ i is the mass density, (λ i , µ i ) are Lamé constants, t is the time and i = 1, 2 indicates the inner (i = 1) and the external (i = 2) cylinders, respectively. For the frequency analysis considered in this paper, in the absence of body forces, the displacement field is governed by the classical Navier's equation

ρ i ∂ 2 u (i) (r, θ, z, t) ∂t 2 = (λ i + 2µ i )∇∇ • u (i) (r, θ, z, t) + µ i ∇ 2 u (i) (r, θ, z, t) (2)
Here,

u (i) = {u (i) r , u (i) θ , u (i)
z } is the vector displacement that can advantageously be expressed as a sum of the gradient of a scalar potential and the curl of a vector potential [START_REF] Morse | Methods of Theoretical Physics[END_REF]]

u (i) = ∇φ (i) + ∇ × (ψ (i) êz ) + ∇ × ∇ × (χ (i) êz ) ( 3
)
where êz is the unit vector along the axial coordinate directions and φ (i) , ψ (i) and χ (i) are the Helmholtz potentials. Substituting Eq. (3) into Eq. (2) leads to the following second-order partial differential equation with constant coefficients:

∇ 2 φ (i) = 1 c 2 Li ∂ 2 φ (i) ∂t 2 , ∇ 2 ψ (i) = 1 c 2 T i ∂ 2 ψ (i) ∂t 2 , ∇ 2 χ (i) = 1 c 2 T i ∂ 2 χ (i) ∂t 2 (4)
where c Li = (λ i + 2µ i )/ρ i and c T i = µ i /ρ i are the propagation velocities of dilatational and shear waves in the elastic solid, respectively. Furthermore, using Eq. ( 3), the relevant displacement components in cylindrical coordinates in terms of compressional and shear wave potentials may simply be written as

u (i) r = ∂φ (i) ∂r + 1 r ∂ψ (i) ∂θ + ∂ 2 χ (i) ∂r∂z , u (i) θ = 1 r ∂φ (i) ∂θ - ∂ψ (i) ∂r + 1 r ∂ 2 χ (i) ∂θ∂z (5) u (i) z = ∂φ (i) ∂z + ∂ 2 χ (i) ∂z 2 -∇ 2 χ (i)
Also, the pertinent stress components associated with the scalar potentials can be derived using the three-dimensional constitutive equation (1), which leads to the following expressions:

σ (i) rr = µ i 2 ∂ 2 φ (i) ∂r 2 + λ µ ∇ 2 φ (i) + 2 r 2 ∂ ∂θ r ∂ψ (i) ∂r -ψ (i) + 2 ∂ 3 χ (i) ∂r 2 ∂z σ (i) rθ = µ i 2 r 2 ∂ ∂θ r ∂φ (i) ∂r -φ (i) + ∇ 2 ψ (i) -2 ∂ 2 ψ (i) ∂r 2 - ∂ 2 ψ (i) ∂z 2 + 2 r 2 ∂ 2 ∂θ∂z r ∂χ (i) ∂r -χ (i) (6) σ (i) rz = µ i 2 ∂ 2 φ (i) ∂r∂z + 1 r ∂ 2 ψ (i) ∂θ∂z + ∂ ∂r 2 ∂ 2 χ (i) ∂z 2 -∇ 2 χ (i)
Now, to yield a solution for the proposed problem, the specific boundary and continuity conditions must be satisfied at the interface of two-layer isotropic cylinder:

• at the solid-fluid interface, the continuity of the stress and velocity of the solid and fluid is given by

σ (1) rr r=a = -p r=a , σ (1) rθ r=a = σ (1) rz r=a = 0, jωu (1) r r=a = v r r=a (7)
where v r is the fluid velocity components along radial direction and p is the pressure. • at the internal interface of layers (i = 1) and (i = 2), there is equality of each displacement and of stress:

σ (1) rr r=a1 = σ (2)
rr r=a1

, σ

(1) rθ r=a1

= σ

(2) rθ r=a1

, σ (1)

rz r=a1 = σ (2) rz r=a1 u (1) r r=a1 = u (2) r r=a1 , u (1) θ r=a1 = u (2) θ r=a1 , u (1) z r=a1 = u (2) z r=a1
(8)

• the outside (i = 2) surface of the two cylinders are completely stress-free

σ (2) rr r=a2 = σ (2) rθ r=a2 = σ (2) rz r=a2 = 0 (9)
We therefore examine the equations of motion for the fluid. Following the standard methods of theoretical acoustics, the field equations for an inviscid and ideal compressible medium that cannot support shear stresses may conveniently be expressed in terms of a scalar velocity potential as [START_REF] Pierce | Acoustics: An Introduction to Its Physical Principles and Applications[END_REF]

v = -∇Φ, p = -jωρ 0 Φ, ∇ 2 Φ + ω 2 c 2 0 Φ = 0 (10)
where c 0 is the sound speed in the fluid medium and ρ 0 is the ambient density.

In this way, inspection of the considered wave dispersion problem is reduced to the study of the eigenvalue problem formulated through Eqs. ( 4) and ( 10) and conditions ( 7)-( 9).

Solution method and frequency equation

Solution method

For the investigation of the dynamical problem of two-layer isotropic cylinder filled with a compressible fluid, various solution methods, including numerical (FEM and BEM for example) and analytical methods are used. In the present study, for the solution of the problem under consideration, we use both analytical and numerical methods. The resolution procedure starts with finding the exact solution of the problem and then the eigenvalue equation is obtained from the corresponding contact and boundary conditions. Finally, the solution of this eigenvalue equation is numerically computed using the Comsol Multiphysics software [Comsol Multiphysics, 2016]. In order to solve the free vibration problem of simply supported circular cylindrical isotropic structures, a time harmonic dependence exp(jωt) is assumed, where j is the imaginary unit, ω the circular frequency and t the time. The three subproblems consist of determining the scalar potentials φ (i) (r, θ, z), ψ (i) (r, θ, z) and χ (i) (r, θ, z) satisfying the Helmholtz equations (4). By going through the usual separation of variable technique, the general solution associated with an axial wavenumber k z and circumferential modal parameter n in cylindrical coordinates can be expressed as

       φ (i) ψ (i) χ (i)        (r, θ, z, t) =        φ (i) (r) sin(nθ) sin(k z z) ψ (i) (r) cos(nθ) sin(k z z) χ (i) (r) sin(nθ) cos(k z z)        exp(jωt) (11)
where φ (i) (r), ψ (i) (r) and χ (i) (r) are defined as

φ (i) (r) = A (i) 1 J n (k φ i r) + A (i) 2 Y n (k φ i r) ψ (i) (r) = A (i) 3 J n (k ψ i r) + A (i) 4 Y n (k ψ i r) χ (i) (r) = A (i) 5 J n (k ψ i r) + A (i) 6 Y n (k ψ i r)
In the above equations, k φ i is the modified wavenumber associated with the dilatational wave, k ψ i is the modified wavenumber associated with the shear wave, J n and Y n are the nth order standard Bessel function of the first and the second kind, respectively; whereas

A (i) 1 , A (i) 2 , A (i) 3 , A (i) 4 , A (5) 1 , and
A (i)
6 are unknown wave propagation coefficients whose solutions are determined below. Note that the right axial wavenumber k z is found by satisfying the required boundary conditions at the two ends of the hollow cylinder. The modified dilatational and shear wavenumbers are calculated using

k φ i = ω 2 c 2 Li -k 2 z , k ψ i = ω 2 c 2 T i -k 2 z (12)
in which k z is the wavenumber in the z-direction. At each end of the two layered cylinder, two boundary conditions must be specified. For the simply supported cylinders, the following boundary conditions are imposed:

u (i) r z=0 = u (i) θ z=0 = u (i) r z=L = u (i) θ z=L = 0
For these two boundary conditions, it can be easily shown that the axial wavenumber k z is a function of axial mode m and cylinder length L as k z = mπ/L. The term k z c Ti defined by Eq. ( 12) is a frequency parameter for torsional modes when the shear wavenumber k ψ i is taken to be equal to zero. These frequencies can be plotted as the functions of the ratio L/a for different values of the longitudinal wavenumber m. The torsional mode is characterized by a continued rotation in a single direction. This rotation direction is continued through the length of the longitudinal wave. Therefore, the torsional wave problem does not involve the Bessel equations. Direct substitution of field expansions (11) into the field equation ( 5) leads to

       u (i) r u (i) θ u (i) z        (r, θ, z, t) =        u (i) r (r) sin(nθ) sin(k z z) u (i) θ (r) cos(nθ) sin(k z z) u (i) z (r) sin(nθ) cos(k z z)        exp(jωt) ( 13
)
where

u (i) r (r), u (i) θ (r) and u (i)
z (r) are defined as

u (i) r (r) = φ (i) (r) - n r ψ (i) (r) -k z χ (i) (r) u (i) θ (r) = n r φ (i) (r) -ψ (i) (r) - nk z r χ (i) (r) u (i) z (r) = k z φ (i) (r) + k ψ 2 i χ (i) (r)
where the prime signifies the differentiation with respect to the argument r. Then, the stress σ (i) can be readily obtained using Eqs. ( 6) and ( 11)

         σ (i) rr σ (i) rθ σ (i) rz          (r, θ, z, t) =          σ (i) rr (r) sin(nθ) sin(k z z) σ (i) rθ (r) cos(nθ) sin(k z z) σ (i) rz (r) sin(nθ) cos(k z z)          exp(jωt) (14)
where

σ (i) rr (r), σ (i) rθ (r) and σ (i)
rz (r) are defined as

σ (i) rr (r) = µ i 2φ (i) (r) - λ i µ i ω 2 c 2 Li φ (i) (r) - 2n r 2 [rψ (i) (r) -ψ (i) (r)] -2k z χ (i) (r) σ (i) rθ (r) = µ i 2n r 2 [rφ (i) (r) -φ (i) (r)] -2ψ (i) (r) -k ψ 2 i ψ (i) (r) - 2nk z r 2 [rχ (i) (r) -χ (i) (r)] σ (i) rz (r) = µ i 2k z φ (i) (r) - nk z r ψ (i) (r) + (k ψ 2 i -k 2 z )χ (i) (r)
Similarly, the solution for Eq. ( 11) which belongs to the longitudinal and shear waves inside the compressible fluid is represented in cylindrical coordinates by

Φ(r, θ, z, t) = AJ n (k Φ r) sin(nθ) sin(k z z) exp(jωt), k Φ = ω 2 c 2 0 -k 2 z ( 15
)
where J n (k Φ r) is the Bessel function of the first kind and A is the unknown modal coefficient. Note that φ (i) (or ψ (i) , χ (i) , Φ) has components that are symmetric or antisymmetric in θ and z. Following standard practice, the solutions with symmetric (antisymmetric) axial φ (i) are called the antisymmetric (symmetric) axial modes, respectively, with k a z and k s z denoting the corresponding eigenvalues. Thus, we have k a z = (m -1 2 )π/L and k s z = mπ/L. However, the azimuthal modes corresponding to cos(nθ) and sin(nθ) are really the same, due to periodicity in the azimuthal direction, i.e., there is no distinction in the values of n for the two families.

Frequency equation

Now, we shall derive the frequency equation for the isotropic circular cylinder. Combining the continuity and boundary conditions ( 7)-( 9) with Eqs. ( 13)-( 15) yields for each vector mode number (n, m) to the following linear system:

NY = 0 Y = {A (1) 1 , A (1) 2 , A (1) 3 , A (1) 4 , A (1) 5 , A (1) 6 , A (2) 1 , A (2) 2 , A (2) 3 , A (2) 4 , A (2) 5 , A (2) 6 , A} T (16)
The matrix N is defined as follows:

N =                                    N (1)
11 N

(1)

12 N

(1)

13 N

(1)

14 N

(1)

15 N

(1)

16 N

(2)

11 N

(2)

12 N

(2)

13 N

(2)

14 N

(2)

15 N

(2)

16 N 1 N (1) 21 N (1) 22 N (1) 23 N (1) 24 N (1) 25 N (1) 26 N (2) 21 N (2) 22 N (2) 23 N (2) 24 N (2) 25 N (2) 26 0 N (1) 31 N (1) 32 N (1) 33 N (1) 34 N (1) 35 N (1) 36 N (2) 31 N (2) 32 N (2) 33 N (2) 34 N (2) 35 N (2) 36 0 N (1) 41 N (1) 42 N (1) 43 N (1) 44 N (1) 45 N (1) 46 N (2) 41 N (2) 42 N (2) 43 N (2) 44 N (2) 45 N (2) 46 N 4 N (1) 51 N (1) 52 N (1) 53 N (1) 54 N (1) 55 N (1) 56 N (2) 51 N (2) 52 N (2) 53 N (2) 54 N (2) 55 N (2) 56 0 N (1) 61 N (1) 62 N (1) 63 N (1) 64 N (1) 65 N (1) 66 N (2) 61 N (2) 62 N (2) 63 N (2) 64 N (2) 65 N (2) 66 0 N (1) 71 N (1) 72 N (1) 73 N (1) 74 N (1) 75 N (1) 76 N (2) 71 N (2) 72 N (2) 73 N (2) 74 N (2) 75 N (2) 76 0 N (1) 81 N (1) 82 N (1) 83 N (1) 84 N (1) 85 N (1) 86 N (2) 81 N (2) 82 N (2) 83 N (2) 84 N (2) 85 N (2) 86 0 N (1) 91 N (1) 92 N (1) 93 N (1) 94 N (1) 95 N (1) 96 N (2) 91 N (2) 92 N (2) 93 N (2) 94 N (2) 95 N (2) 96 0 N (1) 101 N (1) 102 N (1) 103 N (1) 104 N (1) 105 N (1) 106 N (2) 101 N (2) 102 N (2) 103 N (2) 104 N (2) 105 N (2) 106 0 0 0 0 0 0 0 N (2)
111 N

(2)

112 N

(2)

113 N

(2)

114 N

(2)

115 N

(2)

116 0 0 0 0 0 0 0 N (2)
121 N

(2)

122 N

(2)

123 N

(2)

124 N

(2)

125 N

(2)

126 0 0 0 0 0 0 0 N (2)
131 N

(2)

132 N

(2)

133 N

(2)

134 N

(2)

135 N

(2)

136 0                                   
and whose elements are given in the Appendix. For nontrivial solution, the determinant of the coefficient matrix in Eq. ( 16) must be zero:

det(N) = 0 (17)
For given material and geometric properties, the expansion of Eq. ( 17) will give the frequency equation for a particular value of n. For any value of k z , the frequency equation will yield an infinite number of values of ω, each corresponding to a different mode of wave propagation. Note that the solutions for n = 0 correspond to the longitudinal modes, together with the torsional modes, for n = 1 to the flexural modes and for n > 1 to the lobar modes.

Longitudinal vibration

The longitudinal mode refers to the absence of torsional displacement u θ = 0. In this case, the displacement field is independent of the θ-coordinate and is of the form (u

(i) r , 0, u (i) z
). This mode of wave propagation corresponds to n = 0 in Eq. ( 17) and results in ψ (i) = 0 [START_REF] Morse | Methods of Theoretical Physics[END_REF]. Therefore, the boundary condition Eqs. ( 7)-( 9) become

σ (1) rr r=a = -p r=a , σ (1) rz r=a = 0, jωu (1) r r=a = v r r=a σ (1) rr r=a1 = σ (2) rr r=a1 , σ (1) rz r=a1 = σ (2) rz r=a1 u (1) r r=a1 = u (2) r r=a1 , u (1) z r=a1 = u (2) z r=a1 σ (2) rr r=a2 = σ (2) rz r=a2 = 0
Then, Eq. ( 16) becomes

N L Y L = 0, Y L = {A (1) 1 , A (1) 2 , A (1) 5 , A (1) 6 , A (2) 1 , A (2) 2 , A (2) 5 , A (2)
6 , A} T where

N L =                        N (1) 11 N (1) 12 N (1) 15 N (1) 16 N (2) 11 N (2) 12 N (2) 15 N (2) 16 N 1 N (1) 31 N (1) 32 N (1) 35 N (1) 36 N (2) 31 N (2) 32 N (2) 35 N (2) 36 0 N (1) 0 N (1) 101 N (1) 102 N (1) 105 N (1) 106 N (2) 101 N (2) 102 N (2) 105 N (2) 106 0 0 0 0 0 N (2)
111 N

(2)

112 N

(2)

115 N

(2) 116

0 0 0 0 0 N (2)
131 N

(2)

132 N

(2)

135 N

(2) 136

0                       
Solving det(N L ) = 0 gives the frequency equation for longitudinal modes.

Finite element modeling

In this section, the commercially available FEA package COMSOL Multiphysics was used to develop the finite element model of an elasto-acoustic coupling between two circular cylinders and dense fluid. This modeling procedure requires two modules: one for simulating the elastic tube and the other for the dense fluid. Each module provides a wide range of equations, which was needed in specifying subdomains and boundaries. The theories and equations behind this model are based on the governing equations in Sec. 2. The cylinders were considered to be a threedimensional solid, stress-strain problem. The elastic cylinders and fluid were simulated using Lagrange-Quadratic element and Lagrange-P 2 P 1 element, respectively. The cylinders and fluid elements at the interface shared the same nodes and had extremely fine meshes to capture the details during the coupled vibrations. The applied boundary conditions were as follows. The cylinders were simply supported at two ends (u x = u y = 0). The fluid pressure nodes at the two ends were constrained (p = 0). Moreover, for simulation of the coupled system, two modules (Fluid Flow and Structural Mechanics) were used. Two elastic cylinders were used to envelope the fluid domain. On the outer cylindrical perimeter of the solid domain, a free constraint condition was used. The coupling of elastic cylinders' vibration with the dense fluid inside the cylinders was performed using both Fluid Flow and Structural Mechanic modules. For this purpose, some variables were set to make the connection between these two modules. At fluid-structure interface, kinematic and dynamic continuity have to be ensured. The complete coupled problem has to fulfil the condition that the location of the fluid-structure interface coincides for both the fields. Thus, the fluid-structure interaction boundary condition concerning the fluid is of an inhomogeneous Dirichlet type (kinematic condition) which represents the mass conservation across the interface, whereas the fluid-structure boundary condition for the solid is given by an inhomogeneous Neumann condition (dynamic condition) which describes the equivalence of fluid stresses and solid stresses in normal direction. For the modeling of the coupled system, eigenfrequency solver was selected to solve the model to give the first eight eigenvalues.

Results and Validation

The frequency equation ( 17) gives the frequency equations of the lobar, flexural, longitudinal and torsional modes in terms of its elastic constants. We consider twolayer circular cylinder filled with a compressible fluid as shown in Fig. 1 and assume that the radius of the internal circle of the inner cylinder is a and the thickness of the middle and outer cylinders is h 1 = a 1 -a and h 2 = a 2 -a 1 , respectively. The resulting frequency equation given in Eq. ( 17) is transcendental in nature with respect to the frequency ω and wavenumber k z . To obtain the roots of the frequency equation, the Mathematica software was used [Mathematica, 2015]. To validate the analytical results, the natural frequencies and mode shapes are also computed using Comsol

Multiphysics FEM Simulation Software [Comsol Multiphysics, 2008]. To simplify the presentation of results and comparison with other established results for specific cases, natural frequencies are normalized and introduced as the nondimensional frequency, which is defined by Ω = ωh/(πc T i ), where h = h 2 + h 1 . In order to compare the results, we maintain the same frequency factor as the one defined by [START_REF] Armenakas | Free Vibrations of Circular Cylindrical Shells[END_REF] and [START_REF] Hamidzadeh | Vibrations of Thick Cylindrical Structures[END_REF] in the case of uncoupled model (cylinders in the vacuum). The focus of this paper is to validate the accuracy of the method for calculating the coupled natural frequency of two circular cylinders filled with a compressible fluid.

The same geometric and material properties used for the finite circular cylinders in [START_REF] Hamidzadeh | Vibrations of Thick Cylindrical Structures[END_REF] were taken to construct this numerical example. The following values of the isotropic layer (middle layer i = 2) were assumed: ρ i = 7747.6 (kg/m 3 ) (density), ν i = 0.32 (Poisson ratio) and E i = 2.1 • 10 11 (Pa) (Young's modulus). The two circular cylinders are simply supported at both ends and subjected to a compressible fluid with sound speed c 0 = 1500 (m/s) and mass density ρ 0 = 1000 (kg/m

3 ). Uncoupled and coupled natural frequencies are obtained analytically by Mathematica [2010] and numerically by Comsol Mathematics [2008].

Table 1 shows the comparison of the first eight uncoupled frequencies of the two thick cylinders by [START_REF] Armenakas | Free Vibrations of Circular Cylindrical Shells[END_REF] and in [START_REF] Hamidzadeh | Vibrations of Thick Cylindrical Structures[END_REF], FEM, and the present method. The uncoupled frequencies predicted by the present method are identical to those calculated by FEM and by [START_REF] Armenakas | Free Vibrations of Circular Cylindrical Shells[END_REF] and in [START_REF] Hamidzadeh | Vibrations of Thick Cylindrical Structures[END_REF], indicating the validity of our model. The comparisons show that the present method is correct with good accuracy. For an uncoupled shell's example, exact and approximate results for a cylindrical shell, simply supported at both ends, have been reported by [START_REF] Oliazadeh | Exact analysis of resonance frequency and mode shapes of isotropic and laminated composite cylindrical shells; Part I: Analytical studies[END_REF]. This shell has the following properties: E = 68.2 (GPa), ν = 1/3, ρ = 2700 (kg/m 3 ), h/a = 1/20 and L/a = 3. Table 2 shows a comparison of exact results obtained by using the present study with those of [START_REF] Oliazadeh | Exact analysis of resonance frequency and mode shapes of isotropic and laminated composite cylindrical shells; Part I: Analytical studies[END_REF] and calculated using FEM (Comsol Multiphysics). It can be seen generally that exact results obtained using the present study agree with the exact results using different shell theories (see [START_REF] Oliazadeh | Exact analysis of resonance frequency and mode shapes of isotropic and laminated composite cylindrical shells; Part I: Analytical studies[END_REF]). In the case of two cylinders with a compressible fluid (Table 3), a very good agreement was observed between the results of the present method and those of FEM throughout the very low relative difference ((FEM-Present)/Present) which is actually ≤ 0.01%. Figures 2-5 illustrate the first three coupled modal shapes. The modal shape can be regarded as the mode (n, m), where n is the modal number in the circumferential direction and m is the modal number in the axial direction. The modal shapes in Figs. 2 and4 are not in order with the parameters n and m in the case of thin cylinders (short or long). The frequency of mode (5, 1) is lower than that of mode (4,1), for example. This feature of thin cylinders' vibration is different from that of the beam vibration in which the order increases with the modal parameter. Therefore, in the vibration of thin cylinders, one should be careful to find the right mode of the vibration. On the other hand, in Figs. 3 and5, the modal shapes are in order with the parameters n and m in the case of thick cylinders (short or long). The frequency of mode (2, 1) is lower than that of mode (1,1), for example.

For the free vibration analysis of short (or long) and thick (or thin) two elastic and isotropic cylinders subjected to compressible fluid, the effects of cylinder parameters on the coupled natural frequencies are presented in the present method. In these studies, investigations were carried out to study the effects of circumferential mode n, axial wavenumber k z , thickness-to-radius parameter h/a and length-toradius parameter L/a on the coupled frequencies. We begin by examining how the coupled frequencies for the given h/a and L/a ratios vary with the axial mode m and circumferential mode n. Figure 6 exhibits that for h/a = 0.02 (case of thin structures), with L/a = 1.33 and L/a = 6, the frequencies first decrease and then increase as the circumferential mode n increases. This behavior is also found for short thick structures (Fig. 8). Figure 7 indicates that both the m frequency curves converge when n is large which means that for large circumferential mode n, the frequencies are only determined by the n and regardless of the axial mode m. This behavior is also found for long thick structures (Fig. 9). For a given n, the frequency increases as m increases. Therefore, the fundamental frequency is always with the m = 1 curve. In the following, we study the frequency curves of m = 1 (k z L = π) with the change of the other cylinder parameters. For better understanding how the dense fluid (added mass) affects the frequencies, Figs. 10-13 have been plotted. Theses figures display the coupled and uncoupled frequencies for h/a = 0.02 and h/a = 0.7, with L/a = 1.33 and L/a = 6, varying with the circumferential mode n. As n increases, the difference between the coupled and uncoupled frequencies reduces at 0 < n < 3, but increases after the fundamental frequency at n = 3 for long thin structures (Fig. 11). For the thick structures, as n increases, the difference between the coupled and uncoupled frequencies reduces at 0 < n < 6, but intensifies after the fundamental frequency at n = 6 for long thick structures (Fig. 13). However, for the short thick structures (Fig. 12), it is observed that as n increases, the difference between the coupled and uncoupled frequencies reduces at 0 < n < 4, but amplifies after the fundamental frequency at n = 4.

A final inspection into cylinder vibrations is given by investigating how the coupled and uncoupled frequencies vary with L/a ratios. Figures 10-12 show the uncoupled and coupled frequencies, with h/a = 0.02 and h/a = 0.7, varying with the L/a parameter. The frequency curves of n = 0 (longitudinal modes), n = 1 (flexural modes) and n = 2 (lobar modes) are presented. Due to the coupling effects, all the frequencies decrease as the L/a ratios raise for the curve of n = 0 (Figs. 14 and15). It is noted that the frequency spectrum for the vibration modes n = 1 and n = 2 for thin structures has similar behavior to that of n = 0 (Figs. 16 and18). It can be observed in Figs. 16-19, respectively, that the coupled and uncoupled frequencies, with h/a = 0.02 and h/a = 0.7, are varying with the parameter L/a for n = 1 and n = 2. For all the above figures, both the coupled and uncoupled frequencies decrease as the L/a ratios increase, but at different rates, the difference between the coupled and uncoupled frequencies reduces as the L/a ratios increase. This shows the different coupling effects on the different mode vibrations n with the L/a ratios. As expected, these figures also register that in the case of thick and long structures, the added mass has no significant effect (Figs. 17 and 19).

Conclusion

The characteristics of vibration analysis of two circular cylinders subjected to compressible fluid based on the linear three-dimensional elasticity theory are investigated. The spatial displacement components in the three coordinate directions are described by Helmholtz potential functions. The governing equations are solved to obtain the exact frequency equation. Frequency equations are presented for uncoupled and coupled problems. The nondimensional frequency parameters obtained by the present approach agree very well with those available in the literature in the case of uncoupled problem. However, in the case of coupled problem, the results obtained in this paper can serve as benchmark solutions and as a reference to assess the accuracy of approximate methods. High accuracy have been observed in the analysis. For demonstration, the natural frequencies of different aspect ratios have been studied in detail, and the first eight frequency parameters for various mode categories are given. The validity of the present solution is solidly confirmed numerically. Additionally, the studying of the thermal coupling effect should be explored in the future work.

  Fig. 1. The geometry of the cylinder with two layers filled with compressible fluid.

Fig. 2 .

 2 Fig. 2. First three coupled Ωn,m and associated mode shapes (h/a = 0.02 and L/a = 1.33).

Fig. 6 .

 6 Fig. 6. Variation of coupled frequency parameter Ω versus n for h/a = 0.02 and L/a = 1.33.

Fig. 8 .

 8 Fig. 8. Variation of coupled frequency parameter Ω versus n for h/a = 0.7 and L/a = 1.33.

Fig. 9 .

 9 Fig. 9. Variation of coupled frequency parameter Ω versus n for h/a = 0.7 and L/a = 6.

Fig. 11 .

 11 Fig. 11. Variation of frequency parameter Ω versus n for h/a = 0.02 and L/a = 6.

Fig. 13 .

 13 Fig. 13. Variation of frequency parameter Ω versus n for h/a = 0.7 and L/a = 6.

Fig. 14 .

 14 Fig. 14. Frequency curves for longitudinal modes (n = 0), as Ω versus L/a with h/a = 0.02.

Fig. 15 .

 15 Fig. 15. Frequency curves for longitudinal modes (n = 0), as Ω versus L/a with h/a = 0.7.

Fig. 17 .

 17 Fig. 17. Frequency curves for flexural modes (n = 1), as Ω versus L/a with h/a = 0.7.

Fig. 18 .

 18 Fig. 18. Frequency curves for lobar modes (n = 2), as Ω versus L/a with h/a = 0.02.

Table 1 .

 1 Comparison of the first eight Ωn,m = hωn,m/(πc T i ) in the case of uncoupled vibration for h/a = 0.66 and L/a = 1.33.

	(n, m) Present	FEM	Hamidzadeh et al.	Armenakas et al.	Mode shape
	(1, 1)	0.361	0.361	0.359	0.359	Flexural
	(0, 1)	0.371	0.371	0.369	0.369	Longitudinal
	(2, 1)	0.375	0.375	0.373	0.373	Lobar
	(3, 1)	0.444	0.444	0.441	0.441	Lobar
	(0, 1)	0.499	0.499			Torsional
	(1, 1)	0.536	0.536	0.536	0.536	Flexural
	(4, 1)	0.553	0.553	0.550	0.550	Lobar
	(2, 1)	0.614	0.614	0.614	0.614	Lobar

Table 3 .

 3 Comparison of the first eight Ωn,m = hωn,m/(πc T i ) in the case of coupled vibration for h/a = 0.66 and L/a = 1.33.

	(n, m)	Present	FEM	Mode shape	Order
	(0, 1)	0.2225	0.2225	Longitudinal	First
	(1, 1)	0.2776	0.2776	Flexural	First
	(2, 1)	0.3422	0.3422	Lobar	First
	(0, 1)	0.3655	0.3655	Longitudinal	Second
	(1, 1)	0.3678	0.3678	Flexural	Second
	(2, 1)	0.4016	0.4016	Lobar	Second
	(3, 1)	0.4194	0.4194	Lobar	First
	(0, 1)	0.4407	0.4407	Longitudinal	Third

Acknowledgments

We would like to sincerely thank Dr. F. Razafimahery for his useful comments, suggestions, agreement and helpful hints to use FEM Comsol Multiphysics.

N (2) 86

Appendix A A.1. For i = 1, 2, the following formulas may expediently be utilized to obtain the elements of the matrix given in Eq. (17) :

Coefficients of the corresponding displacement components in the radial, circumferential and axial directions (Eq. ( 13))

A.3. Coefficients of the pertinent stress components (Eq. ( 14))

A.4. Elements of the matrix given in Eq. ( 17)

4 (a 1 )

4 (a 1 )


Adil El Baroudi
email: adil.elbaroudi@ensam.eu

Fulgence

Fulgence Razafimahery

Transverse vibration analysis of Euler-Bernoulli beam carrying point masse submerged in fluid media

Keywords: Euler-Bernoulli beam, Fluid-structure interaction, Finite element method, Frequency equation, Inertial coupling

may come

Transverse vibration analysis of Euler-Bernoulli beam carrying point masse submerged in fluid media

Introduction

Systems of a beam or rod carrying masses are frequently used as design models in engineering. The operation of the machine may introduce severe dynamic stresses on the beam. It is important, then, to know the natural frequencies of the coupled beam-mass and beam-mass-fluid system, in order to obtain a proper design of the structural components. Numerous papers have been published on the vibration analysis of beams carrying concentrated masses at arbitrary positions [START_REF] Low | Natural frequencies of a beam-mass system in transverse vibration: Rayleigh estimation versus eigenanalysis solutions[END_REF][START_REF] Stanisic | On the response of beams to an arbitrary number of concentrated moving masses[END_REF][START_REF] Laura | Optimization of eigenvalues in the case of vibrating beams with point masses[END_REF][START_REF] De Rosa | On the dynamic behaviour of slender beams with elastic ends carrying a concentrated mass[END_REF][START_REF] Low | On the eigenfrequencies for mass loaded beams under classical boundary conditions[END_REF]. The problem of a vibrating simply supported beam carrying a concentrated mass at its center was solved analytically by [START_REF] Chen | On the Vibration of Beams or Rods Carrying a Concentrated Mass[END_REF]. The Laplace transformation technique to formulate the frequency equation for beams carrying intermediate concentrated masses have been used in [START_REF] Liu | Free vibration of beams with elastically restrained edges and intermediate concentrated masses[END_REF]. In [START_REF] Laura | A note on transverse vibrations of continuous beams subject to an axial force and carrying concentrated masses[END_REF] Rayleigh-Ritz method was used to study continuous beams subjected to axial forces and carrying concentrated masses. Laura [START_REF] Laura | A note on the vibration of a clamped-free beam with a mass at the free end[END_REF] studied the cantilever beam carrying a lumped mass at the top in order to obtain analytical solution by introducing the mass in the boundary conditions. Other studies of the influence of these factors on natural frequencies of beams and rods are given in [START_REF] Naguleswaran | Transverse vibrations of an Euler-Bernoulli uniform beam carrying several particles[END_REF][START_REF] Rossit | Transverse vibrations of a cantilever beam with a spring mass system attached on the free end[END_REF][START_REF] Su | Exact natural frequencies of structures consisting of two part beam-mass systems[END_REF][START_REF] De Rosa | The use of Mathematica in the dynamics analysis of a beam with a concentrated mass and dashpot[END_REF][START_REF] Addessi | Free in-plane vibrations of highly buckled beams carrying a lumped mass[END_REF][START_REF] Sundara Raja Iyengar | Free Vibration of Rectangular Beams of Arbitrary Depth[END_REF]. In most of the studies mentioned above the obtaining the natural frequencies of beams carrying point masses is performed in vacuum. This paper presents an aspect of fluid-structure interaction in order to investigate the influence of the added mass due to fluid on the natural frequencies of a beams carrying point masses submerged in fluid (water). Fluid-structure interaction problems since long have attracted the attention of engineers and applied mathematics. The most important applications of this theory, probably, structural acoustics [START_REF] Ohayon | Structural Acoustics and Vibration[END_REF], vibrations of fluid-conveying pipes [START_REF] Paidoussis | Fluid-Structure Interactions Slender Structures and Axial Flow[END_REF] and biomechanics. In fact many engineering fields have related interest, including naval architecture, offshore structures, hydrodynamics, dam-reservoir systems under seismic, noise control and vibration isolation, as well as nuclear reactor plants. As these problems are rather complicated, some simplifications are typically adopted to facilitate their solving. In particular, it's quite typical to ignore viscosity effects (especially in structural acoustics) or to use local theories of interaction, such as, the one referred to as thin layer or plane wave approximation. Therefore vibrations of fluid-coupled beams have received a great attention due to their importance in various engineering applications. As a consequence of the fluid-structure interaction, there are coupled dynamical equations to be solved simultaneously. A number of papers has been published in recent years investigating the added mass effects of the interacting fluids in structures. Xing [START_REF] Xing | Natural vibration of a beam-water interaction system[END_REF] examined the dynamical behaviour of a flexible beam-water interaction system and showed that in the water domain, the natural frequencies of the coupled dynamic system are lower than those of the flexible dry beam, indicating that the influence of water on the beam has the effect of an additional mass. The present manuscript deals with the coupling effect of liquid on the free vibration characteristics of a Euler-Bernoulli beams carrying point masse submerged in water contained in a rigid rectangular container. The originality of this study is to investigate the effects of both pointed mass and fluid interaction on the dynamics of beams under the fixed-fixed boundary condition at the beam edge. The governing equations describing the behavior of the system are analyzed using the separation of variables method and their solutions presented. The eigenvalue equation of the natural vibration of the beam-water system is derived and exact solutions are obtained. Natural frequencies and modes shapes of the beam-masse-water are investigated. The theoretical results from the present formulation are compared to those from a finite element analysis using a commercial finite element analysis. Calculations show that, in the water domain, the concentrated mass become less influential on the natural vibration behavior of the coupled system. Finally, one investigates how the free surface affects the frequencies. Finding a general analytical solution for Euler-Bernoulli beams carrying point masse submerged in fluid remains a challenge.

Governing equations of fluid domain

First of all, the fluid motion must be defined in order to find the coupled natural frequencies and wet mode shapes of the beams carrying point masses in contact with the fluid. To account for the effect of the added mass in the beam response, the acoustic mode shapes and natural frequencies are first computed for the rectangular domain fluid assuming rigid and soft boundary conditions. The fluid is assumed non-viscous and isotropic which satisfies the acoustic wave equation. Since no source is considered in the acoustic domain, the Helmholtz equation, which describes a harmonic wave equation propagating in medium while neglecting dissipation, is represented as

∇ 2 p + ω 2 c 2 p = 0 (1)
where

∇ 2 = ∂ 2 ∂x 2 + ∂ 2
∂y 2 is the Laplacian operator in a cartesian coordinate system. p is the acoustic pressure in a fluid medium, ω is the angular frequency and speed of sound c. The Eq. ( 1) indicates the distribution of the acoustic pressure p in space and can be solved if the boundary condition is knew well. The natural frequencies of the acoustical system are obtained by assuming that the boundaries of the enclosure are hard except the boundary Γ p (zero free surface wave disturbance), hence the pressure gradients on the boundaries Γ 1 , Γ 2 and Γ 3 are set to zero :

∂p ∂x (0, y) = 0 , ∂p ∂x (L, y) = 0 , ∂p ∂y (x, 0) = 0 , p(x, H) = 0 (2)
The solution of (1) with boundary conditions introduced in ( 2) is given as p(x, y) = P cos (α n x) cos (α q y) , n, q = 0, 1, 2, ...

where α n = nπ/L, α q = (2q + 1) π/2H and P is unknown coefficient which will be determined later by imposing the appropriate boundary conditions. 1) and ( 2)) carrying a point masse (denoted m) submerged in a liquid-filled rectangular geometry.

Theoretical modal analysis of beams

Modal analysis of an elastic submerged structures is needed in every modern construction and have wide engineering application especially in ocean engineering. In this study, modal analysis is important to predict the dynamic behavior of the submerged beams. It is well known that the natural frequencies of the submerged elastic structures are different from those in vacuum. The effect of fluid forces on the submerged beams is represented as added mass, which decreases the natural frequencies of the submerged beams from those which would be measured in the vacuum. This decrease in the natural frequencies is caused by the increase of the kinetic energy of the fluid-beams system without a corresponding increase in strain energy.

The model allows to analyze the influence of the added masse on the dynamic behavior of beams with two classic boundary conditions. The Euler-Bernoulli beams carrying point masse is submerged inside rectangular fluid domain where the lower liquid region is represented by h. The interaction between the fluid and the Euler-Bernoulli beams carrying point mass is taken into account to calculate the natural frequencies and modes shapes of the coupled system. The dynamics of each beam portion are treated separately. It is assumed that the beam has aligned neutral axis. The equation of motion for the transverse deflection of a uniform elastic beam ignoring shear deformation and rotary inertia effects can be written in the form

EI ∂ 4 v i (x, t) ∂x 4 + ρS ∂ 2 v i (x, t) ∂t 2 = p(x, h, t) , i = 1, 2 (4)
where v i (x, t) is the lateral deflection at distance x along the length of the beam, EI, ρ and S are the flexural rigidity, the mass per unit volume and the cross-sectional area of the beams. Assuming time harmonic motion at angular frequency of the form exp (jωt), but neglecting it in the subsequent expressions for clarity, the general solutions of the ordinary differential Eq. ( 4) for the beams system, as shown in Figure 1, can be written in different segments in terms of evanescent and propagating waves as

v i (x) = A i φ(x) + B i ψ(x) + C i χ(x) + D i ϕ(x) + P cos (α m h) EI (α 4 n -k 4 ) cos (α n x) (5)
where the trigonometric functions

φ(x) = sin (kx) , ψ(x) = cos (kx)
represent the propagating waves, and the hyperbolic functions

χ(x) = sinh (kx) , ϕ(x) = cosh (kx)
represent the evanescent waves. A i is a wave propagating to the left, B i is a wave propagating to the right, C i is an evanescent wave decaying to the left and D i is an evanescent wave decaying to the right; k is the flexural wavenumber and is given by

k = √ ω ρS EI
The nine constants P, A i , B i , C i and D i (i = 1, 2) can be found by imposing the following boundary conditions. Eq. ( 5) is the general solution for the vibration modes of beams submerged in fully liquid. In the case of a coupled system, the effect of sound field on the flexible beams carrying concentrated masse must be considered. On the fluid-beams interface, the normal acceleration must be continuous. Therefore, the pressure in the fluid and the deflection of a beams v i satisfy the relation :

∂p(x, h) ∂y = ρ f ω 2 v i (x) (6)
where ρ f is the density of the fluid. Inserting Eq. (3) into Eq. ( 6) and substituting the outcome into Eq. ( 5) and after some manipulations, leading to the following equation :

A i φ(x) + B i ψ(x) + C i χ(x) + D i ϕ(x) = -P α m sin (α m h) ρ f ω 2 + cos (α m h) EI (α 4 n -k 4 ) cos (α n x) (7)
Making use of the orthogonality of trigonometric functions, both sides of the Eq. ( 7) are multiplied by cos(α j x) and then integrated over 0 < x < L to yield the following equation :

A i U + B i X + C i Y + D i Z = - P L 2 α m sin (α m h) ρ f ω 2 + cos (α m h) EI (α 4 n -k 4 ) (8)
where the coefficients U , X, Y and Z are defined as follows :

U = L 0 φ(x) cos(α n x)dx , X = L 0 ψ(x) cos(α n x)dx Y = L 0 χ(x) cos(α n x)dx , Z = L 0 ϕ(x) cos(α n x)dx
Now the expression of the lateral movement of the beams, fully in contact with the fluid v i (x) can be formulated using Eq. ( 5) and taking into account Eq. ( 8) as :

v i (x) = A i [φ(x) -βU cos (α n x)] + B i [ψ(x) -βX cos (α n x)] + C i [χ(x) -βY cos (α n x)] + D i [ϕ(x) -βZ cos (α n x)] (9)
where

β = 2 L EI α 4 n -k 4 α m tan (α m h) ρ f ω 2 + 1
From equation ( 9), it can be seen that, when β → 0, the coupled system beams-mass-water is reduced to the case of beams motion in vacuo. Now, to derive the frequency equation of beams submerged in fully fluid, one assumes that on the interface between each portions of the beam at the position l 1 , the deflection, the rotation angle, the internal shear force and bending moment of the beam must be continuous. This is satisfied when

v 1 (l 1 ) = v 2 (l 1 ) v 1 (l 1 ) = v 2 (l 1 ) v 1 (l 1 ) = v 2 (l 1 ) EIv 1 (l 1 ) + mω 2 v 1 (l 1 ) = EIv 2 (l 1 ) (10)
where primes denote differentiation with respect to the position variable x. To complete the formulation of the boundary-value problem, the four boundary conditions (left and right beam ends) for the beam ends considered in this work are specified as follows :

v 1 (0) = v 1 (0) = v 2 (l) = v 2 (l) = 0 (11)
The boundary conditions (Eqs. 10 and 11) can be evaluated, giving rise to a linear homogeneous system of eight equations :

Nx = 0 , x = [A 1 , B 1 , C 1 , D 1 , A 2 , B 2 , C 2 , D 2 ] T (12)
This system can have non trivial solutions only if the determinant of the matrix N is zero, leading to the frequency equation that takes the following form N = [N 11 , N 21 , N 31 , N 41 , N 51 , N 61 , N 71 , N 81 ] T where the matrix N is written down explicitly in the Appendix.

Results and discussion

In this paper, the results obtained by the present analytical solution for the beams carrying point masse coupled with fluid are compared with those acquired by the finite element method (FEM) using Comsol Multiphysics FEM Simulation Software to show the applicability, reliability and effectiveness of the presented formulation. The mechanical and geometrical properties of the beams used in numerical computation are : Young's modulus E = 13 • 10 In Tables 1 and2, natural frequencies of beams in (Hz) unit carrying or not point masse obtained by the present theory are presented and compared with numerical data. As can be seen, there is a very good agreement between the present results and those of FEM and the relative difference ((FEM-Present)/Present) is 0.1%. This shows that the algorithm implemented in Comsol Multiphysics [START_REF]Comsol Multiphysics 3.5a, User's Guide and Reference Guide[END_REF] software for numerical computation is highly reliable and accurate. Is the attention to use the numerical formulation in future for more general geometries and other kind of fluid-structure interactions such as those arising in aeroelasticity for example. This present work, the effects of concentrated mass and add mass on the coupled natural frequencies of beams are presented with the present method. First, one investigates how the dense fluid (added mass) affects the natural frequencies. Due to the practical applications, the natural frequencies of beams coupled with fluid are listed in Table 2. It is evident from Table 2 that when the beams is submerged fluid, the natural frequencies of the beams takes lower values. Also, tables 1 and 2 show that when the concentrated mass is placed at x = l/2 (symmetric problem), only the symmetric modes (odd modes) are affected.

Secondly, in order to see the influence of concentrated mass m on the vibration characteristics of the beams, the first six natural frequencies of beams are listed in Table 1. It is seen that the concentrated mass decreases the natural frequencies in vacuo (higher 20%). In water, it is interesting to note that the concentrated mass become less influential on the natural frequencies of the beams (Table 2). Figures 234567show also that the concentrated mass modifies the modal shapes in vacuum. The presence of the fluid has low influence on the modal shapes because the added mass is diagonal (Figures 8910111213). The coupling is light in such a case. Thirdly, one investigates how the free surface affects the frequencies. The above mentioned discussions are obtained in the case of zero free surface wave disturbance (p = 0 at Γ p in Figure 1). In order to investigate the influence of the free surface on the natural frequencies of beams, the free surface wave disturbance Γ p is governed by the equation

∂p ∂y = - 1 g ∂ 2 p ∂t 2 (13)
and g = 9.81 [m • s -2 ] is the acceleration due to gravity. Figures 141516171819shows the natural frequencies and mode shapes of the submerged beams without point mass, varying with the free surface. The free surface decreases significantly the natural frequency. This behavior is also found for submerged beams with point mass.

The study also shows that the concentrated mass do not affect the coupled natural frequencies in the case of free surface wave disturbance. Comparing different boundary conditions, it is observed that added mass factor has a greater effect with free surface on the coupled natural frequencies. Therefore, the effect of free surface is very significant. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

Conclusion

This paper presents an analytical solution for the dynamic system of Euler-Bernoulli beams carrying point mass submerged in fluid media. Separation variable technique and orthogonality of trigonometric functions are used

Figure 1 :

 1 Figure 1: Sketch of beam composed of 2 uniform beam segments (denoted by (1) and (2)) carrying a point masse (denoted m) submerged in a liquid-filled rectangular geometry.

8 [

 8 Pa], material density ρ = 2000 [kg • m -3 ], Poisson's ratio ν = 0.3, length l = 0.6 [m], width a = 0.006 [m] and thickness e = 0.006 [m]. The fluid density is taken to be ρ f = 1000 [kg • m -3 ] and the speed of sound in a fluid c = 1500 [m • s -1 ]. The modeling was set up by placing the beams in a rectangular reservoir with dimensions L = 3 [m] and H = 1 [m] and lower fluid region h = 0.5 [m]. The mass of attached masse m = 0.04 [kg].

  (a) f = 13.81 Hz (b) f = 8.57 Hz (c) f = 7.50 Hz

Figure 2 :

 2 Figure 2: The first natural frequency and mode shape of beams in "vacuo": the colors pertain to the displacement field of beams. (a) without point mass. (b) with point mass at l/3. (c) with point mass at l/2.

  (a) f = 38.07 Hz (b) f = 29.10 Hz (c) f = 38.07 Hz

Figure 3 :

 3 Figure 3: The second natural frequency and mode shape of beams in "vacuo": the colors pertain to the displacement field of beams. (a) without point mass. (b) with point mass at l/3. (c) with point mass at l/2.

  (a) f = 74.64 Hz (b) f = 73.37 Hz (c) f = 59.36 Hz

Figure 4 :

 4 Figure 4: The third natural frequency and mode shape of beams in "vacuo": the colors pertain to the displacement field of beams. (a) without point mass. (b) with point mass at l/3. (c) with point mass at l/2.

Figure 5 :

 5 Figure 5: The fourth natural frequency and mode shape of beams in "vacuo": the colors pertain to the displacement field of beams. (a) without point mass. (b) with point mass at l/3. (c) with point mass at l/2.

Figure 6 :

 6 Figure 6: The fifth natural frequency and mode shape of beams in "vacuo": the colors pertain to the displacement field of beams. (a) without point mass. (b) with point mass at l/3. (c) with point mass at l/2.

  (a) f = 257.43 Hz (b) f = 255.06 Hz (c) f = 257.43 Hz

Figure 7 :

 7 Figure 7: The sixth natural frequency and mode shape of beams in "vacuo": the colors pertain to the displacement field of beams. (a) without point mass. (b) with point mass at l/3. (c) with point mass at l/2.

Figure 8 :

 8 Figure 8: The first natural frequency and mode shape of beams in water in the case of zero free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

  (a) f = 8.34 Hz

Figure 9 :

 9 Figure 9: The second natural frequency and mode shape of beams in water in the case of zero free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

Figure 10 :

 10 Figure10: The third natural frequency and mode shape of beams in water in the case of zero free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

  (a) f = 38.49 Hz

Figure 11 :

 11 Figure 11: The fourth natural frequency and mode shape of beams in water in the case of zero free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

  (a) f = 63.84 Hz

Figure 12 :

 12 Figure12: The fifth natural frequency and mode shape of beams in water in the case of zero free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

  (a) f = 97.68 Hz

Figure 13 :

 13 Figure13: The sixth natural frequency and mode shape of beams in water in the case of zero free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

Figure 14 :

 14 Figure 14: The first coupled natural frequency and mode shape of beams in water in the case of free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

Figure 15 :

 15 Figure 15:The second coupled natural frequency and mode shape of beams in water in the case of free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

Figure 16 :

 16 Figure16: The third coupled natural frequency and mode shape of beams in water in the case of free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

  (a) f = 0.032 Hz

Figure 17 :

 17 Figure 17: The fourth fourth coupled natural frequency and mode shape of beams in water in the case of free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

  (a) f = 0.036 Hz

Figure 18 :

 18 Figure18: The fifth fourth coupled natural frequency and mode shape of beams in water in the case of free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

Figure 19 :

 19 Figure 19: The sixth coupled natural frequency and mode shape of beams in water in the case of free surface wave disturbance. The colors pertain to the displacement field of beams and to the absolute acoustic pressure.

Table 1 :

 1 The first six natural frequencies in "vacuo" of the beams with point mass neglected and considered.

	Modes	Neglected		Considered			Diference
				at l/3		at l/2
		Present FEM	Present FEM	Present FEM	(%)
	1	13.81	13.81	8.57	8.57	7.50	7.50	0
	2	38.07	38.07	29.10	29.10	38.07	38.07	0
	3	74.64	74.64	73.37	73.37	59.36	59.36	0
	4	123.38 123.38	110.07 110.07	123.38 123.38	0
	5	184.31 184.31	162.20 162.20	156.94 156.94	0
	6	257.43 257.43	255.06 255.06	257.43 257.43	0

Table 2 :

 2 The first six natural frequencies in "water" of the beams with point mass neglected and considered in the case of zero free surface wave disturbance.

	Modes	Neglected		Considered			Diference
				at l 1 = l/3	at l 1 = l/2
		Present FEM	Present FEM	Present FEM	(%)
	1	1.76	1.76	1.74	1.74	1.73	1.73	0
	2	8.34	8.34	8.00	8.00	8.34	8.34	0
	3	19.82	19.82	19.82	19.82	17.99	17.99	0
	4	38.49	38.49	35.71	35.71	38.49	38.49	0
	5	63.84	63.84	58.18	58.18	58.07	58.07	0
	6	97.68	97.68	97.34	97.34	97.68	97.68	0

Appendix B. Matrix elements

The matrix N in Eq. ( 12) is defined as follows :

N 11 = {a 11 , a 12 , a 13 , a 14 , 0, 0, 0, 0.} N 21 = {a 21 , a 22 , a 23 , a 24 , 0, 0, 0, 0.} N 31 = {a 31 , a 32 , a 33 , a 34 , a 35 , a 36 , a 37 , a 38 .}