Mamadou Diaga Seck

Emmanuel Keita

Pamela Françoise Faure

Pierrick Cavalié

Maxime Van Landeghem

Stéphane Rodts

Philippe Coussot

Paméla Faure

Pierrick Cavalie

Maxime Van Landeghem

Subflorescence and plaster drying dynamics

Keywords: Plaster, drying, subflorescence, precipitation, dewetting

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

To obtain plasterboards so useful for modern constructions, there are two important phenomena that must take place during the manufacturing process: hydration, through which plaster paste passes from a pasty state to a solid state, and drying, through which the excess of the mixing water added for workability purposes is removed. The hydration reaction can be described as the successive dissolution of hemihydrates (CaSO 4 , ½ H 2 O, initial water-to-powder ratio of 0.8) and precipitation of dihydrates (CaSO 4 , 2 H 2 O) in water [START_REF] Neville | Adsorption and Reaction I[END_REF][START_REF] Ridge | Calcium Sulphate Hemihydrate and Its Hydration[END_REF][START_REF] Hansen | Setting and hardening of gypsum plaster[END_REF][START_REF] Vignes | Une vie de plâtre[END_REF]. Detailed information on these processes may be found in [START_REF] Singh | Calcium sulphate hemihydrate hydration leading to gypsum crystallization[END_REF][START_REF] Lemarchand | Plaster Hydration at Different Plasterto-Water Ratios: Acoustic Emission and 3-Dimensional Submicrometric Simulations[END_REF]. After this chemical reaction which lasts about 40 min without accelerators [START_REF] Jaffel | Probing Microstructure Evolution during the Hardening of Gypsum by Proton NMR Relaxometry[END_REF], a significant excess of water must be removed from the material, an operation which requires a lot of energy. Understanding the drying processes in plaster pastes should contribute to develop more efficient drying techniques and procedures, and to consume less energy. Here we focus on one of the technics used in manufacture to remove the residual water: by a convective air flow [START_REF]Plaster casts: making, collecting and displaying from classical antiquity to the present[END_REF].

Existing knowledge on drying characteristics of plasters essentially concerns materials already in solid state which are imbibed then dried and tends to suggest [START_REF] Petkovic | Salt transport in plaster/substrate layers[END_REF] that as long as pure water is used such materials exhibit drying properties close to usual simple homogeneous porous media. A specific aspect has been studied in more details in recent years, namely the transport of salts during drying of solid plaster samples. It has been shown that salts are transported towards the free surface and then can either crystallize inside the plaster or move outside, towards the substrate [START_REF] Petkovic | Salt transport in plaster/substrate layers[END_REF]. It was also suggested that this salt crystallization leads to the formation of a structure inside the plaster, whose characteristics depend on the drying rate and in turn impact the drying rate [START_REF] Vosten | The drying of crust-forming materials as exemplified by the drying of gypsum wallboard[END_REF].

The drying characteristics of plaster pastes just after their preparation has not been studied in depth so far whereas it is a critical step of the process of plasterboard production. Previous studies concerning plaster properties during its preparation focused on the structure evolution in time, in particular with the help of NMR, showing that the pore distribution evolves as a result of hydration [START_REF] Song | Magnetic resonance studies of hydration kinetics and microstructural evolution in plaster pastes[END_REF][START_REF] Song | Simultaneous monitoring of hydration kinetics, microstructural evolution, and surface interactions in hydrating gypsum plaster in the presence of additives[END_REF]. A recent work finally showed that the drying rate of plaster pastes is significantly lower than that expected for a pure liquid evaporating from a simple homogeneous porous medium [START_REF] Seck | The mechanisms of plaster drying[END_REF]. This effect was shown to be enhanced by the air flow velocity and the initial solid/water ratio. From further tests under various conditions it was suggested that this effect is due to the crystallization of gypsum ions below the sample free surface, which creates a dry region and decreases the drying rate by increasing the length of the path the vapor has to follow before reaching the free surface. The specificity of plaster drying is indeed that at the end of the initial setting process (due to precipitation of dihydrates) the interstitial solution is saturated with ions (Ca 2+ , SO [START_REF] Vignes | Une vie de plâtre[END_REF] 2-), which can then precipitate as gypsum crystals and accumulate in some regions as soon as some water is extracted by evaporation. Further work is needed to understand the physical mechanisms and then to be able to predict the drying rate as a function of material characteristics and boundary conditions.

The understanding of this phenomenon is interesting for other reasons. Salt transport and crystallisation in building materials are responsible for major damages of modern structures [START_REF] Liu | The damage of calcium sulfoaluminate cement paste partially immersed in MgSO4 solution[END_REF] and cultural heritage [START_REF] Sawdy | A review of salt transport in porous media, assessment methods and salt reduction treatments[END_REF]. Similar effects in rocks can alter CO2 injectivity [START_REF] Peysson | Permeability alteration due to salt precipitation driven by drying in the context of CO2 injection[END_REF] and are at the origin of various geomorpholocial processes [START_REF] Rodriguez-Navarro | Salt Weathering: Influence of Evaporation Rate, Supersaturation and Crystallization Pattern[END_REF]. For example ions dissolved during imbibition can then, during drying, reach the solubility limit and precipitate as crystals below (subflorescence) or above (efflorescence) the sample free surface. These crystals may clog the porous structure [START_REF] Espinoza-Marzal | Impact of in-pore salt crystallization on transport properties[END_REF], lead to flaking and crumbling [START_REF] Kramar | Assessment of limestone deterioration due to salt formation by micro-Raman spectroscopy: application to architectural heritage[END_REF] or reduce the compression strength of the material [START_REF] Foraboschi | Experimental investigation on bricks from historical Venetian buildings subjected to moisture and salt crystallization[END_REF]. The mechanisms at the origin of damages have been discussed: crystallization pressure, supersaturation, solid volume change from thenardite to mirabilite [START_REF] Scherer | Stress from crystallization of salt[END_REF][START_REF] Flatt | Salt damage in porous materials: how high supersaturations are generated[END_REF][START_REF] Schiro | Damage Mechanisms of Porous Materials due to In-Pore Salt Crystallization[END_REF][START_REF] Flatt | Chemo-mechanics of salt damage in stone[END_REF][START_REF] Tsui | Crystallization damage by sodium sulfate[END_REF]; but it was also acknowledged that the remaining liquid films can play a significant role on crystallization processes [START_REF] Espinosa-Marzal | Advances in Understanding Damage by Salt Crystallization[END_REF]. Extensive studies [START_REF] Rodriguez-Navarro | Salt Weathering: Influence of Evaporation Rate, Supersaturation and Crystallization Pattern[END_REF][START_REF] Schultz | Influence of additives on crust formation during drying[END_REF] showed the variety of trends observed depending on solution properties (salt type and additives), but it was finally pointed out [START_REF] Rodriguez-Navarro | Salt Weathering: Influence of Evaporation Rate, Supersaturation and Crystallization Pattern[END_REF][START_REF] Shahidzadeh-Bonn | Damage in porous media du to salt crystallization[END_REF][START_REF] Desarnaud | Impact of the kinetics of salt crystallization on stone damage during rewetting/drying and humidity cycling[END_REF] from detailed information concerning salt and damage localization, that possible damages critically depend on the localization of salt crystals which requires the description of the dynamics of salt transport, spatial distribution and crystallization and its coupling with drying mechanisms.

For efflorescence the basic theoretical concepts of salt transport resulting from water flow during drying were confirmed by following the salt distribution by MRI [START_REF] Pel | Efflorescence pathway diagram: understanding salt weathering[END_REF][START_REF] Petkovic | Moisture and salt transport in three-layer plaster/substrate systems[END_REF]. It was then shown that crystals deposited at the sample free surface behave as a porous medium which pumps the liquid and brings ions around the most external liquid-air interface to form new crystals [START_REF] Veran-Tissoires | Discrete Salt Crystallization at the Surface of a Porous Medium[END_REF][START_REF] Sghaier | Effect of efflorescence formation on drying kinetics of porous media[END_REF]. A recent study showed in addition that the structure of the efflorescence layer can evolve with successive imbibition-drying cycles, which successively modifies the drying rate, increasing or decreasing it depending on the crystallization pathways or on the humidity of the crust [START_REF] Desarnaud | Drying of salt contaminated porous media: Effect of primary and secondary nucleation[END_REF]. Microtomography observations of the three phases around the sample top provided a detailed description of the process, which makes it possible to establish a model of crust formation dynamics as a function of drying rate [START_REF] Rad | Effects of grain and pore size on salt precipitation during evaporation from porous media: A pore-scale investigation[END_REF]. Besides it was shown that some unexpected effects (in drying rate variations) can result from the interplay between the crust and the evaporation conditions [START_REF] Gupta | How ferrocyanide influences NaCI crystallization under different humidity conditions[END_REF][START_REF] Gupta | Paradoxical drying of a fired-clay brick due to salt crystallization[END_REF].

