L Mónica Bravo-Anaya

Denis C D Roux
email: denis.roux@univ-grenoble-alpes.fr

J Félix

Armando Soltero Martí Nez

Francisco Carvajal Ramos

Frédéric Pignon
email: frederic.pignon@univ-grenoble-alpes.fr

Oonagh Mannix
email: oonmannix@gmail.com

Marguerite Rinaudo
email: marguerite.rinaudo@sfr.fr

Role of Electrostatic Interactions on Supramolecular Organization in Calf-Thymus DNA Solutions under Flow

Keywords: DNA, external ionic concentration, rheology, birefringence, SAXS, electrostatic interactions

Previous investigations were conducted on two concentrations of DNA solution: 4 mg/mL, for which it has been shown that no supramolecular organization is induced under flow at low shear rates; and 10 mg/mL, in which a liquid crystalline-type texture is formed under flow at low shear rates, attesting to an orientation of pre-organized chains. Rheological experiments are discussed and their results supported by small-angle X-ray scattering (SAXS) and flow birefringence visualization experiments. Scattering from polyelectrolytes has a characteristic signal, which is here observed in SAXS, showing a strong correlation peak between charged chains in water, for both concentrations. This peak is weaker in the presence of 0.01 M NaCl and suppressed in salt excess at 0.1 M NaCl. No plateau in the σ(𝛾 ´) plot was observed in analysis of rheological experiments on low DNA concentration (4 mg/mL). As typically observed in polyelectrolyte systems both the dynamic moduli and shear viscosity were higher in water as electrostatic forces dominate, than in the presence of salt, especially at low shear rates. The rheological results for concentrations of 0.01 M NaCl are lower than in water as expected due to partial screening of electrostatic repulsions.

Rheological data for concentrations of 0.1 M NaCl are unexpected. Electrostatic forces are partially screened in the low salt concentration, leading to a drop in the rheological values. For high salt concentration there are no longer interchain repulsions and so steric interactions dominate within the entangled network leading to the subsequent increase in rheological parameters. Regardless of the solvent, at high shear rates the solutions are birefringent. In the 10 mg/mL case, under flow, textures are formed at relatively low shear rate before all the chains align going to a pseudo-nematic liquid crystalline phase at high shear rate. The electrostatic repulsion between semi-rigid chains induces a correlation between the chains leading to an electrostatic pseudo-gel in water and loosely in 0.01M NaCl at low stress applied. To our knowledge, this is the first time that such behavior is observed. In 0.1M NaCl, DNA behavior resembles the corresponding neutral polymer as expected for polyelectrolyte in salt excess, exhibiting a yield stress. When texture appears in water and in 0.01M NaCl, a critical transition is observed in rheological curves, where the viscosity decreases sharply at a given critical shear stress corresponding to a plateau in the σ(𝛾 ´) plot also observed in creep transient experiment.

Introduction

DNA is a semi-rigid polyelectrolyte with a high persistence length (Lp~50 nm) [START_REF] Maret | Magnetic birefringence study of the electrostatic and intrinsic persistence length of DNA[END_REF][START_REF] Odijk | On the Ionic-Strength Dependence of the Intrinsic Viscosity of DNA[END_REF][START_REF] Manning | The Persistence Length of DNA Is Reached from the Persistence Length of Its Null Isomer through an Internal Electrostatic Stretching Force[END_REF].

Understanding electrostatic interactions between stiff chain segments in DNA, and polyelectrolyte solutions is important, both fundamentally and for biological applications [START_REF] Odijk | Possible scaling relations for semidilute polyelectrolyte Solutions[END_REF]. These interactions depend on the charge density of the polyelectrolyte favoring parallel alignment of the chains, and on the external salt concentration screening the electrostatic repulsions [START_REF] Malovikova | Viscosimetric behavior of Na-polygalacturonate in the presence of low salt content, Macro-ion characterization from dilute solutions to complex fluids[END_REF]. Due to condensed counterions in 1-1 electrolyte, there exists a range of distances less than the Debye length where polyions attract each other [START_REF] Ray | Formation of Loose Clusters in Polyelectrolyte Solutions[END_REF][START_REF] Ray | An Attractive Force between Two Rodlike Polyions Mediated by the Sharing of Condensed Counterions[END_REF][START_REF] Manning | Counterion condensation theory of attraction between like charges in the absence of multivalent counterions[END_REF][START_REF] Wissenburg | Multimolecular Aggregation in Concentrated Isotropic Solutions of Mononucleosomal DNA in 1 M Sodium Chloride[END_REF][START_REF] Wissenburg | Multimolecular Aggregation of Mononucleosomal DNA in Concentrated Isotropic Solutions[END_REF]. Birefringence has previously shown that textures corresponding to orientated, ordered domains are induced under shear flow for DNA concentrations higher than 5 mg/mL at low ionic concentration [START_REF] Bravo-Anaya | Supramolecular organization in calf-thymus DNA solutions under flow in dependence with DNA concentration[END_REF]. A schematic representation of the evolution of these solutions at rest and under flow was proposed [START_REF] Bravo-Anaya | Supramolecular organization in calf-thymus DNA solutions under flow in dependence with DNA concentration[END_REF]. To complete our previous work, the rheological behavior of DNA in aqueous solution is studied: at two polymer concentrations (4 mg/mL and 10 mg/mL); in water, and in the presence of two ionic concentrations: 0.01M and 0.1M NaCl; to investigate the influence of ionic concentration on rheological behavior and on birefringence induced under flow.

The polymer concentrations are chosen such that DNA is in the semi-dilute regime with a low degree of entanglement. SAXS experiments complete the data studying in the same solvent conditions, longrange interactions between chains.

Materials and Methods

Materials and solutions preparation

High-molecular weight calf-thymus DNA (MW=6 560 000 g/mol) was used to prepare solutions at concentrations of 4 and 10 mg/mL [START_REF] Bravo-Anaya | Supramolecular organization in calf-thymus DNA solutions under flow in dependence with DNA concentration[END_REF]. Ultrapure water and NaCl anhydrous were used to obtain different external ionic concentrations: 0, 0.01 and 0.1 M NaCl. The double helix conformation of DNA is stabilized under these conditions. Vials were closed and sealed with Parafilm ® to prevent water evaporation and concentration changes. Solutions were stored at 4 °C to prevent degradation.

Sigma-Aldrich Company supplied all reagents.

Rheometric measurements.

DNA rheology was studied at 25 ºC using a DHR-3 rheometer from TA Instruments Company and a MCR-501 rheometer from Anton-Paar Company. Two geometries were used: a steel cone-plate geometry with 40 mm diameter and an angle of 0.035 rad was used in DHR3 rheometer; and a glass plate-plate geometry with a 50 mm diameter was employed with the MCR501 rheometer to observe birefringence under crossed polarizers. Birefringence visualizations are integrated through the gap with a wave vector in the direction of the velocity gradient. A solvent cap was used to prevent solvent evaporation.

To define the linear viscoelastic regime oscillatory strain sweeps were undertaken at an angular frequency of 10 rad/s in the strain range from 0.01 to 100 %. Frequency sweeps were performed at 1% strain within the linear viscoelastic regime of G' and G'' in a frequency range ω from 0.01 to 100 rad/s, with at least 5 points per decade. Simple steady state shear was measured at shear rates from 10 -3 to 1 000 s -1 with at least 5 points per decade. Due to thixotropic behavior of 10mg/mL DNA solution, in order to reduce the time to obtain the steady state, increasing/decreasing ramps were imposed.

Firstly, a rapid shear rate increase from 0 to 100 s -1 was applied in 160 s followed directly by a decreasing ramp in 3000 s.

Transitory experiments in creep flow consist of applying constant stresses from 10 -2 to 20 Pa and observing the creep flow over a time up to 10 3 s for each applied stress. In the first part of the shear rate (time response), the transient regime is characterized by damped oscillations due to a coupling between tool inertia and solution viscoelasticity. Final constant shear rate corresponds to the steady state.

Polymers 2018, 10, x FOR PEER REVIEW 3 of 20

Flow birefringence visualizations