Although it is a more damaging effect than efflorescence much less is known concerning the local mechanisms of subflorescence. It was shown from X-Ray microtomography that subtle local effects occur essentially at the scale of the pore [START_REF] Shokri | Pore-scale dynamics of salt transport and distribution in drying porous media[END_REF], and tend to induce pore clogging which reduces the drying rate, but these effects depend of the type of salt at work [START_REF] Espinoza-Marzal | Impact of in-pore salt crystallization on transport properties[END_REF]. From original experiments in 2D system it was also shown that a specific subflorescence pattern can grow because the crystal structure form a porous media which can drain a saturated salty solution up to the evaporating front, which eventually causes the drying rate to increase [START_REF] Sghaier | Evaporation-driven growth of large crystallized salt structures in a porous medium[END_REF]. Recently the work of Derluyn with the help of neutron radiography and X-ray microtomography provided new insights on salt transport, accumulation, and induced damage in limestone samples [START_REF] Derluyn | Characterizing saline uptake and salt distributions in porous limestone with neutron radiography and X-ray micro-tomography[END_REF][START_REF] Derluyn | Deformation and damage due to drying-induced salt crystallization in porous limestone[END_REF]. However our knowledge of the relationship between subflorescence and drying dynamics in porous media is still limited, essentially because it is more difficult to get a direct view of the processes at a local scale. Yet it governs the spatial distribution of crystals which in turn determines the position of possible damages, and this knowledge would open the way to the development of relevant tools for controlling or avoiding such effects. Taking advantage of MRI (Magnetic Resonance Imaging) data such an approach was recently developed for predicting the accumulation of initially suspended colloidal particles and its impact on drying dynamics [START_REF] Keita | MRI evidence for a receding-front effect in drying porous media[END_REF].

Here we aim at studying in more detail the process of subflorescence in plaster. From repeated imbibition-drying cycles on the same sample we observe that the drying rate is dramatically decreased as a result of crystal accumulation inside the sample. Through MRI and X-Ray microtomography we show that crystals are deposited around the first (i.e. most external) air-liquid interface, which induces its recession in the sample at a rate only depending on the current saturation (water to pore volume ratio) and the current total liquid volume in the sample. From this information we deduce a simple phenomenological model which predicts the drying dynamics of plaster under different conditions.

Materials and methods

We used a commercial α-plaster (Saint-Gobain Recherche, France) with a purity higher than 95% wt, a specific surface of 7500 g/m² (measured with the Blaine air permeability test), and a grain size between 10 and 20 m. The impurities are mainly clay (colloidal) particles. Our plaster paste samples are prepared according to a well-defined mixing protocol [START_REF] Seck | The mechanisms of plaster drying[END_REF][START_REF] Jaffel | Multi-scale characterization of evolving porous media[END_REF]. The initial water to powder ratio (w/p), i.e. the ratio of the mass of water to the mass of powder, was between 0.5 and 0.8. At the end of the hydration reaction (after about 25 min, as found from NMR measurements [START_REF] Seck | The mechanisms of plaster drying[END_REF]) which consumed about 20% of the initial water mass, the solid structure is filled with a calcium -sulphate solution at the gypsum solubility limit ( =2.65 g/l). This porous medium has initially a homogeneous porosity of 58% and a bimodal pore size distribution (from Mercury Intrusion Porosimetry) with peaks at 0.6 and 5 microns or so.

The plaster pastes are finally poured into a cylindrical container (35 mm or 54 mm diameter, and different heights from 12 to 40 mm) just at the end of the mixing protocol (i.e. initial time). The sample set up is then immediately dried by exposing its (open) top free surface to a tangential stream of dehydrated air (relative humidity < 1%) at ambient temperature (

1 C 23  
) through a long rectangular conduit of 60 x 3 mm² section. The sample has its free surface appearing from a hole in the continuity of the channel wall, and lies on a scale so that its mass is followed in time. The saturation  is computed from the ratio of the current water content to the initial water mass. This initial water mass is deduced from complete drying at 45°C after the test. Note that we choose this temperature of drying to avoid any possibility of decomposition of gypsum as it was indeed shown [START_REF] Jaffel | Multi-scale characterization of evolving porous media[END_REF] that at such temperature the chemical structure of the hydrates is unchanged. The drying dynamics is well described through the drying rate, i.e.

t d d 
, as a function of  . Such data are reproducible with an uncertainty of about 20% on the drying rate [START_REF] Seck | The mechanisms of plaster drying[END_REF] in particular due to the slight variations of the surface shape of the samples and the data processing to compute

t d d  from experimental data.
In order to look at the water distribution inside the sample in time during the test we inserted the whole drying set-up (without the scale) at the center of the magnet of a vertical proton MRI (Magnetic Resonance Imaging) system (Avance 24/80 DBX by Bruker, 0.5T superconducting magnet by Oxford, 20 cm inner diameter). The distribution of apparent water along the sample axis was measured using a one-dimensional double spin-echo sequence, i.e. two first echoes of the so-called Carr-Purcell-Meibom-Gill sequence [START_REF] Carr | Effects of Diffusion on Free Precession in Nuclear Magnetic Resonance Experiments[END_REF][START_REF] Meiboom | Modified SpinEcho Method for Measuring Nuclear Relaxation Times[END_REF][START_REF] Mcdonald | Surface relaxation and chemical exchange in hydrating cement pastes: a two-dimensional NMR relaxation study[END_REF] with extrapolation to compensate spin-spin relaxation.

Each measured value of this distribution corresponds to the total amount of water in a thin crosssectional layer 98 microns thick situated at a given position along the sample vertical axis. To get a satisfactory signal to noise ratio, the measurement time for each profile was two hours, so that a profile represents the average water distribution over that duration. Profile measurements were repeated every two hours.

In order to know how the internal structure evolves as a result of drying we used a RX-Solutions EasyTom tomograph equipped with an X-ray polychromatic microfocus source with a spot size of 5 µm (6,4W and 80 Kev) allowing a 3D imaging of the structure with cubic volume element of 5.8 µm edge.

Results

During the first drying after preparation the drying rate progressively decreases at the beginning of the desaturation then more rapidly decreases below 15 . 0   (see Figure 1). As already explained in [START_REF] Seck | The mechanisms of plaster drying[END_REF] this result contrasts with the well-known effect of constant drying rate down to a saturation of the order of 10% (at sufficiently low air flow velocity) for a simple porous medium filled with a pure hydrophilic liquid [START_REF] Brakel | Mass transfer in convective drying[END_REF][START_REF] Coussot | Scaling convective drying[END_REF]. Indeed, during drying simple materials exhibit a first regime, called the CRP (Constant Rate Period), during which the drying rate is approximately constant. During this period the saturation remains uniform in the material and the evaporation occurs from the upper liquid films very close to the free surface of the sample, while the liquid is drained by capillary effects through the film network throughout the sample. Then the drying rate starts to decrease significantly, this is the Falling Rate Period (FRP). In that regime the liquid flow driven by basic capillary effects is unable to provide liquid at a sufficient rate at the sample free surface, so that a thin dry region forms in this region which lowers the drying rate (since now the vapor has to diffuse through this region before reaching the sample free surface). The size of this dry region slightly increases during this regime. Finally a last regime takes place when the dry region now increases widely because there is no longer any possibility for the liquid to be drained towards the sample free surface. . A drying process should thus be followed from right to left.

0.0 0.5 1.0 10 -2 10 -1 10 0 (4) (3) (2)
Note that this effect of decrease of the drying rate observed with a plaster paste is not due the hydration reaction since after a relatively short time with regards to the total duration of the drying test (typically 40 hours) the hydration reaction ends, and the solid structure is fixed [START_REF] Seck | The mechanisms of plaster drying[END_REF]. Here the main difference between our sample and a simple porous medium as above described is the presence of ions (Ca 2+ , SO 4 2-) in the interstitial solution, which can precipitate as gypsum crystals and accumulate in some regions, thus reducing the permeability and the drying rate. Since the sample thickness remains constant and we do not see significant accumulation of crystals above the free surface, this means that we are dealing here with subflorescence. During most of the initial period of drying the liquid flows towards the regions of evaporation, which are mainly around the air-liquid interfaces the closest to the sample free surface because the air in the regions below is almost saturated in vapour [START_REF] Coussot | Scaling convective drying[END_REF]. The ions contained in the liquid are thus globally transported towards the free surface where the slight concentration increase due to each evaporation increment will readily induce crystallisation since the solution is at its solubility limit. This implies that we have a reduction of the pore size or the porosity at least in the upper regions of the sample. In porous media however, such effects usually only slightly reduce the duration of the initial constant rate period [START_REF] Brakel | Mass transfer in convective drying[END_REF][START_REF] Coussot | Scaling convective drying[END_REF]; indeed in this regime the rate of evaporation is imposed by the external air flow as long as the liquid can flow sufficiently rapidly through the medium towards the sample free surface as a result of capillary effects. This is so because the thickness of such a layer of small pores is so small that the pressure gradient imposed by capillary effects are easily extremely large and can easily move water through it. This means that in the present case something different, or more precisely some additional effect, happens.

In order to better appreciate the effect we repeated the test several times with the same sample after the initial drying. Our procedure consisted to first immerse the sample in a calcium -sulfate saturated solution, then put this system under vacuum. Under such conditions it seems unlikely that the initially solid gypsum can decompose, as can be observed under other conditions [START_REF] Mantellato | Reliable specific surface area measurements on anhydrous cements[END_REF], since it is not directly in contact with a gas at low pressure and instead is in contact with a saturated solution. Thus the initial system at the beginning of each new drying process was a sample filled with a calcium -sulphate saturated water solution. This situation differs from that of Desarnaud et al. [START_REF] Desarnaud | Drying of salt contaminated porous media: Effect of primary and secondary nucleation[END_REF] in which pure water is re-introduced in the porous sample (sandstone) and possible salt dissolution depended on the exact impregnation procedure.