Birefringence was tested using a glass plate-plate configuration with an adapted gap, of at least 0.2 mm, such that the geometry was filled completely with DNA solution. This geometry is placed between crossed polarizers, one below and another just above the sample (Figure 1). The sample was illuminated vertically by a white source of light reflected by a mirror placed at 45º below the first polarizer under the geometry. A second mirror, placed at 45º above the upper polarizer, collected pictures at different imposed shear rates with CMOS device cameras (Canon 100 D or Sony DSLR-A580) and adapted commercial zooms. (1) white source, (2) lens, (3) and (8) mirrors, ( 4) and ( 7) crossed polarizers, ( 5) and ( 6) plate-plate geometry and (9) camera.

Small-angle X-ray scattering (SAXS) measurements

The SAXS experiments were performed on ID02 high brilliance beamline at European Synchrotron Radiation Facility (ESRF, Grenoble, France). DNA solutions were analyzed under controlled temperature (25 ± 1 °C) in a flow-through capillary cell (diameter ~ 2 mm). An incident monochromatized X-ray beam at wavelength of 0.995 Å was used. Scattered intensities were recorded on a two-dimensional CCD detector at sample-to-detector distances of 0.75 m and 8 m. The scattering intensity distribution as a function of scattering vector was obtained by radial integration of 2D scattering pattern. Background subtraction was done using a cell filled with water for all samples. This is not a perfect background for solutions with high salt concentration, as seen in the increased intensity in the high-q region. As scattering from the salt ions makes no contribution in the low-q region this does not change outcome of the analysis.

Results and discussion

X-Ray scattering

Both X-ray scattering and rheological experiments were performed in identical experimental conditions. The influence of ionic concentration is examined by adding external salt. A correlation peak is observed in the SAXS data, this is evidence of long-range electrostatic repulsions. These repulsions are caused by interactions between polymer chains forming a local network with a preferential inter-chain distance.

The correlation peak for polyelectrolyte solution in absence of external salt is located at q* such as 2 π/q*= 21nm at 10mg/mL DNA solution in water. It is progressively depressed in the presence of external salt as shown for 0.01M NaCl and is suppressed in 0.1M NaCl. Correlation peak extinction in salt excess is typically observed in polyelectrolyte solutions and also observed in neutron and light scattering experiments [START_REF] Nierlich | Small angle neutron scattering by semi-dilute solutions of polyelectrolyte[END_REF][START_REF] Morfin | Further evidence of liquid-like correlations in polyelectrolyte solutions[END_REF][START_REF] Borsali | Light scattering and small-angle neutron scattering from polyelectrolyte solutions: the succinoglycan[END_REF][START_REF] Milas | Influence of the shear rate on the small-angle neutron scattering pattern of polyelectrolyte solutions: the xanthan example[END_REF]. At high q values (q> 3 nm -1 ), the form factor of DNA is independent of the salt content (Figure 2) and agrees with the model given by Pringle and Schmidt [START_REF] Pringle | Small-angle X-ray scattering from helical macromolecules[END_REF]. is presented in the absence of external salt at 4 and 10 mg/mL. q*, corresponding to the maximum in the SAXS intensity as a function of q, decreases when the polymer concentration (Cp) decreases as found in the literature. The values of q* vary following a scaling law q*~Cp 1/2 in agreement with a hexagonal packing with a correlation distance between packed strands, indicating a preferential distance between chains. At 4mg/mL, q*=0.17nm -1 which gives a predicted value of q*=0.27nm -1 at 10 mg/mL using this scaling law. This value is in agreement with the experimentally determined value of q*= 0.28 nm -1 .

These results confirm our hypothesis on the formation of organized chains. 11 The SAXS values obtained indicate a nearly homogeneous distribution of stiff interacting chains throughout the solution. This agrees with the schematic representation given by Odijk for the case where polymer persistence length is greater than the distance between two chains [START_REF] Odijk | Possible scaling relations for semidilute polyelectrolyte Solutions[END_REF]. Figure 3. SAXS intensity I(q) as a function of the wave vector q for 4mg/mL and 10 mg/mL DNA concentrations at 25 °C.