All cycles share the same initial drying rate (see Figure 1) since each sample is saturated at the beginning of the test so that at the initial conditions of evaporation under the same flux are similar.

On the contrary drying rates observed after a short while are significantly lower for an increasing number of cycles, so that after five imbibition-drying cycles, for  in the range [0.15-0.8] the drying rate reaches values almost two orders of magnitude lower than the initial one (see Figure 1). It is remarkable that this dramatic decrease in drying rate is obtained with a total (after the five cycles) precipitated gypsum to sample mass ratio of 0.24 %. We also remark that after the first imbibitiondrying cycle the shape of the drying curve tends to a strong initial decrease followed by a plateau and finally a new strong decrease (see Figure 1). This suggests that the conditions of evaporation widely vary, which should be visible in the liquid distribution.

We can get further information from MRI experiments providing the evolution of the water distribution along the sample axis. At first sight (see Figure 2a) the results for the first cycle do not differ from those typically observed in simple porous media [START_REF] Seck | The mechanisms of plaster drying[END_REF]. The water distribution apparently remains homogeneous in the sample during a long first period (down to saturation of the order of 30%) (see Figure 2a), then we can see a dry front receding inside the sample while the plateau level behind goes on decreasing. Similar results are obtained for the subsequent cycles (see Figure 2b), although the time to reach each curve increases with the cycle number and some region of lower saturation forms at the beginning of drying below the sample free surface. This region of thickness much larger than the zone of accumulation of crystals is likely due to some heterogeneity of the sample. It does not evolve further during the subsequent cycles, and does not play any role in the evolution of the drying rate discussed in this paper.

The second stage of evolution of these profiles, i.e. large receding front, corresponds to the saturation range for which we have the second fast decrease of the drying. In this regime the liquid flow under capillary effects is not sufficiently rapid compared to the evaporation imposed by air flow along the sample free surface. The point is that in contrast, in the first stage the homogeneous saturation is usually associated with the constant rate period [START_REF] Faure | Drying of a model soil[END_REF], which is not the case here, in particular with the dramatic decrease in drying rate for the last cycles. There is thus an effect, which is not just a dry front of significant thickness, which remains to be found to explain this drying rate decrease. This suggests that some evolution of the porous structure plays a role. The aspect of the structure at different depths from the free surface as observed by X-ray microtomography is shown in Figure 3. These images are maps of the attenuation coefficient for X-Rays in a thin layer of thickness 5.8 µm situated at a given position along the sample axis. In such pictures, since initial solid matrix and subsequently deposited crystals are of the same chemical composition and density they exhibit the same X-ray attenuation coefficient, so that the grey scale is linearly related to the porosity, i.e. black for a void volume of thickness larger than 5.8 µm, white for a pure solid volume. At depths beyond a few hundreds of microns the structure keeps the aspect shown in Figure 3c, with a large number of relatively dark regions of various size, which indicates a large number of voids. At the approach of the free surface the number of dark regions decreases, while the amount of bright regions is higher; this indicates that the becomes globally less porous (see Figure 3a,b), the pores apparently become smaller and less numerous, an effect which obviously results from the deposition of gypsum crystals.

In order to quantify this effect we computed the average porosity for each image along the vertical direction. The examples of Figure 4 show that after a rapid decrease, which corresponds to the crossover of the sample roughness, the porosity reaches a minimum value associated with a high concentration of crystals, then it increases rapidly and reaches a plateau. The deposition of crystals occurs preferentially where the evaporation occurs, since, there, by water removal, the critical dissolution concentration tends to be overcome first, and from MRI we know that most of the evaporation occurs close to the free surface so that most crystals are formed around the sample top. Some crystallization may occur while the receding front develops over large depths inside the sample but in that case the remaining ion content is small and the crystals will be dispersed over a large distance. We deduce that the plateau essentially corresponds to the unchanged porous structure.

The porosity curves vary from one cycle to another essentially by a variation of the peak (minimum) size, which increases for increasing cycle number. This result is qualitatively close to that obtained by Derluyn et al. [START_REF] Derluyn | Characterizing saline uptake and salt distributions in porous limestone with neutron radiography and X-ray micro-tomography[END_REF] with sodium sulfate crystalizing in drying limestone. This means that the pores are filled more and more at each step, essentially below the sample free surface. Note that the curve for the third cycle does not exhibit a shape quite consistent with the general evolution of the other curves. This is likely due to the uncertainty on data (which is hard to quantify). More generally, considering the (complex) way these data are obtained we cannot consider that they provide precise values of the effective porosity, the most relevant results concern the variations of this apparent porosity with the depth in the sample. As a consequence, in the following we will essentially focus of these variations (through a probability density function on crystal location, see below). Figure 4: Average porosity (as deduced from X-Ray tomography) in elementary cross-sectional layers along the sample axis as a function of the depth in the sample for the successive dryingimbibition cycles (same conditions and colors as in Figure 1).

Discussion

General data analysis

In order to look consistently at the evolution of the profiles of porosity as obtained in Figure 4 we can withdraw the plateau level from each of these curves, and rescale them by their integral. Note that in this representation we choose the initial abscissa as the mid position in the sample roughness. The resulting curves (see Figure 5) are proportional to the probability density function ( p ) on crystal location, formed since the beginning of the series of cycles, along the downward vertical direction ( x ). Thus the probability to have a crystal between x and