To reflect the role of inter-chain electrostatic repulsion the DNA Debye-lengths were estimated. This includes the DNA counterion contribution (Cp (ion/L) = [(0.12/370)*C (g/L)] where 0.12 is the activity coefficient of counterions and 370 the average molar mass of the nucleic repeat unit; and the ionic contribution of external salt (Table 1). In 0.1M NaCl the behavior should approach the corresponding neutral polymer. In these experimental conditions, the influence of the ionic concentration, as demonstrated by the estimated Debye length, is low compared to the intrinsic value of the persistence length, around 50 nm. Finally, it is important to mention the large increase of I(q) at low q values indicating the formation of supramolecular structure as proposed in our model [START_REF] Bravo-Anaya | Supramolecular organization in calf-thymus DNA solutions under flow in dependence with DNA concentration[END_REF]. Moreover, the intensity profile at low q values (q~10 -2 nm -1 ) agrees with the existence of aggregates with constant dimension for H2O and 0.01M NaCl [START_REF] Potemkin | Associating polyelectrolytes: finite size cluster stabilization versus physical gel formation[END_REF]. In salt excess condition, I(q) slope is lower indicating a significant change of structure (Figure 2). Unfortunately, this does not allow direct determination of the persistence length from experimental data as obtained previously on another polyelectrolyte [START_REF] Moan | Study of the conformational rigidity of polyelectrolytes by elastic neutron scattering: 1. Carboxymethylcelluloses in the intermediate momentum range[END_REF].

Birefringence and texture under shear

Data previously obtained for DNA solution at 10 mg/mL, in TE buffer at ionic concentration around 0.01M, were completed with birefringence observations, as a function of shear rate. The two polymer concentrations are chosen as at 4 mg/mL the absence of birefringent texture was previously demonstrated, while texture was observed at 10 mg/mL. This enables investigation of influence of texture on rheological behavior.

In Table 2 it is shown that salt has a critical role in observing birefringence, and birefringent textures. Figures 4 and5 present examples of birefringence and the textures observed.

Table 2. Influence of polymer and ionic concentration on the critical shear rate for birefringence and texture formation. The role of polymer concentration on the supramolecular organization of DNA in solution is demonstrated for the 10 mg/mL DNA concentration. Birefringence occurs at low shear rates, but larger shear stress due to the high solution viscosity. At 10mg/mL, liquid crystalline like textures occur over shear rates of 1s -1 confirming our previously published data [START_REF] Bravo-Anaya | Supramolecular organization in calf-thymus DNA solutions under flow in dependence with DNA concentration[END_REF]. In the presence of salt at 0.1M NaCl, heterogeneous birefringent elongated domains appear under shear. The density of these domains progressively increases to form an almost homogeneous nematic phase. In all these experiments, a slightly negative normal force is observed when chains are orientated. It increases when the textures appear and become slightly positive.

Rheological behavior

To complete our previous rheological experiments on DNA solutions [START_REF] Bravo-Anaya | Supramolecular organization in calf-thymus DNA solutions under flow in dependence with DNA concentration[END_REF], we examined the rheological behavior at different external salt contents: in the absence of salt, and at 0.01M and 0.1M

NaCl. Electrostatic repulsions are maximum in the absence of salt and progressively screened with the addition of salt. At 0.1 M in which the Debye length is small (Table 1), the behavior goes to that of the corresponding neutral polymer involving steric hindrance between semi-rigid chains. The influence of external salt, on viscometric measurements, is well described in literature [START_REF] Malovikova | Viscosimetric behavior of Na-polygalacturonate in the presence of low salt content, Macro-ion characterization from dilute solutions to complex fluids[END_REF]. Two different polymer concentrations are tested to evidence clearly the influence of texture on the rheological behavior in the semi-dilute regime, 4 mg/mL being the model for absence of texture.

Behavior at 4 mg/mL

This solution is considered as a model for an entangled semi-dilute system in which the role of salt concentration, intrinsic stiffness, and absence of "liquid crystalline" ordered domains or texture are investigated as a function of the ionic concentration. 3. It is seen that the bulk relaxation time is longer in water, as is the stiffness of the system, due to stronger inter-chain repulsion requiring a stronger stress to perturb the supramolecular structure composed of correlated chains. The G'-G'' crossing transition ω0 does not increase linearly with salt concentration. Firstly, ω0 increases upon addition of salt (0.01M NaCl). This demonstrates a decrease in electrostatic forces, further evidenced by the reduction of the correlation peak in SAXS. Then ω0 decreases upon further addition of salt (0.1M NaCl). This state resembles that of the corresponding uncharged polymer with stronger entangled contacts than in the case of simple electrostatic repulsion (seen for the 0.01 M NaCl concentration) leading to an increase in G'. 7. The well-known behavior of polyelectrolyte viscosity is observed in water. In absence of external electrolyte, the viscosity increases rapidly when the shear rate decreases due to strong electroviscous effects [START_REF] Yamanaka | Revisit to the intrinsic viscosity-molecular weight relationship of ionic polymers. 2. Viscosity behavior of salt-free aqueous solutions of sodium poly (styrenesulfonates)[END_REF]. In 0.01M NaCl, the screening gives a pseudo plateau at low shear rate when the ionic contribution of DNA is in the same range as the solvent ionic concentration as discussed previously [START_REF] Malovikova | Viscosimetric behavior of Na-polygalacturonate in the presence of low salt content, Macro-ion characterization from dilute solutions to complex fluids[END_REF].

 (s -1 )
. In excess of salt at 0.1M NaCl, the DNA solution behaves as the corresponding neutral polymer with a viscous plateau below 0.02 s -1 . The viscosity is lower than DNA solution in water but larger than in 0.01M NaCl due to screening of electrostatic long-range interactions. This is directly related to the dynamic behavior discussed previously (Figure 6). Over 𝛾 ´ > 1 s -1 , the viscosities are nearly the same in both NaCl-containing solvents corresponding to progressive flow alignment. Birefringence measurements shows a progressive increase interpreted as an alignment of domains culminating in the alignment of all chains forming a pseudo-nematic phase (see pictures in Figure 4). To demonstrate the validity of the dynamic and shear viscosities, the Cox Merz representation is given in Figure 7. s -1 , with flow viscosities being lower than the complex viscosities. This is related to the screening of electrostatic repulsions and subsequent increase in the elastic contribution of the chains.

Observing a decrease in viscosity in polyelectrolyte solutions is a normal behavior upon increasing salt concentration due to a decrease in the electroviscous effects. However, upon complete electrostatic screening, as in the case of 0.1 M NaCl, the viscosity increases at low shear rate compared with DNA solution in 0.01M NaCl. This is due to the large steric interactions present in the semidilute entangled solution.

In Figure 9, the viscosity is plotted as a function of the applied shear stress. At low shear stresses viscosity of the 0.1M NaCl solution is lower than that of water but higher than the 0.01M NaCl solution. A plateau is observed when NaCl is present but not in water alone. In water the viscosity increases rapidly at very low shear stress this is related to the existence of a yield stress caused by an electrostatic network formation. This point will be discussed later for higher polymer concentration. corresponding to the appearance of birefringence, and for stress greater than 4 Pa complete orientation of chains and a pseudo-nematic phase is observed. Figure 9 shows a smooth transition in η (𝛾 ´) curve confirming the absence of a plateau on 𝜎(𝛾 ˙)representation as shown in our previous work.

Behavior at 10 mg/mL

Under the same conditions as for DNA 4mg/mL solution, dynamic experiments were performed in the linear domain. The important difference concerns the lower stress at which birefringence appears, and the presence of texture in water and low ionic concentration (Table 2, Figure 5).

Dynamic measurements

At low frequencies the behavior, as observed by dynamic rheological measurements, is nearly identical in the three solvents. An original behavior is observed in Figure 10.

Flow measurements

It is difficult to reproduce viscosity measurements at low shear rate for concentrated systems [START_REF] Barnes | A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure[END_REF][START_REF] Barnes | The yield stress-a review or 'παντα ρει'-everything flows?[END_REF][START_REF] Auffret | Stress and strain controlled rheometry on a concentrated lyotropic lamellar phase of AOT/Water/Iso-octane[END_REF][START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF][START_REF] Cheddadi | Steady Couette flows of elastoviscoplastic fluids are nonunique[END_REF]. This is illustrated, in Figure 11 in 0.01M NaCl. Firstly, an increasing shear rate is imposed resulting in a stress versus shear rate dependence with a slope of 1 up to around 0.1 s -1 , corresponding to an apparent Newtonian plateau. Above 0.1s -1 , shear stress follows a smooth evolution up to 100 s -1 corresponding to a bulk flow regime. On the decreasing shear rate ramp, the shear stress decreases taking values larger than those obtained during the increasing ramp at low shear rate. The increase behavior is characteristic of slipping at the walls and thixotropy occurring in more concentrated solutions [START_REF] Barnes | A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure[END_REF][START_REF] Meeker | Slip and flow in pastes of soft particles: direct observation and rheology[END_REF]. In concentrated system, those generic flow curves are the signature of an apparent yield stress fluid with slipping as confirmed hereinafter. In fact, the yield stress is proportional to the elastic shear modulus [START_REF] Barnes | A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: its cause, character, and cure[END_REF], which is nearly the same in the three solvents with a small influence of electrostatic interactions due to the presence of external salt (Figure 10). . From these data, the complex viscosity fairly agrees with the shear viscosity measurements following the Cox Merz rule (Figures 10 and12) for which G' and G'' overlap over ω=10 -1 rad/s. At low frequency, in the range of gel-like behavior, the Cox Merz rule is no more valid in presence of external salt probably due to slip at the wall in absence of any supramolecular structural change as indicated in Table 2.

To relate the rheological behavior with supramolecular organization, viscosity is plotted as a function of the applied stress in Figure 14. It is shown that the stress needed to orientate the electrostatic network with flow direction in water is around 10 Pa much larger than for 4 mg/mL DNA solution (around 2 Pa, Figure 9). In these curves, two series of values are plotted and are found in good agreement: values from flow experiments on the decreasing step and values from the stationary state from creep step (see later in Figures 161718). In the presence of 0.01M NaCl, the orientation stress increases up to almost 15 Pa. A higher energy is required to align the system due to partial screening of the electrostatic interactions (Figure 14). This is in agreement with the lower order parameter as shown in X-rays and the appearance of a plateau in the σ(𝛾 ´) representation (Figure 15). In the presence of 0.1M NaCl, a smooth variation of the viscosity as a function of stress is obtained and is in good agreement with absence of a plateau in σ(𝛾 ´) representation (Figure 15). .

(10 Pa) at which the viscosity drops (Figure 14) giving a plateau in the σ(𝛾 ´) curve (Figure 15). In 0.01M

NaCl DNA solution a transition occurs around 5 Pa. This corresponds to the appearance of an equivalent texture (Figure 5 b). The transition is completed at a constant shear stress (12 Pa) before the viscosity drops. This drop in viscosity corresponds to a plateau in the σ(𝛾 ´) curve. The vertical viscosity drop is longer in water than in 0.01M NaCl. This is due to the decrease of electrostatic interchain interactions in the presence of external salt. In 0.1MNaCl DNA solution, the viscosity decreases progressively as a function of stress applied (Figure 14). Unlike in water and 0.01M NaCl solutions, no texture is observed (Figure 5c). Above shear stress of 8 Pa, an heterogenous birefringent solution organization across the entire volume is obtained indicating the orientation of large associated domains going to fully birefringent solution at high shear stress. Irrespective of the solvent, for shear stresses greater than 15-20 Pa, complete orientation of DNA forms a pseudonematic phase as suggested (birefringence is maximum) [START_REF] Bravo-Anaya | Supramolecular organization in calf-thymus DNA solutions under flow in dependence with DNA concentration[END_REF].

Transient regime

The transitory state is defined when constant shear stress is applied. Creep experiments enable this regime to be studied. In water, at 10 -2 Pa, no flow occurs (Figure 16). After transient inertia viscoelastic coupling, the shear rate decreases until 10 -4 s -1 for 10 3 s, indicating a solid like behavior in water (as indicated in Figure 10). The sample starts flowing at 10 -1 Pa corresponding to 3 × 10 -3 s -1 .

Flow is accompanied by a change in the polymer structure when the network is ruptured as evidenced by a rapid decrease in viscosity (Figure 12). The first oscillations in Figure 16 occur at a characteristic time, which is progressively shifted for stress applied higher than 5 Pa. At 7 Pa the oscillations disappear giving a bump around 10 s corresponding to the transition in shear stress (Figure 14). Over 9 Pa, corresponding to the pseudo plateau in shear stress versus shear rate representation, when the texture is observed the supramolecular structure is stabilized very rapidly. the sample behaves as a solid (Figure 17). The sample begins to flow at a lower yield stress than the sample in water and is faster to reach the stationary state. For samples in 0.01 M NaCl solution flow occurs at circa 1 Pa, and the stationary state is obtained after a maximum of 100 s. When stress above 10 Pa is applied the sample exhibits a textured structure under birefringence. Physically this behavior is caused by a looser electrostatic network providing less hindrance to movement than for DNA in water solution corresponding to lower viscosity (Figure 12). In 0.1M NaCl, the flow occurs at very low shear stress corresponding to an apparent yield stress lower than 0.1 Pa (Figure 18). In this salt concentration, the supramolecular organization changes to a looser structure.

Irrespective of solvent conditions, the DNA solutions at 10 mg/mL present a yield stress at very low shear stress in agreement with dynamical and flow measurements (Figures 10 and12).

Conclusions

The influence of the external salt concentration on the rheology and supramolecular structure of DNA solutions was discussed in this paper. Rheological measurements were performed with SAXS and birefringence experiments. In water, strong interchain electrostatic repulsions induced the existence of a correlation peak between chains in solution. This correlation peak is progressively depressed in presence of external salt. Two DNA concentrations in the semi entangled regime were chosen: 4mg/mL which exhibits birefringence at high shear rates, but no particular texture is observed; and 10 mg/mL where clear supramolecular structure was evidenced previously in 0.01M buffer solution [START_REF] Bravo-Anaya | Supramolecular organization in calf-thymus DNA solutions under flow in dependence with DNA concentration[END_REF]. This work clearly describes the influence of polymer concentration and external ionic concentration on the rheological behavior and supramolecular structure of this original stiff DNA polymer. It is demonstrated that at higher polymer concentration, a plateau is obtained in the σ(𝛾 ´) representation as soon as long-range electrostatic interactions and texture exist. An important point is that the rheology is dependent on the history of the samples at higher polymer concentration (as seen in the example of DNA solution at 10mg/mL). The experimental conditions need to be precisely controlled, especially the ramp speed in viscosity and shear rate experiments, to obtain robust data.

For the first time, the presence of an electrostatic network is demonstrated, especially in water at 10 mg/mL DNA concentration where electrostatic repulsions are at a maximum. In addition, for all samples studied it is shown that the rheological behavior is particularly original at low flow shear rate and low dynamic frequency applied.

Figure 1 .

 1 Figure 1. Schematic description of the home-made birefringence set-up on the rheometer MCR-501.

Figure 2 .

 2 Figure 2. SAXS intensity I(q) as a function of the wave vector q at different external salt concentrations. DNA solution at 10 mg/mL at 25°C.

Figure 4 .

 4 Figure 4. Birefringence of DNA solution at 4 mg/mL in 0.1M NaCl under shear a) at 10 -3 s -1 b) at 500 s -1 .

Figure 5 .

 5 Figure 5. Supramolecular organization of 10 mg/mL DNA solution sheared at 25°C: a) Texture obtained in water; b) in 0.01M NaCl; c) heterogeneous birefringence in 0.1M NaCl. The white bars represent the scale of 10 mm. Two videos for solutions in water and in 0.1M NaCl under shear are joined in SI.

11

 11

3. 3

 3 .1.1. Dynamic measurements Rheological experiments were performed in the linear domain. The influence of external salt concentration (Cs) is shown in Figure 6 and summarized in Table

Figure 6 .

 6 Figure 6. Elastic (G') and viscous (G") moduli as a function of the angular frequency at 25°C for different solvent ionic concentrations 0, 0.01M and 0.1M NaCl for ADN 4 mg/mL.

Figure 7 .

 7 Figure 7. Flow viscosities of 4 mg/mL DNA solutions as a function of the shear rate at 25°C solvents with different ionic concentrations: 0, 0.01M and 0.1M NaCl.

Figure 8 .

 8 Figure 8. Cox Merz overlay for 4mg/mL DNA solution at 25°C in different ionic concentrations: 0, 0.01M and 0.1M NaCl. Continuous lines represent the flow viscosities and the dots represent dynamic oscillatory viscosity measurements.

Figure 8

 8 Figure 8 shows the complex viscosity derived from dynamic experiments and viscosity from flow behavior are in agreement, for the frequency range covered here, corresponding to the Cox Merz rule. In the presence of external salt, there is dissociation of the two types of viscosity curve over 1

Figure 9 .

 9 Figure 9. Viscosity as a function of the applied shear stress on 4mg/mL DNA solutions in different solvents: 0, 0.01M and 0.1M NaCl at T=25°C.

Figure 10 .

 10 Figure 10. Elastic and viscous moduli obtained by dynamic oscillatory rheology for 10 mg/mL DNA solution at different ionic concentrations: 0, 0.01M and 0.1M NaCl (T=25°C).

Figure 11 .

 11 Figure 11. Shear stress as a function of the shear rate obtained for an increasing followed by a decreasing shear rate ramp on DNA solution in 0.01M NaCl at T=25°C.

Figure 12

 12 Figure 12 is a plot of viscosity against shear rate. The points were obtained while decreasing the shear rate. It can be seen that the values are lower for solutions with added salt. This is expected behavior for polyelectrolytes. In all three cases, viscosity increases at low shear rate in agreement with the existence of a yield stress. Over 10 -1 s -1 the viscosities in the three solvents are overlapping in agreement with the behavior in dynamic rheological measurements. The flow viscosity and complex viscosity from the oscillatory behavior are shown in Figure 13.

Figure 12 .

 12 Figure 12. Flow viscosities as a function of the shear rate for DNA 10 mg/mL solutions in the 3 solvents: H2O, 0.01M and 0.1M NaCl. Values are taken on the decreasing shear rate ramp.



  (rad/s),  (s -1 )

Figure 13 .

 13 Figure 13. Cox Merz overlay for 10 mg/mL DNA solution at 25°C in different ionic concentrations: 0, 0.01M and 0.1M NaCl. Continuous lines represent the flow viscosities and the dots represent dynamic oscillatory viscosities measurements.

Figure 14 .

 14 Figure 14. Viscosity as a function of the shear stress at 25°C in H2O, 0.01M NaCl and 0.1M NaCl. Measures were performed on decreasing shear rate sweep (open symbols) and in the stationary state from creep step (filled symbols).

Figure 15 .

 15 Figure 15. Shear stress as a function of the shear rate for DNA 10 mg/mL solution at 25°C in H2O, 0.01M NaCl and 0.1M NaCl.

Figure 16 .

 16 Figure 16. Transient behavior for DNA solution 10mg/mL in water: creep experiment giving the shear rate as a function of time for each stress applied from 0.1 to 10 Pa.

1 )Figure 17 .

 117 Figure 17. Transient behavior for DNA solution 10mg/mL in water: creep experiment giving the shear rate as a function of time for each stress applied from 0.1 to 10 Pa.

Figure 18 .

 18 Figure 18. Transient behavior for DNA solution 10mg/mL) in 0.1M NaCl: creep experiment giving the shear rate as a function of time for each stress applied from 0.1 to 14 Pa.