dx x  is dx x p ) ( .
Finally, it appears that except in the sample roughness, the density function can be very well represented by a decreasing exponential function:

      x p   exp 1 ( 1
)
where x is the depth in the sample (distance from the bottom of the "roughness" zone at free surface) and  a characteristic length. This similar distribution of crystals in the sample for the successive imbibition-drying cycles is rather surprising since for each cycle the porous structure evolves and tends to clog further as a result of the deposition of these crystals. From the ion concentration (leading to total crystal mass of 0.055 g per cycle) and the sample porosity we can compute the effective volume occupied by the crystals. We find that had all the crystals be compacted in the pores they would fill the sample over a thickness of about 40 m after the first drying. This means that for the distribution observed (see Figure 5), extending over a much larger distance, the porosity is far from being filled after the first drying. However it then progressively fills during the next cycles. Indeed after the last cycle the total amount of crystals deposited would saturate the pores over a sample thickness of 200 m. Since a large fraction of crystals are in fact distributed in a region with a thickness of this order of magnitude close to the free surface of the sample (see Figure 5) this means that at that time the porosity is very low in this region, which is qualitatively consistent with the results of Figure 4.

Since crystals are expected to deposit mainly around the top liquid-air interface along which most of the evaporation occurs, and since the saturation in the sample decreases as long as crystals deposit, the crystal density decreasing downwards suggests that the liquid interface withdraws (i.e. advances downwards) progressively while leaving crystals behind. It is possible to check this assumption by coming back to the MRI profiles and zooming on the sample top. We effectively see the first air-liquid interface slightly receding in the sample (see insets of Figure 2a and2b). This effect will induce a decrease of the drying rate with some analogy with the impact of a dry crust formed after efflorescence [START_REF] Desarnaud | Drying of salt contaminated porous media: Effect of primary and secondary nucleation[END_REF]. Note that we cannot know the exact extent of the dry region which forms since the profiles correspond to time and space averages, but the progression of the front anyway means that on average the distance between the top liquid-air interfaces and the sample free surface increases approximately as deduced from these profiles. For the following computations we arbitrarily consider that the position of the first air-liquid interface is the profile position for an NMR signal of in our arbitrary units (this is approximately twice the value of the noise on the data). Also note that this recession of the top air-liquid interface does not correspond to a situation of receding front as it occurs in the last stage of drying with hydrophilic material, or immediately after the beginning of drying in a hydrophobic material, when the saturation decreases strictly as a result of the front advance. Here, even if we locally have a slight receding front effect we still have a major effect of homogeneous decrease of the saturation throughout the sample, as shown for example in Figure 2. Finally it is remarkable that the saturation decreases with the front position in a way very similar to the crystal density distribution (see Figure 5), and this evolution may be well described by a similar exponential model:

) exp(   x   (2)
The most striking phenomenon is that all these characteristics are similar for each new cycle: after rescaling them by their integral over the whole thickness (i.e. the expected total amount of crystals at each step) the crystal distribution curves very well superimpose and may be represented by the same exponential model (see Figure 5). That means that successive imbibition-drying cycles simply lead to add a new distribution of crystals of same concentration at each depth in the sample. This is quite surprising since in parallel the apparent drying conditions, namely the drying rate, strongly change (drop of the drying rate in the first regime after each cycle). However if we now turn to the evolution of the interface motion as seen from MRI and represented in terms of saturation as a function of the depth we discover that it is also similar for each cycle although its dynamics obviously strongly differs from one cycle to another and finally covers a two decade range of drying rate. This demonstrates that the position of the receding front only depends on the saturation, and the above law (equation ( 2)) describing the relationship between the saturation and the front position is general for our material, thus the parameter  is constant in our range of experimental conditions. This phenomenon is likely explained by a kind of dewetting resulting from crystal formation. It was indeed suggested that different types of interactions of the solution and the crystals with the solid surface can induce different types of processes, such as a recession of the line of contact leading to subflorescence [START_REF] Shahidzadeh-Bonn | Salt crystallization during evaporation: impact of interfacial properties[END_REF]. More precisely it was shown [START_REF] Shahidzadeh | Salt stains from evaporating droplets[END_REF] that the drying of NaCl solution droplets significantly differs from calcium sulfate solution. For the former the crystals formed around the contact line tend to progressively move towards the interior of the droplet as a result of capillary effects, whereas for the latter the crystals forms in the bulk [START_REF] Shahidzadeh | Salt stains from evaporating droplets[END_REF][START_REF] Wang | Confinement leads to control over Calcium Sulfate Polymorph[END_REF] but for sufficiently large contact angles of the solution on the solid surface, the crystals remain stuck to the solid surface which finally induces a retraction of the liquid and leads in the end to a "coffee stain effect". Moreover, it was shown, in consistence with our observations, that in drying droplets the crystallization patterns are independent of the drying rate [START_REF] Shahidzadeh | Salt stains from evaporating droplets[END_REF]. These observations might be the key explanation of the front recession associated with crystal formation inside plaster.

A consequence of these results is that the front position does not govern by itself the drying rate, since for the same position of this front we observe widely different drying rates for different cycles. It is more likely a change of the diffusion coefficient in the dry region which plays the major role: the vapour coming from the top liquid-air interface and going to the free surface has to find a path through a more tortuous and with smaller pore sizes porous medium as more crystals are deposited in the structure.

Modelling plaster drying

From these data a simple description of the process emerges. Note that a sophisticated theoretical description allowing to predict salt deposits, water evaporation, stress distribution and damages during this kind of experiment was developed in [START_REF] Derluyn | Deformation and damage due to drying-induced salt crystallization in porous limestone[END_REF]. Here our approach focuses on salt deposit and water transport and only relies on basic physical assumptions in consistency with above observations. Let us assume that during drying the sample progressively desaturates in the form of a "rectangular" liquid distribution with a decreasing amplitude and decreasing extent: a dry region forms down to an increasing depth x , while the saturation ( ) in the wet region decreases homogeneously. Note that from our data (see insets of Figure 2) we see that this is an approximate description of reality: the front of the liquid distribution is not perfectly straight. For a step of advance of the front dx associated with a change of the mean saturation in the wet region from  to   d



we have a loss of water volume equal to

   dxS d S x L d      ) (
, in which L is the sample height, S its section and  its porosity. This leads to the formation of a mass of crystals

 d pdx   (3)
It is interesting to note that this equation predicts our experimental result concerning the link between p and  , as it appears from equations ( 1) and (2) fitted to data, i.e.

dx d p    .
In order to find the law of variation of the saturation or the density function of the crystals we now have to predict the advance of the top liquid-air interface. The above analysis from our data showed that this process mainly depends on the local conditions and not on the drying rate, and it was suggested that it is associated with some dewetting effect. Such a phenomenon depends on the exact shape of the liquid-solid-air interface, which is a priori very complex and in particular depends on the crystal shape which is hardly predictable. We nevertheless assume that this effect can be described in a statistical way in a cross-section of the porous medium so that on average some characteristic step of dewetting would occur after the deposition of some amount of crystals, which may be seen as a number of opportunities of front withdrawal. On the other hand, as the saturation decreases the amount of interfaces to recede decreases. Finally we expect to have an effective withdrawal dx proportional to the volume of crystals deposited between x and dx x  , i.e. , where k is a factor which depends on the physical characteristics of the crystals and the porous structure. We then deduce :

 pdx

 p L k     (4)
From ( 3) and ( 4) we obtain    Ld k dx   which, taking into account the boundary condition

1 ) 0 (   x 
, can be integrated and give the expression of the type (2) with

   L k  (5)
Our model thus predicts the exponential variation of the saturation and the value for the coefficient k may be found from our data and we obtain

7 . 3  k
. It is interesting to note that  is proportional to the sample height. This means that, although we did not explore that aspect in the present study, the spreading of the crystals in the sample should be larger for longer samples.

In this context it remains somewhat surprising that we obtain the same law (i.e. equation ( 2) with the same value for  ) for the successive drying-imbibition cycles whereas the porous structure likely changes as more and more crystals are deposited in a limited volume (which tends to be filled just below the sample free surface). This suggests that k essentially depends on the crystal characteristics rather than on the porous structure. Let us emphasis that k=1 would correspond to the formation of a compact crystal layer with the same thickness as that of the current wetting liquid film. By comparison, k=3.7 here describes 3.7 times less dense of sparser crystallization.

We can elaborate a simple model predicting the drying rate. As usual we assume that air is saturated with vapor below the first air-liquid interface from the free surface of the sample situated at CaSO solution [START_REF] Berdugo | On the behaviour of the Ca-SO4-H2O system[END_REF]. The solution for the vapor flux (V ), expressed in terms of the equivalent volume of liquid water crossing the sample section per unit time, is

     0 0 0 d 1 1 1 x x D n V   (6)
For the factor  which expresses the effect of the new crystals on the tortuosity and porosity, and as a first approximation is inversely proportional to the characteristic length of the path vapor has to follow before reaching the sample free surface, we make the following simple assumption: each new crystal formed (associated to a mass increment dm ) tends to clog a possible path for vapor which now has to find another path through the rest of the porous structure. This means that the characteristic path length is increased proportionally to the current path length: dm d

) 1 ( ) 1 (    . We deduce   ) ( exp 0 x p m     
, with 0 m the total (additional) crystal mass in the sample. Note that this approach can be used only when the volume of additional crystals is smaller than the pore volume.

This model appears capable to well predict the shape of the variations of the drying rate in the first regime (i.e. before the strong drying rate decrease) for different initial water to solid ratios by adjusting the only parameter  (see Figure 6), and taking  and  as constant.  appears to decrease when the water to solid ratio increases but we have no explanation on the physical origin of this variation.

It is remarkable that with the parameters  and  fitted above the model finally reproduces rather well the drying curves (see Fig. 1), not only the decrease during the first stage but also the different plateau levels after each cycle. Thus this model captures the most important physical characteristics of the process. Finally for a low saturation a receding front, of thickness larger than that of the region of accumulation of crystals, starts to develop as a result of a balance capillary and hydrodynamic effects [START_REF] Coussot | Scaling convective drying[END_REF] and the second regime of strong decrease of the drying rate is reached (see Fig. 1). It also appears possible to predict the shape of the drying curve for a sample submitted to different air fluxes (see Figure 7). Note that in that case the initial drying rate taken into account in the model was fitted to the first value in the regime of slow decrease observed just after the initial drop. This initial drop is usually observed in such drying process and is likely due to the disappearance of the thin film of liquid initially covering the sample free surface and the first step of desaturation of the porous medium positioning the first air-liquid interfaces at a depth of the order of a particle radius. However the comparison of model predictions with data for sample of different heights does not show an excellent agreement (see Figure 8). In that case the model parameters are those determined from data of Figure 6 for a sample of height 12 mm which are directly used to compute model predictions for other heights under similar experimental conditions. Although the shape and decrease of the drying rate curves are qualitatively well predicted there remains a significant quantitative discrepancy (see Figure 8). A possible explanation is that fitting the model to data obtained with small sample heights, and for the first drying process only, might be too sensitive to specific values obtained in that case. Indeed in that case the total crystal amount is much lower than that involved in successive drying-imbibition cycles for large sample heights, and the model fitting is likely more robust. ) for an initial water to powder ratio of 0.6 (data from [START_REF] Seck | The mechanisms of plaster drying[END_REF]) at different velocities of the air flux (sample height 12 mm):

Conclusion

This work shows that it is possible to quantify, on average, the effects which occur at a local scale inside the sample in terms of water transport and crystal formation, which then allow to understand the dynamics of drying coupled with crystallization in porous media. The simple physical mechanics observed provides a full frame for describing the drying properties of plaster, and are likely applicable to other situations. However the parameters of the model, i.e.  ,  and  depend on the characteristics of the ions and the porous structure, and their exact physico-chemical origin would need to be studied with the help of experiments at a local scale. Our approach nevertheless provides some basic understanding of subflorescence dynamics in building materials, and suggests that it should be possible to avoid subflorescence by maintaining a wetting up to the sample free surface with appropriate additives in the liquid.

Figure 1 :

 1 Figure 1: Drying rate (rescaled by the initial one, i.e. 2.8 g/hr) as a function of saturation for the plaster for an increasing number (from top to bottom, numbers in brackets) of imbibition-drying cycles after the paste preparation (sample height 40 mm, w/p=0.8, air flux velocity: 7.3 m/s): (continuous lines) experiments, (dotted lines) model (equation 2) with 4 10 48    and

Figure 2 :

 2 Figure 2:Water distribution profiles in a plaster at successive times (every two hours from top to bottom) during the first drying after preparation (a) and after the first cycle (b). The experimental conditions are those described in Figure1. The insets show a zoom on these profiles around the sample free surface.