At 4

 4 mg/mL in water, the G' modulus and the viscosity under flow are higher than in presence of salt. Viscosity decreases steeply around a stress of 5 Pa and birefringence occurs over 50 s-1 corresponding to alignment of all the chains in the flow. Salt in solution changes the rheological behavior due to screening of electrostatic repulsions. The local organization is depressed but it is shown that the dynamic modulus and flow viscosity are lower in 0.01M NaCl than in 0.1M NaCl. In this case, 0.1 M NaCl concentration, long range electrostatic interactions are absent and steric interactions dominate in the entangled semi-rigid chain assembly, inducing higher G' and viscosity in the Newtonian plateau.Irrespective of ionic concentration, increasing the stress progressively allows orienting the interacting chains and the curves obtained for the different solvents superimpose over a stress around 2 Pa (Figure9). Larger stresses correspond to the entrance in the birefringent regime where all the chains orient in the flow over 4 Pa. Clearly, from Figure9, a smooth transition is obtained confirming the absence of plateau on σ(𝛾) representation. At 10 mg/mL, an electrostatic network with G'(ω) = G"(ω) is demonstrated from dynamic rheology at frequency lower than 0.08 rad/s in water; electrostatic interactions decrease progressively in presence of external salt being progressively substituted by steric interactions. Then, an apparent yield stress is found from stress measurement at low shear rate and confirmed in transient experiments in water and to a lower degree in NaCl. Due to this mechanism, viscosity measurements as a function of the shear rate are delicate and the values must be taken on the decreasing step in shear rate. Due to this gelation mechanism, observation time dependence involved in experiments under stress implies change in the supramolecular structure and dynamic of the systems. The curves giving the viscosity as a function of the shear stress applied presents a two steps transition: i) first one around 4 Pa (over 1 s -1 ) where texture starts; the second one at 10 Pa in water; ii) at 4 Pa and around 12 Pa in 0.01M NaCl respectively corresponding to texture appearance with large decrease in viscosity at constant shear stress. The formation of the texture corresponds to a plateau in the shear stress (shear rate) representation. Textures in 0.01M NaCl are less visible than in water corresponding to a lower organization of stiff chains. Nevertheless, in 0.1M NaCl, smooth variation of the viscosity with shear stress with looser organization (no texture as shown in water and 0.01M NaCl) is obtained with no plateau in the stress versus shear rate representation. In 0.1M NaCl, heterogeneous birefringent solution formed by large domains orientated in the flow. Whatever the ionic conditions, at higher shear stress, all the chains orient in the flow with a low viscosity and the solutions are birefringent corresponding a pseudo-nematic supramolecular organization of the system.

Table 1 .

 1 Estimation of the Debye screening length κ -1 expressed in nm. .

	Polymer concentration (mg/mL)	κ -1 water (nm)	κ -1 0.01M NaCl (nm)	κ -1 0.1M NaCl (nm)
	4	8.42	2.83	0.94
	10	5.36	2.61	0.93

Table 3 .

 3 Characterization of the elastic (G') and loss (G'') moduli at the transition for different ionicViscosities as a function of the shear rate are plotted in Figure

	concentrations.
		H2O	0.01M NaCl	0.1M NaCl
	G0 (Pa)	0.84	0.68	0.78
	ω0 (rad/s)	0.31	0.82	0.59
	η0 (Pa.s) for 𝛾 ˙→ 0 𝑠 -1	18 at 𝛾 ´=10 -3 s -1	4.1	7.2
	3.3.1.2. Flow measurements

  Supplementary Materials:The following are available online at www.mdpi.com/link, Video S1, S2: DNA solutions 10 mg/mL in water and in 0.1M NaCl under shear are joined in SI.

	527	conditions	of	the	Creative	Commons	Attribution	(CC	BY)	license
	(http://creativecommons.org/licenses/by/4.0/).

Acknowledgments:

We acknowledge the European Synchrotron Radiation Facility for providing synchrotron radiation facilities and we would like to thank Theyencheri Narayanan and Sylvain Prévost for assistance in using beamline ID02. L.M. Bravo-Anaya acknowledges the fellowship grant given by CONACYT (CVU 350759). Laboratoire Rheéologie et Procédés is a part of the LabEx Tec21 (Investissements d'Avenir-Grant agreement no. ANR-11-LABX-0030), of PolyNat Carnot Institut (Investissements d'Avenir-Grant agreement no. ANR-11-CARN-030-01) and of the program ANR-15-IDEX-02. We acknowledge M. Karrouch and E. Faivre from Laboratoire Rhéologie et Procédés, Grenoble, for their valuable help and technical support.


¸ois Bouchon
email: francois.bouchon@uca.fr

Thierry Dubois
email: thierry.dubois@uca.fr

Nicolas James
email: nicolas.james@math.univ-poitiers.fr

AN IMMERSED METHOD BASED ON CUT-CELLS FOR THE SIMULATION OF 2D INCOMPRESSIBLE FLUID FLOWS PAST SOLID STRUCTURES

Keywords: Immersed boundary methods, Cutt-cell methods, Incompressible viscous flows

We present a cut-cell method for the simulation of 2D incompressible flows past obstacles. It consists in using the MAC scheme on cartesian grids and imposing Dirchlet boundary conditions for the velocity field on the boundary of solid structures following the Shortley-Weller formulation. In order to ensure local conservation properties, viscous and convecting terms are discretized in a finite volume way. The scheme is second order implicit in time for the linear part, the linear systems are solved by the use of the capacitance matrix method for non-moving obstacles. Numerical results of flows around an impulsively started circular cylinder are presented which confirm the efficiency of the method, for Reynolds numbers 1000 and 3000. An example of flows around a moving rigid body at Reynolds number 800 is also shown, a solver using the PETSc-Library has been prefered in this context to solve the linear systems.

Introduction

For some decades, many researchers and engineers have been considering the numerical solution of fluid flows, for different kind of fluids and different geometries. With the increasing performance of super-computers, it has been possible to tackle more and more challenging problems, for higher Reynolds numbers and complex geometries. Several different discretization techniques can be used to consider these problems: Finite Element Methods, Finite Volumes Methods, Spectral Methods, and Finite Difference Methods. The MAC scheme on cartesian grids [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible fl ow of fluid with free surface[END_REF] can be viewed both as a Finite Volume Method or a Finite Difference Method on staggered grids, and is adapted to 2D or 3D flows in simple geometries, for example for lid-driven cavity or backward facing step. To take into account some obstacles in the flows (or complex geometries), immersed boundary techniques have been developped by Peskin in the 80's ( [START_REF] Peskin | The fluid dynamics of heart valves: experimental, theoretical, and computational methods[END_REF][START_REF] Peskin | The immersed boundary method[END_REF]), consisting in using Dirac functions to model the interacting force between the fluid and the solid structure. These methods have inspired many authors in the following years, Mohd-Yusof has combined them with the use of B-Splines ( [START_REF] Mohd-Yusof | Combined immersed-boundary/B-Spline methods for simulations of flow in complex geometries[END_REF]) in his momentum forcing methods to consider complex geometries. The main advantage of these techniques is that the forcing term does not change the spatial operators, making them quite easy to implement (see [START_REF] Mittal | Immersed boundary methods[END_REF] for a review, and refences therein). As an alternative, Bruno et al. have developped penalization techniques to inforce suitable boundary conditions [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF]. Similar techniques have also been investigated by Maury et al. [START_REF] Janela | A penalty method for the simulation of fluid -rigid body interaction[END_REF] and justified from a mathematical point of view in [START_REF] Maury | Numerical analysis of a finite element/volume penalty method[END_REF]. These methods have been shown to be efficient in the context of several particles in a flow [START_REF] Lefebvre | Fluid-particle simulations with freefem++[END_REF], and when considering possible collisions between them [START_REF] Verdon | Contact problems for particles in a shear flow[END_REF].

Arbitrary Lagrangian Eulerian (ALE) methods have been developped for flows in geometries which vary in time (see [START_REF] Richter | A fully eulerian formulation for fluid-structure-interaction problems[END_REF][START_REF] Richter | A monolithic geometric multigrid solver for fluid-structure interactions in ale formulation[END_REF][START_REF] Yang | An ale approach to mechano-chemical processes in fluid-structure interactions[END_REF] where authors use some of the ideas of [START_REF] Belytschko | Fluid-structure interaction[END_REF]). The aim is to formulate the equation in a fixed reference domain, by using a mapping from the reference domain Ω(0) to the domain Ω(t) occupied by the fluid at time t. The position of the moving bodies, which correspond to the boundary of Ω(t), being available, the velocity field of these bodies defined on ∂Ω(t) have to be extended to Ω(t). Once this is done (generally with harmonic extensions), the equations are written in the reference domain by using the chain-rule formula.

For problems involving non-rigid bodies, Roshchenko et al. (see [START_REF] Roshchenko | A time splitting fictitious domain algorithm for fluidstructure interaction problems[END_REF]) have used splitting methods to solve first the evolution of the velocity field in the fluid, and then to consider the deformation of the body. These ideas of splitting the model can be viewed as similar to the projection techniques (see [START_REF] Chorin | Numerical solution of the navier-stokes equations[END_REF], and [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] for a review and references therein).

The method presented in this work joints another family of methods, called cut-cell methods. The idea of these methods is to modify the discretization of the Navier-Stokes Equations in the cells cut by the immersed boundary (see [START_REF] Ye | Numerical simulation of two-dimensional flows over a circular cylinder using the immersed boundary method[END_REF][START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF][START_REF] Tucker | A cartesian cut-cell method for incompressible viscous flow[END_REF][START_REF] Chung | Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape[END_REF]). One can discretize the equations on smallest cells obtained by intersecting the grid-cell with the domain occupied by the fluid, or one can merge these smallest cells with a neighbouring one. These methods can be combined with the levelset methods to track the boundary of the fixed or moving body (see [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF]). These ideas are used in the present work: the body in the fluid is represented by a levelset function, and the location of the velocity components are modified in the cut-cells (see [START_REF] Bouchon | A second-order cut-cell method for the numerical simulation of 2d flows past obstacles[END_REF]), the pressure remaining placed at the center of cartesian grid cells for both fluid-cells and cut-cells. For the Laplacian of the velocity, the classical five-point approximation must be replaced by a local 6-point formula, for which the truncation error is only first order. But as in [START_REF] Matsunaga | Superconvergence of the shortley-weller approximation for dirichlet problems[END_REF], global second order convergence of the method is recovered. This second-order convergence for the velocity and the pressure with our cut-cell scheme has been obtained for the flow past a circular cylinder at Re=40 in [START_REF] James | Assessment of immersed boundary method toward high accuracy[END_REF] by comparing with the reference solution proposed in [START_REF] Gautier | A reference solution of the flow over a circular cylinder at re = 40[END_REF].

The paper is organized as follows: Section 2 is devoted to the presentation of the problem. The Navier-Stokes Equations are considered in a 2D geometry, which is supposed to be fixed for the sake of clarity. We also introduce there the notation for the grids, and detail space discretization. In Section 3, we give some information about computational aspects. In the case of fixed domain, a fast solver adapted from the capacitance matrix method (see [START_REF] Buzbee | The direct solution of the discrete poisson equation on irregular regions[END_REF][START_REF] Buzbee | The direct solution of the biharmonic equation on rectangular regions and the poisson equation on irregular regions[END_REF]) is used to solve the linear systems for both components and for the pressure. We also show that the method can be adapted in the case of moving domain. In this context, the preprocessing step of the capacitance matrix method would need to be done at each time-iteration which would then increase the CPU time. Therefore, we have prefered for this case a parallel version of an algebraic multigrid method (HYPRE BoomerAMG) implemented using the PETSc Fortran library (see [START_REF] Balay | PETSc Web page[END_REF][START_REF] Balay | PETSc users manual[END_REF]).

Section 4 is then devoted to numerical tests, where we compare our results with theoretical predictions and other numerical results available in the literature. A conclusion is given in section 5. 2. The Setting of the Problem 2.1. Flows past obstacles. We consider a flow in a two-dimensional domain Ω = (0, L)×(0, H) which contains a domain Ω S occupied by the solid which is supposed to be fixed for the sake of simplicity. We denote then Ω F = Ω \ Ω S the domain occupied by the fluid (see Fig. 1).

The velocity field in the fluid satisfies the Navier-Stokes equations, with no-slip boundary conditions. We consider then the problem:

∂u ∂t - 1 Re ∆u + ∇ • (u ⊗ u) + ∇p = 0, (1)
∇ • u = 0, (2) u(x, t = 0) = u 0 , (3)
where u(x, t) = (u, v) is the velocity field at x = (x, y) ∈ Ω F at time t > 0, u 0 is the initial condition and Re is the Reynolds number. We impose homogeneous Dirichlet boundary conditions for the velocity field on ∂Ω F :

u = 0 on ∂Ω F (4)
We mention that non-homogenous Dirichlet boundary conditions can also be treated with the method presented here. 2.2. Discretization. For the time-discretization of (1)-(3), we use a second-order backward difference (BDF2) projection scheme. In a first step, the velocity field is advanced in time with a semi-implicit scheme decoupling the velocity and pressure unknowns. Then, the intermediate velocity is projected in order to obtain a free-divergence velocity field.

Let δt > 0 stand for the time step and t k = k δt discrete time values. Let us consider that (u j , P j ) are known for j ≤ k. The computation of (u k+1 , P k+1 ) needs two steps:

(5) 3 u k+1 -4u k + u k-1 2δt - 1 Re ∆ u k+1 + ∇P k = -2 ∇ • (u k ⊗ u k ) + ∇ • (u k-1 ⊗ u k-1 )
with homogeneous Dirichlet boundary condition for u k+1 . Then the intermediate velocity field u k+1 is projected in the free-divergence space to get u k+1 :

u k+1 -u k+1 δt + 2 3 ∇(P k+1 -P k ) = 0, ∇ • u k+1 = 0, (u k+1 -u k+1 ) • n = 0 on Γ. (6)
For the spatial discretization, we modify the MAC scheme near the boundary by changing the location of the unknowns of the velocity components for the cells cut by the solid as depicted on Fig. 2 , the pressure unknowns remaining in their original place (see [START_REF] Bouchon | A second-order cut-cell method for the numerical simulation of 2d flows past obstacles[END_REF] for more details). To discretize the Laplacian in (5), we must replace the five-point formula by a six-point discretization. For the convective terms, the fluxes are computed at the midle of the vertical and horizontal edges (see Fig. 3). For the pressure, linear interpolation are used rather than chang- Although the truncation error is only first order in space for the resulting numerical scheme, the second order accuracy is recovered which is due to a superconvergence phenomena analoguous to those proven in [START_REF] Matsunaga | Superconvergence of the shortley-weller approximation for dirichlet problems[END_REF]. This second order has been observed in [START_REF] James | Assessment of immersed boundary method toward high accuracy[END_REF] by showing results in comparison with those of [START_REF] Gautier | A reference solution of the flow over a circular cylinder at re = 40[END_REF].

Computational Aspects

3.1.

A fast parallel direct solver to treat fixed solid strutures. When considering fixed solid bodies, the use of a direct solver with a preprocessing procedure is efficient. Once the preprocessing computations have been done, the cost paid to solve the linear systems is about twice the case of a numerical simulation in the same computational domain without obstacles. We summarize hereafter the fast direct solver derived from the capacitance matrix method and adapted to the case of non uniform grids (see [START_REF] Bouchon | A second-order cut-cell method for the numerical simulation of 2d flows past obstacles[END_REF] and [START_REF] Bouchon | The immersed interface technique for parabolic problems with mixed boundary conditions[END_REF] for the details) which has been implemented in our code.

After spatial discretization of the Navier-Stokes equations, one linear equation is obtained per node in the part of the computational domain filled by the fluid and per unknown that is u, v and p. We complete these sets of linear equations by adding similar ones for nodes of the cartesian grid lying inside the solid obstacle but with zero as right-hand side. The unknowns corresponding to mesh points in Ω S h are fictitious ones. As in Ω F , the numerical scheme accounts for the boundary conditions on Γ S h , the fluid unknowns are independant to the solid ones. We therefore obtain linear algebraic systems defined on the whole cartesian grid with n x × n y mesh points whose sizes are (n x -1) × n y for u, n x × (n y -1) for v and (n x -1) × (n y -1) for p. All three linear systems are similar in nature : the resulting matrices have similar structures with five or six non-zero coefficients per row.

Let us consider one of these linear system. We denote by N its size and by A ∈ M N (R) its matrix. Then at each time iteration, we have to solve a linear system [START_REF] Bouchon | A second-order cut-cell method for the numerical simulation of 2d flows past obstacles[END_REF] AX = Z with the right-hand side Z computed from the velocity and the pressure at previous time steps.

As it is mentioned above, the matrix A is non-symmetric. Let us consider now the matrix G obtained with the same discretization on the whole computational domain Ω totally filled by a fluid that is with no obstacles. The matrices A and G differ only on rows corresponding to computational meshes for which the five-point stencil interacts with a cut-cell. Let us denote by n c this total number of rows, namely rows such that A -G have non-vanishing coefficients. The efficiency of our direct solver is due to the fact that n c is small compared with N and that the non-zero coefficients on each row of A -G is bounded. The linear system (7) can be rewritten as

(8) GX = Z -QY
where Q is a matrix of dimensions N × n c with one non-vanishing coefficient per column, equal to one, and

Y ∈ R nc such that Q Y = (A -G) X. It can be easily shown, using Q t Q = I nc ,
that Y is solution of the following linear system (9)

I nc + M G -1 Q Y = M G -1 Z with M = Q t (A -G). The matrix I nc + M G -1
Q is a non-singular matrix (see [START_REF] Bouchon | The immersed interface technique for parabolic problems with mixed boundary conditions[END_REF] for a proof) of size n c . Based on these relations, the algorithm implemented to solve [START_REF] Bouchon | A second-order cut-cell method for the numerical simulation of 2d flows past obstacles[END_REF] consists in a preprocessing step where the matrix I nc + M G -1 Q is factorized (we use a LU -factorization) followed by i) Compute Z and solve GW = Z ; ii) Compute M W and solve (9) ; iii) Compute QY and solve GX = Z -QY . Recalling that G is the matrix corresponding to the standard MAC scheme on the whole computational mesh, steps i) and iii) can be performed by using any efficient solvers available on cartesian grids. In the present work, we use Discrete Fourier transforms in the vertical direction (where the mesh is uniform) combined with LU -factorizations of the resulting tridiagonal systems.