Figure 3 :

 3 Figure 3: Views from tomography (with the same contrast) of the structure of the plaster sample after the 4 th cycle at different depths ( x ) from the free surface: (a) 40 m, (b) 190 m, (c) 1.7 mm. The experimental conditions correspond to those of Figure 1.

Figure 5 :

 5 Figure 5: Probability density functions on crystal location estimated from Microtomography as a function of the depth in the sample for a plaster sample after different imbibition-cycles (continuous lines with same colors as in Fig.1). Mean sample saturation (colored squares) as a function of the position of the water front (from MRI profiles, see text). Dashed line: exponential model (equation 1) with μm 150  

  density), divided by the average cross-section area occupied by the liquid in the pores. We thus obtain dx

  in contact with the external air flow and the evaporation may be described in the form of simple diffusion of vapor with a diffusion coefficient D through a boundary layer of thickness [START_REF] Keita | MRI evidence for a receding-front effect in drying porous media[END_REF]of the porous structure which reflects the decrease of porosity and the increase of tortuosity when the crystal concentration increases. In that case the diffusion equation writes 0

Figure 6 :

 6 Figure 6: Drying rate (rescaled by the initial one, i.e. 4 g/hr) as a function of saturation and model (equation 6) (continuous line) fitted to data for plaster (first drying after preparation) at different initial water to powder ratios (sample height 12 mm) (data from [13]): 0.8 (green), 4 10 48    ,

Figure 7 :

 7 Figure 7: Drying rate as a function of saturation and model (equation 6) (continuous line) fitted to data for plaster (first drying after preparation) (with

Figure 8 :

 8 Figure 8: Drying rate (rescaled by the initial one, i.e. 3.4 g/hr) as a function of saturation, and model (equation 6, with


Zhanle Yu

Christophe Eloy

Extension of Lighthill's slender-body theory to moderate aspect ratios

Keywords: Slender-body theory, Vortex panel method, Flutter instability

Calculating the fluid forces acting on a moving body at high Reynolds number is crucial in many fluidstructure interaction problems, such as fish swimming or flutter instabilities. To estimate these forces, Lighthill developed the slender-body theory, which assumes a potential flow and an asymptotically small aspect ratio. Yet, it is still unclear whether Lighthill's theory is still valid for aspect ratios of order one. To address this question, we solve numerically with a panel method the full three-dimensional problem of a rectangular plate deforming periodically in a potential flow. These numerical simulations are used to calculate the pressure jump distribution across the plate for different aspect ratios. We find that numerical simulations and slender-body theory give similar results far from trailing edge. Close to the trailing edge however, there is a discrepancy, which is due to the use of a Kutta condition in the simulations (i.e. no pressure jump at the trailing edge), while, in the slender-body theory, the pressure jump is non zero. We propose a simple extension of Lighthill's slender-body theory that accounts for this discrepancy. The usefulness of this extension is then discussed and illustrated with a generic fluid-structure interaction problem and with the flag instability problem.

Introduction

Slender-body theory has first been developed by [START_REF] Munk | The Aerodynamic Forces on Airship Hulls[END_REF] to study aerodynamic loads exerted on an airship hull. [START_REF] Munk | The Aerodynamic Forces on Airship Hulls[END_REF] considered an elongated rigid body moving in a fluid at rest. This theory is essentially based on potential flow theory and conservation of fluid momentum along the body axis. It implicitly assumes that the flow generated by each section does not affect the other sections of the body, which is guaranteed if the body is sufficiently slender and if the section does not vary abruptly.

In the context of fish swimming, [START_REF] Lighthill | Note on the swimming of slender fish[END_REF] later extended slender-body theory to calculate the fluid forces exerted on a deforming body. In the present paper, we will refer to this extension as Lighthill's slender-body theory, or simply slender-body theory when there is no ambiguity. Its hypotheses are the same as in [START_REF] Munk | The Aerodynamic Forces on Airship Hulls[END_REF]: the flow is assumed to be potential and incompressible, the body aspect ratio, defined as the ratio between crosswise and streamwise dimensions, is assumed to be asymptotically small. With these hypotheses, again, pressure forces exerted on the body can be calculated using an argument of momentum conservation. Remarkably, [START_REF] Lighthill | Note on the swimming of slender fish[END_REF] showed that when a fish deforms periodically, the average total force only depends on the kinematics of the tail.

Slender-body theory is limited to displacements that are small compared to the streamwise body length. To circumvent this limitation, [START_REF] Lighthill | Large-amplitude elongated-body theory of fish locomotion[END_REF] extended slender-body theory, a theory known as Lighthill's elongated-body theory. Recently, this elongated-body theory was itself extended by [START_REF] Candelier | Three-dimensional extension of Lighthill's large-amplitude elongated-body theory of fish locomotion[END_REF] and [START_REF] Candelier | Note on the swimming of an elongated-body in a non-uniform flow[END_REF] to respectively account for arbitrary three-dimensional motions (including all components of rotation and translation) and a non-uniform background flow.

An important field of application of Lighthill's slender-body and elongated-body theories has been the hydrodynamics of fish swimming, as first advocated by [START_REF] Lighthill | Note on the swimming of slender fish[END_REF][START_REF] Lighthill | Large-amplitude elongated-body theory of fish locomotion[END_REF]. In this context, these theories have been used to better understand the kinematics and energetics of swimming, and in particular to assess whether fish employ smart drag reduction mechanisms or if, on the contrary, drag is enhanced by fish motions [START_REF] Alexander | Swimming[END_REF][START_REF] Webb | Hydrodynamics and energetics of fish propulsion[END_REF][START_REF] Videler | Swimming movements, body structure and propulsion in cod, Gadus morhua[END_REF][START_REF] Ehrenstein | Skin friction on a moving wall and its implications for swimming animals[END_REF]. To give an unbiased review of the field, we should also mention that some authors criticized Lighthill's theories, arguing that these potential approaches tend to overestimate forces [START_REF] Hess | Fast continuous swimming of saithe (Pollachius virens): a dynamic analysis of bending moments and muscle power[END_REF][START_REF] Anderson | The boundary layer of swimming fish[END_REF][START_REF] Shirgaonkar | A new mathematical formulation and fast algorithm for fully resolved simulation of self-propulsion[END_REF]). Lighthill's slender-body has also been applied to engineering problems, such as the flutter instability [START_REF] Lemaitre | Instability of a long ribbon hanging in axial air flow[END_REF], or energy harvesting [START_REF] Singh | Energy harvesting from axial fluid-elastic instabilities of a cylinder[END_REF][START_REF] Michelin | Energy harvesting efficiency of piezoelectric flags in axial flows[END_REF]. Note that Lighthill's slender-body theory should not be confused with the low-Reynolds-number slender-body theory (e.g., [START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF], even if -to add to the confusion-Lighthill also contributed to this field [START_REF] Lighthill | Flagellar hydrodynamics: the John von Neumann lecture[END_REF].

A significant limitation of Lighthill's slender-body theory is that it is only applicable to bodies of asymptotically small aspect ratios. In the context of the flutter instability, one of us addressed this limitation by developing a spectral approach [START_REF] Eloy | Flutter of a rectangular plate[END_REF][START_REF] Eloy | Aeroelastic instability of cantilevered flexible plates in uniform flow[END_REF]. In essence, the potential flow problem is first solved in the Fourier space for any aspect ratio; then, an inverse Fourier transform is applied to come back to real space. Unfortunately, this approach is limited to body with a rectangular planform and its implementation is not as simple as Lighthill's slender-body theory. The objective of the present study is to examine a alternative analytical method to treat bodies with moderate aspect ratios.

Another important drawback of Lighthill's slender-body theory is that it is not compatible with a Kutta condition at the trailing edge. Indeed, the Kutta condition [START_REF] Crighton | The Kutta condition in unsteady flow[END_REF], even for unsteady flows, requires that the pressure jump across the lifting surface is zero at the trailing edge. In Lighthill's slenderbody theory, not only the force at the trailing edge is finite, but its value is the only relevant quantity if one is interested in the average total force for a periodic motion.

This paper is organized as follows. In §2, we will recall how the lifting-surface problem is stated and we will describe an asymptotic solution of this problem valid in the limit of small aspect ratios. In §3, we will present a panel method in the frequency domain, aimed at solving the same problem numerically. An alternative analytical approach, also valid in the limit of small aspect ratios and called the approximate solution, will then be presented in §4. The mathematical derivation of this approximate solution may be less rigorous than the derivation of the asymptotic solution, but it is much simpler to implement and it provides an straightforward and useful correction to slender-body theory. The merits of these two different methods will be compared to slender-body theory and to the numerical panel method in §5 through a generic fluid-structure interaction problem and through the analysis of the flutter instability of a cantilever flexible plate. Finally, these results will be summarized and discussed in §6.

Asymptotic solution of the lifting-surface integral

Consider a uniform, incompressible, and irrotional flow of velocity U over a flexible plate deforming periodically with angular frequency Ω (Fig. 1). In a Cartesian coordinate system (OXY Z), the incoming fluid of density ρ is flowing along X, and the flexible surface is deforming around the reference plane Z = 0. The periodic deformation of the plate is described by the complex displacement Z(M 0 , T ) = h(M 0 )e iΩT , where M 0 is a point in the plane Z = 0. The displacement is assumed sufficiently small such that (1) vortex sheets attached to the flexible surface and shed in the wake can be represented by their projections onto the reference plane, respectively S and Σ (Fig. 1); (2) the flow induced by the motion of the flexible plate can be considered as a periodic perturbation of the incident flow. For the sake of simplicity, the complex factor e iΩT will be omitted for all unsteady quantities. The point M 0 will represent an arbitrary point in the plane S, the point M will be in S or Σ, and Q will be a point in the volume, outside S and Σ.

Because we assume a potential flow, a perturbation potential Φ can be defined such that the total velocity field is given by U + ∇Φ. Because the flow is incompressible, Φ obeys a Laplace equation with a kinematic Neumann boundary condition given by the impermeability on the plate

∆Φ(Q) = 0, with ∂ Z Φ(M 0 ) = W (M 0 ) = D(h(M 0 )), and D = iΩ + U ∂ X , ( 1
)
where W is the imposed Z-component of the flow velocity, called the downwash, and D is the convective derivative. Note that the Z-derivative of Φ has no definition on S, hence, we define where M 0 is the projection of Q onto the reference plane. Physically, this definition simply states that the Z-component of flow velocity is a continuous function, including in Z = 0. Equation ( 1) can be solved by means of the Green's representation theorem [START_REF] Morse | Methods of theoretical physics[END_REF]. For any point Q, the potential Φ can be written as

∂ Z Φ(M 0 ) = lim Z Q →0 ± ∂ Z Φ(Q), (2)
X Y Z Σ S 2L 2H U h(M 0 )e iΩT
Φ(Q) = - S+Σ [Φ](M ) ∂ ∂Z M G(M, Q)dS M , (3)
where [Φ] is the jump of Φ across S, and G(M, Q) is the Green's function of the three-dimensional Laplace equation:

[Φ](M ) = Φ(X M , Y M , 0 + ) -Φ(X M , Y M , 0 -), (4)
G(M, Q) = - 1 4π|M Q| . ( 5
)
The Neumann boundary condition is imposed by applying Z-derivative operator onto Eq. (3) and taking the limit Z Q → 0 ± . This yields

W (M 0 ) = lim Z Q →0 ± S+Σ [Φ](M ) ∂ 2 ∂Z Q ∂Z M G(M, Q)dS M , (6)
where M 0 is the projection of Q onto the reference plane. The existence of this limit is guaranteed physically, and it can be evaluated by exchanging the limit and integral operators

W (M 0 ) = × S+Σ [Φ](M )K(M, M 0 )dS M , (7)
with

K(M, M 0 ) = lim Z Q →0 ± ∂ 2 ∂Z Q ∂Z M G(M, Q) = 1 4π|M M 0 | 3 . ( 8
)
The above exchange of operators is only permitted when the integral is taken in the finite-part sense of [START_REF] Hadamard | Lectures on Cauchy's problem in linear partial differential equations[END_REF] (represented by the cross on the integral sign). In this case, the singular integral in Eq. ( 7), usually called the lifting-surface integral exists and can, in principle, be calculated, as it has been shown by [START_REF] Mangler | Improper Integrals in Theoretical Aerodynamics[END_REF] and [START_REF] Brandao | Improper integrals in theoretical aerodynamics -The problem revisited[END_REF].

Assuming the potential jump [Φ], solution of Eq. ( 7), is found, the pressure jump across the plate, [P ], can be calculated with the linearized unsteady Bernoulli equation

[P ](M 0 ) = -ρD[Φ](M 0 ). ( 9
)
The downwash W is related to the unknown potential jump [Φ] through a lifting-surface integral (7). Mathematically, this is a Fredholm integral equation of the first kind with a hypersingular kernel K(|M M 0 |). The inversion of Eq. ( 7) is the essential part of a lifting-surface problem. For a lifting surface of arbitrary geometry, the analytic inversion is generally impossible. In this case, numerical techniques have to be used. Popular approaches include panel methods, such as the vortex-lattice method (Katz and Plotkin, 2001), or the doublet-lattice method [START_REF] Albano | A doublet-lattice method for calculating lift distributions on oscillating surfaces in subsonic flows[END_REF]. However, if the aspect ratio of the lifting surface is asymptotically small or large, specific analytic methods such as the Matched Asymptotic Expansion can be used [START_REF] Guermond | A systematic formula for the asymptotic expansion of singular integrals[END_REF][START_REF] Kida | An alternative treatment of lifting-line theory as a perturbation problem[END_REF][START_REF] Eloy | A unified introduction to fluid mechanics of flying and swimming at high Reynolds number[END_REF].

In the limit of asymptotically small aspect ratios, the leading order and the first order correction of the force (per unit chord length) exerted by the plate onto the fluid have been explicitly given in Eloy et al. ( 2010)

F 0 (X) = D(ρπH 2 W (X)), (10)
F 1 (X) = D ρπH 4 4 ln 4 A - 3 4 W (X) + × ∞ -L W (X M ) |X -X M | 3 dX M , ( 11
)
where H is the semi-span, L the semi-chord, and A = H/L the aspect ratio. In the following, we will refer to this approach as the asymptotic solution. The zero th order expression in Eq. ( 10) is strictly equivalent to Lighthill's slender-body theory, although the derivation proposed by [START_REF] Lighthill | Note on the swimming of slender fish[END_REF] is far more elegant.

The factor (ρπH 2 ) represents the added mass of fluid per unit chord, and Eq. ( 10) simply expresses that the force is the rate-of-change of the fluid momentum. One of the objectives of the present paper is to assess whether the first order correction given by Eq. ( 11) improves the accuracy of the fluid force, especially for lifting surfaces of moderate aspect ratios. To do so, we need a reference value of the lifting force that will be computed with a panel method, as explained in the next section.

Panel method in the frequency domain

The panel method is a numerical method that can be used to invert the lifting-surface integral equation and calculate the aerodynamic loads. Since we assume here that the displacement of the moving surface is periodic and small, some simplifications can be made to reduce the computational costs. More details about the classic panel method can be found in Katz and Plotkin (2001).

First, the boundary layer and wake vortex sheet are modeled by vortex rings placed on S and Σ, respectively. The impermeability condition is enforced on S instead of the actual moving surface. The discretizing grid is shown in Fig. 2. The indices i, j, or sequential counter k = (i -1) × N s + j, are used to track the location of a vortex ring (or a collocation point).

Second, the unsteady quantities are expressed with the complex factor e iΩT . Thus, the normal vector n i,j,T defined at the collocation point can be written as

n i,j,T (- ∂h ∂X e iΩT , - ∂h ∂Y e iΩT , 1). ( 12
)
The circulation of a vortex ring is equal to

Γ i,j,T = Γ i,j e iΩT , 1 ≤ i ≤ N c , 1 ≤ j ≤ N s . (13)
In order to fulfill the Kutta condition and enforce zero pressure jump across the wake vortex sheet Σ, the circulation of a wake vortex ring is given by

Γ i,j,T = Γ Nc,j e iΩ(T -(i-Nc)∆c/U ) , N c + 1 ≤ i ≤ N w , 1 ≤ j ≤ N s , ( 14
)
where ∆c is the chordwise length of vortex rings (Fig. 2). The above equation implies that wake vortices are first generated at the trailing edge and then advected downstream with velocity U without changing circulation.

X Z U O Y N c N w ∆b Γ i j i i-1 n i j ∆c ∆c Γ M j Γ M j e -i Ω ∆ c / U Γ M j e -i Ω 2 ∆ c / U N s j j+1 ∆c/4 ∆c N c N c N c
Figure 2: Representation of the numerical panel method used to solve in the frequency domain the lifting-surface equation.

As a consequence of previous simplifications, the complex factor e iΩT in the linear system can be omitted, and the impermeability condition for all collocation points on S becomes independent of time. The discrete form of the system is

M a Γ = R, (15)
where M a is the influence matrix, the vector Γ contains the unknown circulations of all vortex rings on S, and the elements of R are the values of the downwash at the collocation points (see Eq. 1). Equation ( 15) can be viewed as the discrete form of the lifting-surface integral (7).

The coefficients of the influence matrix M a are given by

a kl = a * kl (if l is not on the trailing edge), a * kl + b kl (if l is on the trailing edge), (16)
where the details of the coefficients a * kl and b kl can be found in Katz and Plotkin (2001). The coefficients a * kl represent the Z-component of the flow velocity at the collocation point k induced by the vortex ring l with unit circulation, and the coefficients b kl are the Z-component of the velocity induced by all wake vortex rings, which result from the advection of the trailing edge vortex ring l (see Eq. 14).

The circulations of vortex rings on S can be computed by inverting the linear system (15). We then calculate the pressure jump by applying a discretized version of the linearized Bernoulli equation ( 9). The advantage of the present approach is that the whole computation needs to be executed only once for a periodic lifting-surface problem (there is no time-marching). The computational costs are thus reduced compared to a classic unsteady panel method (Katz and Plotkin, 2001). This approach will be used to compute the aerodynamic forces exerted on a moving flexible plate. These computations will serve as reference values to assess the validity of our analytic asymptotic approaches.

Approximate solution of the lifting-surface integral

To begin with, we re-formulate the lifting-surface integral in terms of the acceleration potential Ψ = -P/ρ = D(Φ) (we simply use the linearized unsteady Bernoulli equation ( 9)). Applying the linear differential operator D to Eq. (1) yields a Laplace problem for Ψ with Neumann boundary conditions

∆Ψ(Q) = 0, with ∂ Z Ψ(M 0 ) = D(W (M 0 )). ( 17
)
As described in section 2, the solution of Eq. ( 17) can be found by using the Green's representation theorem [START_REF] Morse | Methods of theoretical physics[END_REF], and the Neumann boundary condition can be enforced by applying the Z-derivative operator to Ψ and taking the limit Z Q → 0. Using this method, the lifting-surface integral is expressed in terms of the acceleration potential jump

D(W (M 0 )) = × S [Ψ](M )K(|M M 0 |)dS M , (18)
where [Ψ] is the acceleration potential jump across S, and K is the kernel function given by Eq. ( 8). Note that, contrary to the lifting-surface integral for the pressure jump given by ( 7), here the integral can be taken over S only, and not S + Σ, because the acceleration potential jump is zero in the wake. However, just like Eq. ( 7), there is no general method to invert this Fredholm integral problem. In this section, we will seek an approximate solution instead.

In order to facilitate the solving procedure, the following dimensionless variables are introduced

x = X M 0 L , ξ = X M L , y = Y M 0 H , η = Y M H , [ψ] = [Ψ] U 2 , w * = D(W )L U 2 , t = T U L , ω = ΩL U . ( 19
)
Separating the chordwise and spanwise integration, Eq. ( 18) can be written

w * (x, y) = A 4π × 1 -1 I(x, )dη, (20)
with

I(x, ) = 1 -1 [ψ](ξ, η) [(ξ -x) 2 + 2 ] 3/2 dξ, (21)
where = A(η -y) is asymptotically small. Assume [ψ](ξ, η) is analytic with respect to ξ, it can be approximated by a Taylor series

[ψ](ξ, η) = [ψ](x, η) + ∂ x [ψ](x, η)(ξ -x) + • • • (22)
Because its kernel is large when ξ ≈ x, the integral I(x, ) can be approximated by conserving only the zero th order term in Eq. ( 22), and inserting it into Eq. ( 21). This leads to

I(x, ) = 1 2 1 + x [(1 + x) 2 + 2 ] 1/2 + 1 -x [(1 -x) 2 + 2 ] 1/2 [ψ](x, η). ( 23
)
Noting that the pressure jump must be zero at the side edges (i.e. [ψ](x, -1) = [ψ](x, 1) = 0), an integration by parts of Eq. ( 20) yields

w * (x, y) = 1 4πA - 1 -1 b η -y ∂ η [ψ](x, η)dη, ( 24
) with b = [(1 + x) 2 + A 2 (η -y) 2 ] 1/2 1 + x + [(1 -x) 2 + A 2 (η -y) 2 ] 1/2 1 -x , (25)
where the Cauchy principal value of the integral should be taken (indicated by the horizontal bar across the integral sign). The integral problem given by Eq. ( 24) could be inverted in the limit of small A. However, this would give a solution identical to the asymptotic solution presented in §2. Here, we use a different approach. We approximate this integral by first rewriting ∂ η [ψ](x, η) as The function b * (x) is chosen such that b × b * is approximately equal to 2, which is indeed true far from the leading edge and when (η -y) 2 is approximated by the value 1/4 (Fig. 3). Using this approximation (i.e. b × b * ≈ 2) and inserting it into Eq. ( 24), the integral problem is transformed into

∂ η [ψ](x, η) = b * ∂ η [ψ * ](x, η), with b * = 2/ 1 + [(1 -x) 2 + A 2 /4] 1/2 1 -x . ( 26
w * (x, y) 1 2πA - 1 -1 1 η -y ∂ η [ψ * ](x, η)dη. (27)
The above manipulations and approximations have led to an inverse problem where the kernel no longer depends on x and the unknown no longer depends on y. Equation ( 27) is a Carleman-type singular integral equation [START_REF] Carleman | Sur la résolution de certaines équations intégrales[END_REF] that can be inverted and whose solution ∂ η [ψ * ](x, η) is given by [START_REF] Muskhelishvili | Integrals of the Cauchy type[END_REF][START_REF] Tuck | Application and solution of Cauchy singular integral equations[END_REF])

∂ y [ψ * ](x, y) = - 2A π 1 -y 2 - 1 -1 1 -η 2 w * (x, η) η -y dη + C 0 1 -y 2 , ( 28
)
where the constant C 0 can be found by enforcing that the pressure jump is zero at the side edges. Integrating Eq. ( 28), using Eq. ( 26), and integrating with respect to y, the pressure jump [ψ] can be calculated. This approximate solution of the lifting-surface integral equation will be referred to as the approximate solution.

To give a clearer view of the solving procedure, we will consider a particular case: a flexible surface whose displacement only depends on x (i.e. w * is a function of x only). Using Eq. ( 28) and the boundary condition on the side edges, we obtain

∂ y [ψ](x, y) = 2Ab * w * (x) y 1 -y 2 . ( 29
)
And a simple integration gives the leading-order acceleration potential jump, which can be written in dimensional terms as

[P ](X, Y ) = 2ρb * D(W (X)) H 2 -Y 2 . ( 30
)
Now integrating Eq. ( 30) between -H and H gives the force per unit chord length exerted on the fluid

F a 0 (X) = b * ρπH 2 D(W (X)). ( 31
)
This equation is very similar to the results of Lighthill's slender-body theory, recalled in Eq. ( 10), except for the factor b * , which is plotted in Fig. 3 far from the trailing edge. At the trailing edge, b * is exactly zero, which is what is expected from the Kutta condition. In between, in a region whose chordwise length scales as the inverse of the span, the factor b * varies continuously between 1 and 0. In other words, the present approximate solution, given by Eq. ( 31), is equivalent to Lighthill's slender-body theory far from the trailing edge, but a correction is introduced near the trailing edge to fulfill the Kutta condition. We will now try to assess the benefits of this approach by comparing its results to the asymptotic solution introduced in section 2 and to the numerical calculations obtained through a panel method described in section 3.

Comparison between numerical and analytical results

Generic fluid-structure interaction problem

To assess the validity of the approximate solution, given in Eqs. ( 30) and ( 31), we consider a flexible plate with prescribed lateral motion. The leading edge of the surface is clamped at X = 0, while the trailing edge is free at X = 1 (Fig. 1). We choose the lateral displacement to be the first eigenmode of a clamped-free Euler-Bernoulli beam in vacuo [START_REF] Païdoussis | Fluid-structure interactions: Slender structures and axial flow[END_REF][START_REF] Païdoussis | Fluid-structure interactions: Slender structures and axial flow[END_REF]. The reduced frequency is defined as k = ΩL/U and, for simplicity, the pressure jump will be normalized by ρU 2 . The aspect ratio of the flexible plate is first set to A = 0.1. The distribution of real and imaginary parts of pressure jump derived respectively from the panel method (section 3) and the approximate solution (Eq. 30) are plotted in Fig. 4.

Focusing first only on the real parts of the pressure jumps (panels a and c in Fig. 4), we can see that the numerical and approximate pressure jump both exhibit an elliptical distribution in the spanwise direction, with zero pressure jump at the side edges (Y = 0 or Y = 0.1). A leading edge singularity is clearly visible in the numerical results (Fig. 4a), but not apparent in the approximate solution (Fig. 4c). Apart from this discrepancy near the leading edge (0 ≤ X ≤ 0.1), both approaches give similar results. The same conclusion can be reached for the imaginary parts of pressure jumps (panels b and d in Fig. 4). To reach a more quantitative comparison, we now compare the results of the panel method (section 3), slender-body theory [START_REF] Lighthill | Note on the swimming of slender fish[END_REF], the asymptotic solution (section 2), and the approximate solution (section 4), by plotting the chordwise force distribution (Fig. 5), now normalized by ρπH 2 U 2 /(2L) since it is a force per unit length. For this relatively small aspect ratio, all four approaches give similar results in the interval 0.1 ≤ X ≤ 0.9. Near the trailing edge, we see that the approximate solution gives the best result, with a force that is very close to the numerical results of the panel method (Fig. 5). In this region, slender-body theory is off, in particular because it does not fulfill the Kutta condition (i.e. F (X = 1) = 0), and the asymptotic solution gives a good result except very close to the trailing edge, where it diverges.

Near the leading edge, the asymptotic solution shows the best agreement with the numerical solution, while slender-body theory and the approximate solution give a similar solution, which was expected since b * ≈ 1 there (see discussion in section 4). The singularity at the leading edge results from a suction effect that could be probably be added in the model [START_REF] Lighthill | Aquatic animal propulsion of high hydromechanical efficiency[END_REF][START_REF] Wu | Hydromechanics of swimming propulsion. part 2. some optimum shape problems[END_REF][START_REF] Paraz | Thrust generation by a heaving flexible foil: Resonance, nonlinearities, and optimality[END_REF]. Note that this inverse-square-root singularity is integrable and does not play a major role in most applications since the leading edge is clamped and the aerodynamic loads applied near the trailing edge thus do negligible work. Note also that this singularity does not pose any problem for the numerical panel method thanks to its integrability.

To evaluate the pros and cons of the approximate solution for moderate aspect ratios, we now consider the same lifting-surface problem but with a larger aspect ratio: A = 0.3. The pressure jump distribution across the plate is plotted in Fig. 6, and the force distribution along the chord is plotted in Fig. 7. In Fig. 6, the comparison between the results of the panel method and the approximate solution leads to similar conclusions as for the smaller aspect ratio (Fig. 4). We find a good agreement except close to the leading edge. This is confirmed by comparing the force distributions in Fig. 7. The main difference with the results of aspect ratio A = 0.1 (Fig. 5) is that the leading-edge and trailing-edge regions are now wider, which is what is expected since the length of these regions scales as the plate span.

The average thrust generated by the deforming plate can also be used to compare the different methods. Following [START_REF] Lighthill | Note on the swimming of slender fish[END_REF], the thrust is evaluated by first calculating the mean rate of work W , as the work done by the force F (X) for a displacement of velocity iΩh(X) averaged over one period of time

W = Ω 2π 2π Ω 0 L -L iΩh(X)F (X)dXdT. ( 32
)
We then calculated the rate of kinetic energy lost into the wake as 1 2 ρπH 2 U W 2 (L). When this power lost in the wake is subtracted from W , we obtain the useful work of the thrust force T U , which means that the average thrust force T can be obtained through the equation where W is given by Eq. ( 32). As already noted by [START_REF] Lighthill | Note on the swimming of slender fish[END_REF], in slender-body theory, the thrust has to be calculated through this energy-based approach and cannot be calculated by directly projecting the force along the x-direction. We now assume that the lateral displacement of the plate is again the first mode of a clamped-free bream with a reduced frequency k = 10 (the value k = 0.5 used in Figs. 4 and5 could also be used but would yield a negative thrust). Using the formulation of the force given by Eqs. ( 10) and (31), for respectively the slender-body theory and the proposed approximate method, and inserting it into Eq.( 33), the average thrust can be calculated for any aspect ratio. These results are then compared to the numerical results obtained with the paned method in Fig. 8. Note that the asymptotic solution described in §2 cannot be used to calculate the thrust since the force has non-integrable singularities at the leading and trailing edges. It appears that the approximate solution give better results than slender-body theory. As expected, slender-body theory tends to overestimate the thrust generated because of the finite pressure jump at the trailing edge that contributes to the thrust generation. When Kutta condition is enforced, as it is the case in the approximate solution and in the numerical panel method, the pressure jump at the trailing edge is zero and the thrust is notably reduced.

T = W U - 1 2 ρπH 2 U W 2 (L), (33)
Figure 8: Comparison of the average thrust calculated with the panel method, slender-body theory, and the proposed approximate method. The displacement is the first eigenmode of a clamped-free beam, with reduced frequency k = 10. Note that the asymptotic method cannot be used to calculate the thrust as it diverges.

These results, and other tests that have been made with different reduced frequencies or different eigenmodes (not reproduced here), give credit to the approximate solution proposed in this paper. This approximate solution describes the aerodynamic loads better than Lighthill's slender-body theory, especially for moderate aspect ratios and near the trailing edge. Although the asymptotic solution may be more rigorous from a mathematical point of view, its numerical implementation is difficult because of the presence of finite-part integrals. In addition, the asymptotic solution diverges at the trailing edge, which may be a problem in some applications. This is why we believe the use of the approximate solution is preferable. Its implementing cost is minimal compared to slender-body theory and the Kutta condition is ensured at the c as a function of aspect ratio A for a clamped free plate of mass ratio M * = 1. Crosses are the results from [START_REF] Eloy | Flutter of a rectangular plate[END_REF]. The black curve is the marginal stability curve derived with slender-body theory and the red curve shows the result with the approximate solution.

trailing edge.

Application to the flag instability

To further evaluate the interest of the proposed approximate solution for fluid-structure interaction problems, we will now examine its application to the flag instability [START_REF] Païdoussis | Fluid-structure interactions: Slender structures and axial flow[END_REF][START_REF] Shelley | Flapping and bending bodies interacting with fluid flows[END_REF]. We consider again a clamped-free rectangular plate immersed in a potential flow of incoming velocity U (Fig. 1). Its semi-span is H, its semi-chord is L, its mass per unit area ρ p , and flexural rigidity D. We neglect viscoelastic damping and tension due to the skin friction or nonlinear effects.

Following [START_REF] Eloy | Flutter of a rectangular plate[END_REF], we use 2L, , 2L/U , ρU 2 the characteristic length, velocity, time, and pressure to make the problem dimensionless. The motion of flexible plate is then governed by the linearized Euler-Bernoulli beam equation

∂ 2 z ∂t 2 + 1 U * 2 ∂ 4 z ∂x 4 + M * F a 0 = 0, (34)
where U * is the reduced flow velocity and M * is the mass ratio,

U * = 2 ρ p D LU, M * = 2ρL ρ p , (35)
and F a 0 is the dimensionless force per unit length exerted by the plate on the fluid, which is given by the approximate solution (31).

This linear stability problem can be solved by decomposing the lateral deflection onto the eigenmodes of a clamped-free plate in vacuo. For a mass ratio M * , aspect ratio A immersed in a flow of reduced velocity U * , the plate is stable if the growth rates of all flutter modes are negative; otherwise it is unstable. The critical velocity U * c is the reduced flow velocity at which at least one flutter mode becomes unstable, when increasing U * from zero.

For a plate of mass ratio M * = 1, the critical flow velocity is plotted as a function of the aspect ratio A in Fig. 9. The results obtained using slender-body theory and the approximate solution are compared to the results of [START_REF] Eloy | Flutter of a rectangular plate[END_REF] obtained with a numerical spectral method. This comparison shows that slender-body theory fails for aspect ratios A 0.2, while the approximate solution gives accurate results up to A ≈ 0.7. This shows that the approximate solution, which can be viewed as a simple correction of Lighthill's slender-body theory, allows to extend the domain of validity up to moderate aspect ratios. Note that the discontinuities of U * c in Fig. 9 (when A = 0.54 for the black curve, A = 0.98 for the red one) correspond to shifts of the unstable flutter mode.

Conclusion

In this paper, we have proposed a method to improve Lighthill's slender-body theory for moderate aspect ratios. We started by recalling how the lifting-surface problem is stated and we gave an asymptotic correction of slender-body theory. We then briefly described a panel method in the frequency domain that has been used to obtain numerical results and compare them with the proposed analytical approaches. Finally, we proposed an approximate solution of the lifting-surface problem that can be viewed as correction of slender-body theory near the trailing edge.

To the benefits of the proposed approximate solution, we applied it to two problems: a generic lifting-surface problem and the flutter instability of a flexible plate (the flag instability). In the generic problem, we compared the pressure jump distribution due to an imposed periodic displacement with the four methods: the panel method, slender-body theory, the asymptotic solution, and the approximate solution. For lifting surfaces of small and moderate aspect ratios, the approximate solution gave superior results, except near the leading edge. In particular, the Kutta condition is fulfilled in the approximate solution, but not in the slender-body theory, nor in the asymptotic solution. In the flag instability problem, we calculated the critical flow velocity with slender-body theory and the approximate solution, and compared those results to those of a spectral method [START_REF] Eloy | Flutter of a rectangular plate[END_REF]. Again, the approximate solution appeared superior to slender-body theory for moderate aspect ratios.

The main contribution of this paper is the approximate solution of the lifting-surface problem given in Eqs. ( 30) and ( 31). This approximate solution can be viewed as a correction of Lighthill's slender body theory near the trailing edge that satisfies a Kutta condition. Although the proposed approximate solution is simple to implement and has no extra cost compared to Lighthill's slender-body theory, it appeared to give better results, in particular for moderate aspect ratios. We hope this approach will prove useful and will be applied in fluid structure interaction problems.

Figure 1 :

 1 Figure 1: Schematic representation of the lifting-surface problem. A flexible plate, with chord 2L and span 2H, is deforming periodically in an incoming flow U.

Figure 3 :

 3 Figure 3: The factor b * as a function of chordwise coordinate for various aspect ratios.

  as a function of x. It can readily be seen that b * is approximately 1

Figure 4 :

 4 Figure 4: Real and imaginary parts of numerical (a, b) and approximate (c, d) pressure jump distribution. The aspect ratio is A = 0.1, the reduced frequency is k = 0.5, and the displacement is the first eigenmode of a clamped-free beam in vacuo. The computation takes about 30 CPU minutes, with Nc = 900, Ns = 90, Nw = 1800 (convergence tests have shown that the error is less than 1%).

Figure 5 :

 5 Figure 5: Real (a) and imaginary (b) parts of the force per unit length. The parameters are the same as in Fig. 4.

Figure 6 :

 6 Figure 6: Same as Fig. 4 for an aspect ratio A = 0.3. The numerical parameters used are: Nc = 500, Ns = 150, Nw = 1000.

Figure 7 :

 7 Figure 7: Same as Fig. 5 for an aspect ratio A = 0.3.

Figure 9 :

 9 Figure 9: Critical flow velocity U *c as a function of aspect ratio A for a clamped free plate of mass ratio M * = 1. Crosses are the results from[START_REF] Eloy | Flutter of a rectangular plate[END_REF]. The black curve is the marginal stability curve derived with slender-body theory and the red curve shows the result with the approximate solution.

  . acknowledges support from the China Scholarship Council (No. 201204490206).