The parallel version of this direct solver is based on explicit communications performed by calling functions of the MPI library. The main feature of MPI is that a parallel application consists in running p independant processes which may be executed on different computers, processors or cores. These processes can exchange datas by sending/receiving messages via a network connecting all the involved computing units.

The first step when developping a parallel algorithm is to define a suitable and efficient splitting of the datas among the MPI processes : each MPI process will treat datas associated with a part of the total computational mesh. For our problem, this choice is straightforward and is related to the algorithm used to solve the linear systems. Indeed, it is much easier to implement a parallel resolution of tridiagonal linear systems rather than a parallel version of the DFT. Therefore, the parallel version of the code is based on a splitting of the datas along the horizontal axis, so that each MPI process works with a vertical slice of the computational mesh as it is illustrated on Figure 4.

In the framework of finite volume or finite difference schemes on cartesian grids, the explicit computation of spatial derivatives is local and involves very few communications. The only tricky part concerns the resolution of the linear systems. The step ii) of the direct solver described in the previous section consists in solving a linear system involving the matrix

I nc + M G -1 Q.
As the LU -factorization of this matrix has been computed and stored in a pre-processing step at the beginning of the time iterations, we have to solve two triangular systems which can not be efficiently performed on parallel computers. As n c is small compared to the size of the global problem, we choose to dedicate this task to one given MPI process (fixed in advance) per unknown, that is u, v and p. Once these linear systems are solved, the resulting vectors are scattered from these MPI processes to the other ones.

P k-1 P k P k+1
As it was previously mentioned, linear systems of steps i) and ii) are solved by first applying a DFT in the y-direction : these computations are independant and can be performed without any communications due to the distribution of datas among the MPI processes. This results in a collection, one per grid point in the y-direction, of independant tridiagonal linear systems connecting all nodes of the mesh in the x-direction. A parallel direct solver based on the divide and conquer approach (DAC) for tridiagonal matrices has been implemented (see [START_REF] Bondeli | Divide and conquer: a parallel algorithm for the solution of a tridiagonal linear system of equations[END_REF]). The DAC method, applied to solve one tridiagonal linear system on n p > 1 MPI processes, consists in splitting the tridiagonal matrix into n p independant blocks (one per MPI process). The solutions of these systems have to be corrected in order to recover the solution of the global system. These corrections correspond to 2n p -1 values which are solutions of a tridiagonal linear system of size 2n p -1. This phase of the DAC method is sequential and has to be performed on one process inducing a useless waiting time for the other processes. However, as we have to solve n y such systems simultaneously, this sequential part can be distributed among all the n p processes. In this context, the DAC algorithm leads to an efficient parallel code.

The parallel code has a good level of performance : less than 15% of the CPU time is spent in communications between MPI processes. The sequential part performed on one process represents a negligible amount of CPU time. The computations presented here have been performed on a DELL cluster using up to 32 cores of Xeon processors. A low latency bandwith network connects the cluster nodes.

3.2.

Iterative solver for the case of moving bodies. For solid bodies moving in a computational domain filled by a fluid, as the case considered in Section 4.2, cut-cells may change at each time iteration. Therefore, the coefficients of the matrices of the linear systems for the velocity components, issued from the discretizetion of the momentum equations, and for the pressure increment computed in the projection step of the time scheme, have to be recomputed at each time step. In that context, the use of the fast direct solver described in the previous section is cumbersome and inefficient except on coarse meshes. In order to be able to treat such configurations, we have implemented a PETSc version [START_REF] Balay | PETSc Web page[END_REF][START_REF] Balay | PETSc users manual[END_REF] of our cut-cell scheme. The main advantage of the PETSc Library is that, in a parallel programming environment based on MPI, many iterative solvers combined with different preconditionners can be used. The choice can be made at run time.

Numerical Results

4.1.

Flow past a circular cylinder at Re = 1000 and 3000. In this section, we present numerical simulations, performed with the parallel direct solver method described in Section 3.1. We consider the case of flows past a fixed circular cylinder of diameter D. The Reynolds number is defined based on the diameter D of the cylinder, i.e. Re = U ∞ D/ν where U ∞ is the horizontal free stream velocity. As non-dimensional time, we consider T = 2U ∞ t/D.

A circular cylinder of diameter equals to unity is centered at the origin of the computational domain Ω = (-L x , L x ) × (-L y , L y ). As boundary conditions, a uniform velocity profile u(x = -L x , t) = (1, 0) is imposed at the inflow and a convective boundary condition is applied at the exit, namely the convective equation ( 10)

∂u ∂t + (1, 0) • ∇u = 0 is solved at x = L x .
On the top and bottom boundaries, that is y = ±L y , slip boundary conditions are used that is ∂u ∂n = 0 and v = 0. Finally, no-slip (u| Γ S = 0) boundary condition is applied on the surface of the obstacle.

For this problem important quantities reflecting the dynamics of the vortices formed in the vicinity of the solid boundary and developping at the rear of the cylinder are the pressure drag and lift coefficients. They are derived from the total drag force on the body, which is computed as As whown in Figure 5, the values obtained with the cut-cell scheme on a grid with 4096 × 8192 mesh points discretizing the domain Ω = (-10, 10) 2 perfectly match the theoretical curve drawing [START_REF] Chung | Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape[END_REF] on the time interval T ∈ [0, 0.2] for Re = 1000. Our cut-cell method captures the square-root singularity of the drag coefficient. This simulation has been run using 16 MPI processes. On longer time interval, namely T ∈ [0, 5], the results are in good agreement with those obtained by Koumoutsakos and Leonard in [START_REF] Koumoutsakos | High resolution simulations of the flow around and impulsively started cylinder using vortex methods[END_REF] with a vortex method. In order to test the grid convergence of these results, the same simulation has been conducted on a grid with two times more points in both spatial directions, that is 8192 × 16384 mesh points, in the same computational domain. In that case, 32 MPI processes have been used. Both results are almost indistinguishable on Figure 6 indicating that the coarser resolution is enough to capture the essential features of the flow at this Reynolds number. The development of the flow around the impulsively started cylinder at Re = 1000 can be seen on Figure 7. In the early stage T ≤ 1, a primary vortex develops in the vinicity of the boundary at the rear of the cylinder. Then for T ∈ [1, 2] a secondary vortex appears trying to move insight the primary vortex and to push it away from the solid boundary (T ≥ 4). A tertiary vortex is visible at T = 3 which remains sticked to the boundary constrained by the two other vortices having more strength. These results compare well with the same flow representations shown in [START_REF] Koumoutsakos | High resolution simulations of the flow around and impulsively started cylinder using vortex methods[END_REF]. As expected on short time interval (T ≤ 5) the flow remains symmetric. At Re = 3000, the time evolution of the drag coefficient plotted on Figure 8 exhibits also the square-root singularity on short time interval and is in good agreement with the results of Koumoutsakos and Leonard. As expected, the drag coefficients remains almost constant for T ∈ [START_REF] Balay | PETSc Web page[END_REF][START_REF] Balay | PETSc users manual[END_REF] (see [START_REF] Koumoutsakos | High resolution simulations of the flow around and impulsively started cylinder using vortex methods[END_REF]). Note that on this time interval, a small difference exists between the results computed on the two different grids. However, the coarser simulation is fine enough to capture the flow dynamics at Re = 3000. The mesh size of the coarser grid, which is constant Red solid line: numerical results on a 4096 × 8192 grid in Ω = (-10, 10) 2 ; Green solid line: numerical results on a 8192 × 16384 grid in Ω = (-10, 10) 2 ; Blue dots: results from [START_REF] Koumoutsakos | High resolution simulations of the flow around and impulsively started cylinder using vortex methods[END_REF].

in the vicinity of the solid boundary, is h = 20/8192 ≈ 2.44 × 10 -3 . Note that the size of the computational domain and the boundary conditions imposed at the exit may influence the results. Estimating the values for L x and L y required so that the numerical results being of the order of the numerical scheme error, namely O(h 2 ), is an open question. This will be addressed in further works. At this Reynolds number, a scenario similar than that at Re = 1000 can be observed on Figure 9 with the development of three vortices in the early stage of the flow dynamics. The secondary vortex penetrates further inside the primary vortex aera and the tertiary vortex has more strength as it could be expected with less effects of the viscous forces at Re = 3000. Again, an overall good agreement is found with the vortcity contours shown in [START_REF] Koumoutsakos | High resolution simulations of the flow around and impulsively started cylinder using vortex methods[END_REF] at the same Reynolds number and time T . Red solid line: numerical results on a 4096 × 8192 grid in Ω = (-10, 10) 2 ; Green solid line: numerical results on a 8192 × 16384 grid in Ω = (-10, 10) 2 ; Blue dots: results from [START_REF] Koumoutsakos | High resolution simulations of the flow around and impulsively started cylinder using vortex methods[END_REF].

T = 0 T = 1 T = 2 T = 3 T = 4 T = 5 T = 6 T = 7
As previously mentioned, the flow remains symmetric at the beginning of the simulations for these flows around an impulsively started cylinder. By carrying the time integration over a much longer time interval T ∈ [0, 200] instabilities due to round-off errors and to the nonlinearity of the system develop so that the flow becomes non symmetric for T ≥ 100 at Re = 1000 and T ≥ 50 at Re = 3000 as it can be seen on Figures 10 and11 representing the time history of the drag coefficient. After a transient period, an increase of the drag coeffient is observed which stabilizes and oscillates around a mean value.

4.2.

Flow around moving bodies. The purpose of this section is to show that the present numerical method is also able to simulate incompressible flows around moving bodies. Let us consider a cylinder which starts to move impulsively at t = 0 with the sinusoidal translational motion [START_REF] Gautier | A reference solution of the flow over a circular cylinder at re = 40[END_REF] u body (t) = 2 sin t 2 ; 0 in a fluid initially at rest for Reynolds number 800. We suppose that the fluid is confined within a rectangular computational domain Ω = [-3; 3]×[-1; 1] with no-slip boundary condition on ∂Ω F . The diameter of the cylinder is equal to 1 and it is initially centered at the origin. The boundary condition (13) at the body surface ∂Ω S is enforced through the non exhaustive following right hand side terms which vanish in the case of fixed obstacle : u(κ u,S i,j ), u(ξ S i,j-1/2 , y j-1/2 ), u body (κ S i,j ) and so on. This terms are respectively taking into account in the convective terms, Laplacian operator and continuity equation. More details can be found in [START_REF] Bouchon | A second-order cut-cell method for the numerical simulation of 2d flows past obstacles[END_REF].

As the obstacle moves from one time step to another, we have to update matrices of the linear systems corresponding to Poisson and momentum equations at each iteration. Therefore, in such configuration, the direct solver requires much more CPU time compared to some iterative solver, which does not require a preprocessing step. For large problems, the faster solver we have found is an algebraic multigrid method (HYPRE BoomerAMG) implemented using the PETSc library [START_REF] Balay | PETSc Web page[END_REF][START_REF] Balay | PETSc users manual[END_REF].

A constant mesh size h = 5 × 10 -3 is used in both directions and the value of the time step, satisfying a CFL stability condition, is 10 -3 . As shown in Figure 12, vortices interact

Conclusion

We have presented a cut-cell method for the numerical solution of flows past obstacles. We have detailed the numerical method and the computational aspects for fixed obstacles, and shown numerical results for fixed and moving rigid bodies. The parallel version of the algorithm presented here allows computations of flows at Reynolds number up to 3000. The numerical tests confirm some results of the literature, a good agreement is observed with the numerical simulations in [START_REF] Koumoutsakos | High resolution simulations of the flow around and impulsively started cylinder using vortex methods[END_REF]. We have also shown that the computation of the drag coefficient matches the theoretical square-root singularity predicted by [START_REF] Bar-Lev | Numerical solution of the navier-stokes equations[END_REF]. The choice of the size of the box (compared with the grid size h) is one of the questions that we would like to investigate with this method in further works, we also would like to deal with rigid bodies following the fluid flow.

Figure 1 .

 1 Figure 1. The solid body Ω S with boundary Γ S and the surrounding computational domain Ω F in which the flow is to be simulated.

Figure 2 .Figure 3 .

 23 Figure 2. Location of the unknowns near the solid body

Figure 4 .

 4 Figure 4. Splitting of the computational grid among the MPI processes P k , k = 0, . . . , p -1. The gray zone refers to additional (ghost points) storage used for the MPI communications between neighbooring processes.

  The pressure drag and lift coefficients C p and C are given by C p = 2F b • e x and C = 2F b • e y and the (total) drag coefficient isC d = C p + C .Starting with a flow at rest, the drag coefficient behaves as T -1/2 in the early stage of the development of the vortices. This squareroot singularity has been theoretically predicted by Bar-Lev and Yang in[START_REF] Bar-Lev | Numerical solution of the navier-stokes equations[END_REF]. They have derived the following expression for the (total) drag coefficient

Figure 5 .Figure 6 .

 56 Figure 5. Evolution of the drag coefficient of a circular cylinder at Re = 1000. Solid line (red): theoretical prediction (12); Blue dots : numerical results on a 4096 × 8192 grid in Ω = (-10, 10) 2 .

Figure 7 .

 7 Figure 7. Vorticity of a flow past a circular cylinder at Re = 1000 simulated on a grid with 4096 × 8192 mesh points in the computational domain Ω = (-10, 10) 2 at different times T ∈ [0, 7].

Figure 8 .

 8 Figure 8. Evolution of the drag coefficient of a circular cylinder at Re = 3000. Red solid line: numerical results on a 4096 × 8192 grid in Ω = (-10, 10) 2 ; Green solid line: numerical results on a 8192 × 16384 grid in Ω = (-10, 10) 2 ; Blue dots: results from[START_REF] Koumoutsakos | High resolution simulations of the flow around and impulsively started cylinder using vortex methods[END_REF].

Figure 10 .Figure 11 .

 1011 Figure 10. Evolution of the drag coefficient of a circular cylinder for Re = 1000 simulated on a grid with 4096 × 8192 mesh points in the computational domain Ω = (-10, 10) 2 .

Figure 12 .

 12 Figure 12. Flow around a moving circular cylinder at Re = 800 in the computational domain Ω = (-3, 3) × (-1, 1) discretized with 1200 × 400 mesh points.

Acknowledgements

This research was partially supported by the French Government Laboratory of Excellence initiative n o ANR-10-LABX-0006, the Région Auvergne and the European Regional Development Fund. The numerical simulations have been performed on a DELL cluster with 32 processors Xeon E2650v2 (8 cores), 1 To of total memory and an infiniband (FDR 56Gb/s) connecting network.

with each other and also with the boundaries. The flow remains perfectly symmetric until t = 6π, thereafter the symmetry of the flow is lost due to rounding errors inherent in computer calculations (see Figure 13).