Emilie Ressouche

Sandrine Pensec

Benjamin Isare

Jacques Jestin

Laurent Bouteiller

Two-component self-assemblies: investigation of a synergy between bis-urea stickers

It is of interest to develop two-component systems for added flexibility in the design of supramolecular polymers, nanofibers or organogels. Bis-ureas are known to selfassemble by hydrogen bonding into long supramolecular objects. We show here that mixing aromatic bis-ureas with slightly different structures can yield surprisingly large synergistic effects. A strong increase in viscosity is observed when a bis-urea with the sterically demanding 2,4,6-trimethylbenzene spacer is combined with a bis-urea bearing no methyl group in position 2 of the aromatic spacer (i.e. 4-methylbenzene or 4,6-dimethylbenzene). This effect is the consequence of a change in supramolecular assembly triggered by the composition of the mixture. The mixture of complementary bis-ureas forms rod-like objects that are more stable by about 1 kJ/mol, and that are thicker than the rod-like objects formed by both parent systems.

 and the aim of the present manuscript is to characterize in detail the structure and the stability of these co-assemblies.

Experimental Section

Synthesis. The synthesis of T, 29 M, 27 X, 30 T PIB , 31 and X PIB 28 was previously reported. The synthesis and the characterization of M PIB are reported in Supporting Information.

Introduction

Viscosimetry. Measurements were performed with an Anton Paar AMVn falling ball microviscosimeter with a 0.16 mm Ø capillary, with 3 measurements at an angle of +20° and -20°. Results are reported as an average of those 6 measurements.

Small-angle neutron scattering. Small-angle neutron scattering measurements were made at the LLB (Saclay, France) on the Pace instrument, at two distance-wavelength combinations to cover the 4 10 -3 to 0.24 Å -1 q-range, where the scattering vector q is defined as usual, assuming elastic scattering (q=(4π/λ)sin(θ/2), where θ is the angle between incident and scattered beam).

Data were corrected for the empty cell signal and the solute and solvent incoherent background.

A light water standard was used to normalize the scattered intensities to cm -1 units.

The data was fitted with the DANSE software SasView. The simplest form factor that allowed to fit all the mixtures at both temperatures was the form factor of a cylinder with an elliptical cross-section and a homogeneous contrast. The scattering length densities for toluene-D 8 (5.66 10 -6 Å -2 ) and for the low molar mass bis-urea M (6.78 10 -7 Å -2 ) were calculated from the atomic bound coherent scattering lengths. The other samples contained the polymeric bis-urea X PIB , so that the cylinder can be considered to contain a significant amount of solvent. Therefore, the scattering length density for all the mixtures containing X PIB was adjusted during the simultaneous fit of all the data. The common value that allowed to fit all the systems was 2.00 10 -6 Å -2 . Three parameters were adjusted for each sample: the length of the cylinder and the minor and major radii of the elliptical cross-section. Differential scanning calorimetry. Solution-phase differential scanning calorimetry measurements were performed on a high sensitivity TA Instruments nDSC III system (baseline noise ± 15 nW). The sample cell (0.3mL) was filled with the bis-urea solution and the reference cell with the corresponding solvent. The capillary cells were not capped, and a constant pressure of 5 10 5 Pa was applied. The solutions (2 g/L), were analyzed between 10 and 80 °C, using 3 full heating/cooling cycles, at 1 °C min -1 . The transition temperature, T**, was measured from the second heating scan.

Fourier

Results and Discussion

Comparison of the stickers. Bis-urea T (Figure 1) was shown to self-assemble by hydrogen bonding into two distinct supramolecular structures of high molar masses. In toluene both structures are in competition, and it is possible to switch the assembly from one structure to the other by changing the temperature. [START_REF] Bouteiller | Thickness Transition of a Rigid Supramolecular Polymer[END_REF] The high temperature structure is a long filament with a single bis-urea in the cross-section, [START_REF] Brocorens | Conformational Plasticity of Hydrogen Bonded Bis-Urea Supramolecular Polymers[END_REF] that is responsible for the formation of viscous solutions.

The low-temperature structure is a very long and rigid tube with three bis-ureas in the crosssection, that forms viscoelastic solutions (see Figure 2a for a schematic representation). [START_REF] Shikata | Structure and Dynamics of a Bisurea-Based Supramolecular Polymer in N-Dodecane[END_REF] Bisureas X and M self-assemble into the same supramolecular structures, but their relative stabilities are strongly affected by the substitution on the aromatic spacer. The presence of the methyl groups in positions 4 and 6 stabilizes the tube form of X because the repulsion between the methyl and carbonyl groups enforces a noncoplanar conformation of the urea and phenyl moieties and better preorganizes the monomer for hydrogen bonding in the tube structure. [START_REF] Isare | Conformational Control of Hydrogen-Bonded Aromatic Bis-Ureas[END_REF] In contrast, the additional presence of a methyl group in position 2 destabilizes the tube form of M because this methyl group is placed within the tube cavity, where it competes with the entry of solvent molecules. [START_REF] Isare | Engineering the Cavity of Self-Assembled Dynamic Nanotubes[END_REF] As a consequence, at 10 mM in toluene, the tube structure is stable up to 43°C for bis-urea with tolyl spacer (T), [START_REF] Bouteiller | Thickness Transition of a Rigid Supramolecular Polymer[END_REF] up to 110°C for bis-urea with xylyl spacer (X) [START_REF] Isare | Conformational Control of Hydrogen-Bonded Aromatic Bis-Ureas[END_REF] and only below -5°C for bis-urea with mesityl spacer (M), [START_REF] Isare | Engineering the Cavity of Self-Assembled Dynamic Nanotubes[END_REF] so that at room temperature both T and X form viscoelastic solutions but not M (Figure 1). Mixing bis-ureas T and M significantly improves the stability of the viscoelastic solution (up to 73°C at 10 mM in toluene), [START_REF] Isare | Engineering the Cavity of Self-Assembled Dynamic Nanotubes[END_REF] which shows that the tubes with a mixed composition are considerably more stable than any of the pure tubes. This synergistic effect was attributed to the enhanced preorganization of bis-urea M and to the wider cavity created by bis-urea T: addition of bis-urea M to the T/M mixture is favorable only if its content is low enough, so that the tubular cavity is still large enough to accommodate the solvent molecules. Building on these observations, we first decided to investigate if similar or even stronger synergistic effects could be obtained with the other combinations. Figure 1 shows the appearance of solutions of T, M and X and their equimolar mixtures: although qualitative, this test shows that the composition of the solutions has a strong effect on their flow behavior. The solutions flow in the following order: M < T < T/X << T/M << X << M/X, which means that both T/M and M/X mixtures flow more slowly than their pure components, but not the T/X mixture. This indicates that the presence of the bis-urea with the mesityl spacer (M) is essential to obtain a clear-cut synergy, which is in line with the previous explanation that the methyl group in position 2 of the aromatic spacer allows to stabilize the tube assembly, but only if mixed with a less sterically demanding bis-urea. Among the various systems tested, the M/X mixture forms a particularly strong gel that is worth investigating in detail. However, bis-urea X forms metastable

H N H N H N H O O N H N H N H N H O O 22 22 N H N H N H N H O O N H N H N H N H O O N H N H N H N H O O 18 18
X PIB N H N H N H N H O O 16

M X PIB

solutions that yield a fine precipitate after several weeks, [START_REF]Bis-ureas T and M form homogeneous solutions in toluene by simply stirring at room temperature, and these solutions are stable for unlimited periods of time. In contrast, bis-urea X only dissolves at high temperature and remains optically clear at room temperature for a few weeks only[END_REF] which complicates structural and thermodynamic studies. Therefore, we considered the bis-ureas made of the same stickers linked to small polyisobutene (PIB) side-chains in the hope that the PIB chains would improve the solubility in non-polar solvents: bis-ureas T PIB , M PIB and X PIB (Figure 1) are indeed perfectly soluble in toluene. Table S2 summarizes the appearance of the mixtures combining any low molar mass bis-urea (T, M or X) to any polymeric bis-urea (T PIB , M PIB or X PIB ). This table shows that M/T PIB and M/X PIB are systems of interest since they form gels, unlike their pure components. Some of the other mixtures may also be of interest, but the effect may be hidden by the fact that bis-ureas T and X already form a gel in these conditions.

Therefore, infrared spectroscopy was used to discriminate more sensitively the systems of interest, because the shape of the N-H vibration band is known to be a characteristic signature for the filament and tube structures. [START_REF] Bouteiller | Thickness Transition of a Rigid Supramolecular Polymer[END_REF] Figure 3 shows the relevant part of the spectra for the M/X PIB mixture. The three spectra show two bands at 3330 and 3270cm -1 that are due to hydrogen bonded N-H urea groups (aliphatic N-H and aromatic N-H, respectively). The broadness of these bands for M and X PIB is characteristic for the filament structure, whereas the narrower bands in the case of their mixture is characteristic for the tube structure. [START_REF] Isare | Engineering the Cavity of Self-Assembled Dynamic Nanotubes[END_REF][START_REF] Bouteiller | Thickness Transition of a Rigid Supramolecular Polymer[END_REF] Therefore, the nature of the assemblies can be estimated from the ratio of absorbances between the maximum (at 3330cm -1 ) and the minimum (at 3300cm -1 ). From the data on Figure 3, it is clear that the M/X PIB and M/T PIB mixtures form synergistic co-assemblies because their A max /A min ratio is larger than for both pure components. The effect is less obvious, but may also be present in the case of T/X PIB , T/M PIB , X/T PIB , X/X PIB and X/M PIB mixtures for which the A max /A min ratio is larger than the average values of the pure components, meaning that the component that forms the tube structure on its own (T or X), converts at least part of the other component (T PIB , X PIB or M PIB ) from the filament to the tube structure. characteristic for the filament (resp. tube) structure. [START_REF] Bouteiller | Thickness Transition of a Rigid Supramolecular Polymer[END_REF] See Table S3 for the tabulated values.

While all the previous tests were performed with equimolar amounts of both components, it is of interest to see if unbalanced mixtures can also display some synergistic co-assembly. Figure 4 shows capillary viscosity measurements of T PIB , M PIB or X PIB solutions at a fixed concentration )"

ν"(cm 51 )" ""XPIB" M/XPIB" """"M"

X PIB% M/X PIB% M % A max" A min" no#addi've# T# M# X# 1# 1,1# 1,2# 1,3# 1,4# TPIB# MPIB# XPIB# no#addi've# A max /A min# X PIB% M PIB% T PIB% no#addi've # X % M % T %
no#addi've # (40 g/L) in which a low amount (less than 10 mol%) of T, M or X was added. In these conditions, only two systems display a strong increase in viscosity (M/X PIB and X/X PIB ), and the former system is the most interesting because the viscosity of the mixture is much larger than the viscosity of both pure components. Therefore, this system was selected for further characterization.

1" 10" 0" 0,

η/η 0"

[T]"(mM)"

T/TPIB" T/MPIB" T/XPIB" T" Structural and thermodynamic study of the M/X PIB mixture. First of all, the temperature stability of the supramolecular structures was probed by a combination of VT-FTIR spectroscopy and DSC. Figure 5 shows the typical broadening of the FTIR spectra at high temperatures that is characteristic of the transition from the tube to the filament structure (see above). The transition temperature (T** = 70°C) was measured as the mid-point in the plot of the A max /A min ratio.

T/T PIB% % T/M PIB% " T/X PIB% % T" 1" 10" 0" 0,5" 1" 1,5" 2" η/η 0" [M]"(mM)" *" M/TPIB" M/MPIB" M/XPIB" M" M/T PIB& & M/M PIB& " M/X PIB& & M" 1" 10" 0" 0,5" 1" 1,5" 2" η/η 0" [X]"(mM)" *" X/TPIB" X/MPIB" X/XPIB" X" X/T PIB& & X/M PIB& " X/X PIB&
Figure 6 shows the DSC trace obtained for the same system: the transition is detected through an endothermal peak (T** = 68°C, ΔH** = 6.4 kJ/mol). A significant hysteresis is noticed between the heating and cooling curves, but this transition is perfectly reversible, with heating curves that superimpose over at least four heating-cooling cycles. This calorimetric data was obtained for various compositions and is plotted on Figure 6: it shows that the addition of X PIB to M stabilizes the tube structure by about ΔT** = 75°C. 36 Moreover, the curve has a very broad maximum, which indicates that the optimal composition is close to 1/1, but that deviation from this stoichiometry is not critical. From this calorimetric data, it is also possible to estimate the free energy gain of the M/X PIB tube of 1/1 composition compared to the pure M tubes through equation 1. 37

∆∆𝐺 = ∆𝐻 * * ∆! * * ! * * ~1 𝑘𝐽/𝑚𝑜𝑙 (1)
This value can be compared to the free energy for the formation of such hydrogen bonded supramolecular polymers: in the case of bis-urea T, the association constant for the chain elongation (K) was measured by calorimetry in toluene at 25°C. [START_REF] Isare | Conformational Control of Hydrogen-Bonded Aromatic Bis-Ureas[END_REF] This association constant is related to the free energy of tube growth through equation 2.

∆𝐺 = -RTlnK = -34 𝑘𝐽/𝑚𝑜𝑙 (2)
Therefore, the stabilization induced by the complementarity between stickers represents only a small fraction of the overall system. The evolution of the transition temperature versus composition (Figure 6) is actually confirmed at the macroscopic level by viscosity measurements (Figure 7) that show a strong increase of viscosity for compositions centered around 1/1. The higher the temperature, the narrower the composition range where a gel is formed, and at 80°C (i.e. above the maximum of T**) no increase in viscosity is detected for the mixtures. In order to directly evaluate the enthalpic gain of combining M and X PIB molecules, we attempted an Isothermal Titration Calorimetry (ITC) experiment where a solution of M in toluene is injected into a solution of X PIB at the same concentration. The concentration of the experiment had to be lowered to 0.05 mM to reduce the baseline noise due to the viscosity of the solution. Figure 8 shows the corresponding heat flow measured in comparison with a blank experiment where the M solution is injected into a M PIB solution at the same concentration, i.e.

where no structural transition is expected because the bis-urea sticker is the same and at the same 16" 0" 0,2" 0,4" 0,6" 0,8" 1"

η/η 0"

[M]"/"([M]"+"[X PIB ])" *"""""""""""""""""*"""""""""*"""""""""""""""""*"""""""*" 20°C" 40°C" 60°C" 80°C"

Heat"flow"(µcal/s)" Time"(h)" Bottom: corresponding enthalpograms.

!55# !45# !35# !25# !15# 0# 0,1# 0,2# kJ/mol#of#M#injected# [M]#/#[M PIB #or#X PIB ]#
Finally, the structure of the assemblies was probed directly by Small Angle Neutron Scattering (SANS) to confirm our interpretation of the transition between two supramolecular structures.

Figure 9 shows the scattered intensity for various M/X PIB compositions at 20 and at 80 °C. At 20 °C, all curves are characterized at low angles by a q -1 dependence of the scattered intensity over one decade. This feature is typical of long and rigid objects. In addition, at larger angles, the systems that contain at least 10% of X PIB show a q -2 dependence of the scattered intensity that is typical of solvated polymer chains. The combination of the q -1 dependence at low q and the q -2 dependence at large q is a signature for the presence of rod-like objects that are decorated with solvated chains, i.e., hairy rods. [START_REF] Pedersen | Form Factors of Block Copolymer Micelles with Spherical, Ellipsoidal and Cylindrical Cores[END_REF] Importantly, the fact that the scattered intensity for the mixtures containing 25% or 50% M is higher than the intensities for both pure solutions proves that M and X PIB co-assemble into a different and thicker structure than the pure components.

Indeed, the scattered intensity in the q -1 region depends only on the contrast and on the mass per unit length of the rods. Since the contrast of the mixture is the average of the contrast for the pure components, only a change in mass per unit length of the objects can explain the experimental data. At 80 °C, all curves (except for pure M) reach a constant intensity at low angles showing that the objects are much shorter than at 20°C.

At the concentration of the experiment (6 g/L), the interactions between objects can be neglected. [START_REF] Lortie | Structural and Rheological Study of a Bis-Urea Based Reversible Polymer in an Apolar Solvent[END_REF] Therefore, the fit of the data was attempted with various form factors of rigid cylinders. The simplest one that can describe the whole data set is the form factor for a rigid cylinder with an elliptical cross-section and a uniform contrast. Only three parameters were adjusted for each curve: the length of the cylinder and the minor and major radii of the elliptical cross-section. The other parameters (the scattering length densities and the volume fraction) were a priori evaluated from the known composition, as explained in the Experimental Section. The fits are displayed on Figures S1 to S6 while the evolution of the geometrical parameters deduced from these fits is shown on Figure 10. At 80°C, both the length of the objects and the dimensions of the cross-section evolve gradually with the composition between those of the pure systems.

Therefore, there is no indication of any change in structure with the composition, which is in agreement with the DSC data that showed that 80°C lies above the transition for all compositions. Nevertheless, the addition of the low molar mass bis-urea M to the polymeric bisurea X PIB induces an increase of the length of the objects, but this is probably simply due to a reduced steric hindrance. In contrast, at 20°C, the cross-section of the objects is larger for some intermediate compositions (25 and 50 mol% of M) than for the pure systems. This feature can only be explained if the structure of the assembly at the molecular scale changes with the composition. Following the spectroscopic data discussed above, we propose that all the systems at 80°C and the pure systems at 20°C are assembled into the filament structure that contains a single molecule in the cross-section (Figure 2). In contrast, the 50/50 mixed system at 20°C is assembled into the tube structure that contains three molecules in the cross-section. The coherence with the SANS data is obvious for the filament structure where the elliptical shape of the cross-section comes from the shape of the monomers. Moreover, the major radius measured for M and X PIB correspond to the respective size of the molecules. At first glance, in the case of the tube structure, one could expect a circular cross-section rather than an elliptical one.

However, the presence of two X PIB molecules at the same level along the tube axis is probably disfavored because of steric repulsion between the PIB chains. Therefore, locally, the crosssection is probably not circular. Of course, at larger length scale the orientations of X PIB molecules along the tube axis are not expected to be correlated, so that the projection of the assembly along its axis is probably isotropic. For a given composition, the increase in both major and minor radii on lowering the temperature from 80 to 20°C is then probably due to the increased thickness of the structure at the scale of the urea functions (one versus three bis-ureas in the cross-section) and also to a more extended conformation of the PIB chains because of the shorter distances between X PIB molecules in the tube structure. Top: major radius (full symbol) and minor radius (hollow symbol) of the cross-section. Bottom: length of the cylinder. Blue: data at 20°C; red: data at 80°C. Note that the transition temperature between the tube and filament structure lies between 20 and 80°C only for the mixtures (see Figure 6).

Conclusion

Mixing bis-ureas with slightly different structures can yield surprisingly large synergistic effects, as shown here in the case of bis-ureas with an aromatic spacer. A strong increase in viscosity is observed when a bis-urea with the sterically demanding 2,4,6-trimethylbenzene spacer is combined with a bis-urea bearing no methyl group in position 2 of the aromatic spacer (i.e. 4-methylbenzene or 4,6-dimethylbenzene). This effect is the consequence of a change in supramolecular assembly triggered by the composition of the mixture. The mixture of complementary bis-ureas forms rod-like objects that are more stable by about 1 kJ/mol, and that are thicker than the rod-like objects formed by both parent systems.

  Transform Infrared Spectroscopy. Solution spectra were measured on a Nicolet iS10 spectrometer in KBr cells of 1.0 mm pathlength and are corrected for air, solvent and cell absorption. Isothermal titration calorimetry. Data were recorded on a Microcal VP-ITC calorimeter, injecting a 0.05 mM toluene solution of M into 0.05 mM toluene solutions of M PIB or X PIB . Injections of 10 µL over 20 seconds were performed every 600 seconds at a stirring rate of 260 rpm.

Figure 1 .

 1 Figure 1. Bis-urea structures (top). Appearance of bis-urea solutions and their 1/1 mixtures, at

N

Figure 2 .

 2 Figure 2. a) Schematic supramolecular structures for low molar mass bis-ureas, where the red dotted lines represent the hydrogen bonds between urea groups. 32,33,34 b) Proposed schematic

Figure 3 .

 3 Figure 3. Top: FTIR spectra for solutions of M, X PIB and their 1/1 mixture in toluene (5 mM,

Figure 4 .

 4 Figure 4. Relative viscosity of 40 g/L (16 mM) solutions of T PIB , M PIB or X PIB with low amounts

Figure 5 .Figure 6 .

 56 Figure 5. Top: FTIR spectra of a M/X PIB solution at 4g/L in toluene (1/1 molar fraction).

  "(kJ/mol)" T**"(°C)" [M]"/"([M]"+"[X PIB ])" Page 13 of 25 Langmuir

  concentration in the cell and in the syringe. The blank experiment yields comparably small and constant exothermic heat effects that are probably mainly due to the friction during the injections. In contrast the injection of the M solution into the X PIB solution triggers a significant exothermic signal when the composition in the calorimetric cell is rich enough in M ([M]/([M]+[X PIB ]) > 0.14). The magnitude of this exothermic signal is unfortunately difficult to analyze quantitatively because of the continuously increasing viscosity during the experiment, which entails an increasing contribution of the friction. Nevertheless, the result is in good agreement with the data of Figure6: the addition of M to the X PIB solution produces a significant heat effect only above a critical composition, which probably corresponds to the composition of the filament to tube transition (at 20 °C and 0.05 mM).

Figure 7 .

 7 Figure 7. Relative viscosity of 2 g/L solutions of M/X PIB in toluene, at various temperatures. The

Figure 8 .

 8 Figure 8. ITC data. Top: heat effect produced by injecting 10 μL aliquots of a 0.05 mM solution

Figure 9 .

 9 Figure 9. SANS intensity (I) versus scattering vector (q) for solutions of M or X PIB at a

Figure 10 .

 10 Figure 10. Plot of the parameters obtained by fitting the SANS data of Figure 9 with the form

  "/"([M]"+"[X PIB ])"

& X"

ASSOCIATED CONTENT Supporting Information. Synthetic procedure and additional FTIR and SANS data.

AUTHOR INFORMATION

Corresponding Author

* E-mail: laurent.bouteiller@upmc.fr ACKNOWLEDGMENT LLB (Saclay) is acknowledged for providing beamtime for SANS experiments. We thank BASF company for the gift of a Kerocom PIBA sample. This work benefited from DANSE software developed under NSF award DMR-0520547.


Fatemeh Entezari Heravi

Mohammad Ali Nazari

Franz Chouly

Pascal Perrier

Yohan Payan

Computational fluid dynamics in the upper airway: comparison between different models and experimental data for a simplified geometry with major obstruction

Keywords: Fully coupled fluid-structure interaction, pharyngeal airway, finite element method, laminar fluid flow, turbulent fluid flow

The present study aims at comparing different computational models used for simulating the fluidstructure interaction within an in-vitro setup resembling simplified major obstruction of pharyngeal airway. Due to the nature of the problem, i.e. air flow passing over a deformable latex surface, a fully coupled fluid-structure interaction algorithm is used. A comparison is made between two finite element models for the solid domain, one using shell and the other using volume elements. The material properties of these models follow a hyperelastic behavior. For the fluid part, laminar and various turbulence models such as standard 𝑘 -𝜀, Shear Stress Transport, SSG Reynolds Stress and BSL Reynolds Stress are compared. We evaluate the efficiency of the models and how close to the experimental data are their results. The predictions of the structural model containing volume elements showed better consistency with the experimental data. In addition, the results obtained with the standard 𝑘 -𝜀 turbulence model were the least deviated among all turbulence models.

Introduction

Numerical modelling of flows in human body is of great importance. Various studies have been focusing on the airflow along the respiratory system, particularly in the upper airway. Some of them have modeled only fluid flow assuming rigid boundaries. For instance, [START_REF] Wang | Numerical analysis of respiratory flow patterns within human upper airway[END_REF] [START_REF] Wang | Numerical analysis of respiratory flow patterns within human upper airway[END_REF] have predicted the airflow along the upper airway during a whole respiratory cycle using 𝑘 -𝜀 turbulence model for the transient analysis. In [START_REF] Mylavarapu | Validation of computational fluid dynamics methodology used for human upper airway flow simulations[END_REF] [START_REF] Mylavarapu | Validation of computational fluid dynamics methodology used for human upper airway flow simulations[END_REF], various turbulence models such as unsteady Large Eddy Simulation (LES) model and some of the steady Reynolds-Averaged Navier-Stokes (RANS) models, i.e. 𝑘 -𝜀, 𝑘 -𝜔, 𝑘 -𝜔-based Shear Stress Transport, and Spallart-Alarmas models, are compared on the basis of flow computation in the upper airway. In the aforementioned study, the peak expiratory flow rate goes through the upper airway and the regions of minimum negative pressure are determined. These regions are those where the obstruction is most likely to occur. The characteristics of the flow pattern in the inspiratory phase are discussed in [START_REF] Sung | Customized three-dimensional Computational fluid dynamics simulation of the upper airway of obstructive sleep apnea[END_REF] [START_REF] Sung | Customized three-dimensional Computational fluid dynamics simulation of the upper airway of obstructive sleep apnea[END_REF] and the regions prone for obstruction in nasal cavity and in the pharynx are determined using RNG 𝑘 -𝜀 model, which is a modified version of the standard 𝑘 -𝜀 model that takes account of the contribution of different turbulence scales rather than considering only the length scale. The Spalart-Allmaras model is preferred by [START_REF] Nithiarasu | Steady flow through a realistic human upper airway geometry[END_REF] [START_REF] Nithiarasu | Steady flow through a realistic human upper airway geometry[END_REF] and they claim that this model is commonplace for incompressible flow calculations. Their study resulted in specifying the regions of minimum negative pressure by studying the steady inspiratory flow. In [START_REF] Lucey | Measurement, reconstruction, and flow-field computation of the human pharynx with application to sleep apnea[END_REF] [START_REF] Lucey | Measurement, reconstruction, and flow-field computation of the human pharynx with application to sleep apnea[END_REF], the Shear Stress Transport Turbulence (SST) model is employed to predict the airflow pattern during the inhalation phase of the respiratory cycle within a 3D pharyngeal geometry reconstructed from the data of optical coherence tomography (OCT) which could capture multi-sectional images of a tissue structure [START_REF] Fujimoto | Optical Coherence Tomography: An emerging technology for biomedical imaging and optical biopsy[END_REF] [START_REF] Fujimoto | Optical Coherence Tomography: An emerging technology for biomedical imaging and optical biopsy[END_REF]). Moreover, in this study, the effect of different wall positions is taken into account through a change of the initial rigid wall geometry. Based on a brief comparison of the results of laminar flow model, 𝑘 -𝜀 model, and 𝑘 -𝜔 SST turbulence model they concluded that the SST model would be more appropriate than the other three. The response of the upper airway to the mandibular advancement splint in treating major obstructions of the airway is investigated in Zhao et al. (2013) [START_REF] Zhao | Computational fluid dynamics for the assessment of upper airway response to oral apliance treatment in obstructive sleep apnea[END_REF]. In this study, the 𝑘 -𝜔 based Shear Stress Transport model is considered appropriate for the upper airway flow. The effect of nasal obstruction in positive airway pressure treatment therapy is studied using Spalart-Allmaras model in [START_REF] Wakayama | Effect of nasal obstruction on continuous positive airway pressure treatment: Computational Fluid Dynamics Analyses[END_REF] [START_REF] Wakayama | Effect of nasal obstruction on continuous positive airway pressure treatment: Computational Fluid Dynamics Analyses[END_REF]. In [START_REF] Cisonni | Effect of the velopharynx on intraluminal pressures in reconstructed pharynges derived from individuals with and without sleep apnea[END_REF] [START_REF] Cisonni | Effect of the velopharynx on intraluminal pressures in reconstructed pharynges derived from individuals with and without sleep apnea[END_REF] the pressure drop in the velopharynx of subjects with and without obstructions in the upper airway was computed using fluid flow simulations with 𝑘 -𝜔 SST modeling on the anatomies constructed from optical coherence tomography (OCT).

The effect of fluid pressure on deformable parts of the upper airway has been studied in [START_REF] Pelteret | Development of a computational biomechanical model of the human upper airway soft tissues toward simulating obstructive sleep apnea[END_REF] [START_REF] Pelteret | Development of a computational biomechanical model of the human upper airway soft tissues toward simulating obstructive sleep apnea[END_REF] and Carrigy et al. 2015 [START_REF] Carrigy | Simulation of muscle and adipose tissue deformation in the passive human pharynx[END_REF]. In Pelteret and Reddy (2014) [START_REF] Pelteret | Development of a computational biomechanical model of the human upper airway soft tissues toward simulating obstructive sleep apnea[END_REF] the fluid pressure changes sinusoidally to represent the effect of inhaling and exhaling but its magnitude is considered uniform except in the mouth region where it decreases linearly from the tip of tongue to the pharyngeal part of the tongue to model the inhaling through the mouth. The simulation is then completed using FEM analysis. In [START_REF] Carrigy | Simulation of muscle and adipose tissue deformation in the passive human pharynx[END_REF] [START_REF] Carrigy | Simulation of muscle and adipose tissue deformation in the passive human pharynx[END_REF] a uniform negative pressure is applied on the whole model to study the effect of area changes in velar and oral parts. Then the aforementioned authors determined the most appropriate Young's Moduli of the muscles and the adipose tissues assuming a linear elastic model.

In other studies the assumption of rigid boundaries has been replaced by a more realistic one. There are several works where both the airflow and its deformable pathway are modeled using fluidstructure interaction (FSI) method. In Huang et al. (2008) [START_REF] Huang | Computational simulation of human upper airway collapse using a pressure-/state-dependent model of genioglossal muscle contraction under laminar flow conditions[END_REF] a two-dimensional finite-element model of the upper airway is examined under laminar flow conditions taking into account the muscle's active properties. A 2D FSI simulation is carried out upon a simplified geometry for both normal and abnormal pharyngeal airways in [START_REF] Wang | Numerical analysis of respiratory flow patterns within human upper airway[END_REF] [START_REF] Wang | Numerical analysis of respiratory flow patterns within human upper airway[END_REF]. In Chouly et al. (2008) [START_REF] Chouly | Numerical and experimental study of expiratory flow in the case of major obstructions with fluid-structure interaction[END_REF] the FSI method is used for coupling a 3D deformable solid model and a 2D fluid model of the expiratory flow in the upper airway during major obstructions. In this study, flow features are analyzed with Reduced Navier-Stokes/Prandtl (RNS/P) equations which are described in [START_REF] Lagrée | The RNS/Prandtl equations and their link with other asymptotic descriptions: application to the wall shear stress scaling in a constricted pipe[END_REF] [START_REF] Lagrée | The RNS/Prandtl equations and their link with other asymptotic descriptions: application to the wall shear stress scaling in a constricted pipe[END_REF]. These equations are also used for the analysis of the flow in a 2D model of the oropharynx in [START_REF] Chouly | Simulation of the retroglossal fluid-structure interaction during obstructive sleep apnea[END_REF] [START_REF] Chouly | Simulation of the retroglossal fluid-structure interaction during obstructive sleep apnea[END_REF]. A 3D model of the expiratory airflow is simulated with fully coupled FSI method in [START_REF] Rasani | Simulation of pharyngeal airway interaction with airflow using low-Re turbulence model[END_REF] [START_REF] Rasani | Simulation of pharyngeal airway interaction with airflow using low-Re turbulence model[END_REF]. 𝑘 -𝜔 based Shear Stress Transport model is adopted for this study. In [START_REF] Heil | Flows in deformable tubes and channels[END_REF] [START_REF] Heil | Flows in deformable tubes and channels[END_REF], a review is published that presents numerous models used for investigating general characteristics of airflow in interaction with soft tissues throughout the human body. These models range from zero-equation models to threedimensional Navier-Stokes simulations. This publication emphasizes the necessity and importance of fully-coupled approaches for dealing with such problem. Such a fully coupled FSI simulation is used by [START_REF] Wang | Fluid-structure interaction modeling of upper airways before and after nasal surgery for obstructive sleep apnea[END_REF] [START_REF] Wang | Fluid-structure interaction modeling of upper airways before and after nasal surgery for obstructive sleep apnea[END_REF] for determining the efficiency of nasal surgery, where the standard 𝑘 -𝜔 turbulence model is employed for simulating the flow domain. The same is done by Zhao et al. (2013) [START_REF] Zhao | Simulation of upper airway occlusion without and with mandibular advancement in obstructive sleep apnea[END_REF] for the study on mandibular advancement method, using the 𝑘 -𝜔 based Shear Stress Transport. The deformation along the entire upper airway is studied by transient flow analysis coupled with linear structural analysis in [START_REF] Kim | Numerical investigation of flow-induced deformation along the human respiratory upper airway[END_REF] [START_REF] Kim | Numerical investigation of flow-induced deformation along the human respiratory upper airway[END_REF] using Large Eddy Simulation model for the transient analysis of the fluid domain. In [START_REF] Pirnar | Computational fluid-structure interaction simulation of flow in the human upper airway[END_REF] [START_REF] Pirnar | Computational fluid-structure interaction simulation of flow in the human upper airway[END_REF] both soft palate flutter and the airway narrowing are simulated through fully coupled FSI method. In all of these fluid-structure interaction studies, the solid domain is considered linearly elastic.

Given the diversity of the models available for simulating fluid-structure interaction in the upper airways, and the variety of their predictions in the context of the study of the sleep apnea syndrome, a precise and quantitative assessment of these different approaches is needed. In this aim, we propose to take the measures carried out by Chouly et al. (2008) with a replica of the pharyngeal airway as reliable reference data to which the predictions of the models will be precisely and quantitatively compared. The advantage of using these experimental results is that they were obtained in controlled conditions through precise measurement techniques which could not be achieved while working with patientspecific data. In order for increasing the accuracy of the results, a fully coupled method is chosen as the fluid-structure coupling algorithm, and hyperelastic material properties are assigned to the solid domain. In previous works Chouly et al. (2008) [START_REF] Chouly | Numerical and experimental study of expiratory flow in the case of major obstructions with fluid-structure interaction[END_REF] using RNS/P equations provided predictions of the pressure drop within the constriction of the pharyngeal replica, and of the flow-induced deformation, that fitted well the experimental data. Nevertheless, in order to go further and carry out later on simulations from real and complex patient-specific geometries, other fluid flow models should be considered. Indeed RNS/P equations are mostly limited to simple 2D geometries and induce numerical difficulties in case of heavy recirculation after the constriction. Therefore, complete 3D Navier-Stokes equations with laminar or turbulent closures should be considered, as they are now available in standard finite element packages such as ANSYS/FLUENT(TM). In addition, since various turbulent models are provided in ANSYS/CFX 15.0 and are used in biomedical simulations, their predictive performance and accuracy for studying such a complex phenomenon need to be compared using relevant experimental data. The efficiencies of two different finite element models for the solid domain are examined; one containing shell elements and, the other, volume elements. Some of the most commonplace flow models, i.e. laminar, standard 𝑘 -𝜀, Shear Stress Transport, and BSL and SSG Reynolds Stress models, are compared to each other in this study ("BSL" and "SSG" stand for "Baseline" and "Speziale, Sarkar, Gatski" respectively). The pressure value at the point where the constriction is the narrowest is used as a key variable for comparing the results of the numerical models and the experimental results reported in Chouly et al. (2008) [START_REF] Chouly | Modelling human pharyngeal airway: validation of numerical simulations using in vitro experiments[END_REF].

Theoretical assumptions and materials

Governing Equations

For appropriate simulations of fluid-structure interaction problems such as highly deformable tissues interacting with turbulent flows, an iterative fully coupled transient algorithm is suggested where fluid and solid equations are solved by separate solvers. Especially having large deformations in solid domain would affect the boundary conditions of the fluid flow. In this algorithm, two constraints are set at the interface of the fluid and solid domains for coupling. The first constraint is the relationship between the displacement in the solid domain and the fluid velocity at each interface node as:

t    d u (1)
where 𝐮 is the fluid velocity vector and 𝐝 is the displacement vector of the structure.

The second coupling constraint is the equivalence of the Cauchy stress tensors of the two domains as,

0   s s f f n σ n σ (2)
where 𝐧 𝑠 and 𝐧 𝑓 are the outward unit normal vectors on the solid and fluid surfaces at their interface, f σ and s σ are the Cauchy stress tensors of the fluid and solid domains respectively. Cauchy stress tensor in solid follows a neo-Hookean hyperelastic constitutive law:

I I B B σ ) 1 ( ) ) ( 3 1 ( 3 / 5 1     J K tr J C s ( 3
)
where I is the unit second order stress tensor, B is the left Cauchy-Green strain tensor, and 𝐶 1 , and 𝐾 are shear modulus, and bulk modulus respectively. 𝐽 is the Jacobian of deformation and is equal to √ det (𝐁). This formulation of the stress is derived from the strain energy function for the nearly incompressible isotropic neo-Hookean hyperelastic material defined in [START_REF] Nazari | Simulation of dynamic orofacial movements using a constitutive law varying with muscle activation[END_REF] [START_REF] Nazari | Simulation of dynamic orofacial movements using a constitutive law varying with muscle activation[END_REF] as:

𝑊 = 𝐶 1 2 (𝐽 -2/3 𝐼 1 -3) + 𝐾 2 ( 𝐽 -1 ) 2 (4)
where 𝐼 1 = 𝑡𝑟(𝐁) is the first invariant of the Cauchy-Green strain tensor.

The Cauchy stress tensor in the Stokes fluid flow is:

) ) ( 3
1 ( 2 I V V I σ tr p f      (5)
where𝐕 is the strain rate tensor and 𝑝 is the fluid pressure. 𝐕 is defined as:

𝐕 = 1 2 (∇𝐮 + ∇𝐮 𝐓 ) (6)
where 𝐮 stands for the fluid velocity, and ∇ is gradient.

The general characteristics of the flow are considered to be Newtonian, inviscid, and incompressible. Therefore the incompressible Navier-Stokes equations are,

𝛁 • 𝐮 = 0 (7) 𝜌(𝐮 • 𝛁)𝐮 = -𝛁𝑃 + 𝛁(2𝜇𝐕) ( 8
)
where 𝐮 is the fluid velocity, 𝜌 is the fluid density, 𝜇 is the dynamic viscosity, and 𝐕 is the strain rate tensor.

For turbulent flows, the velocity vector and pressure are considered as the sum of a mean (𝐔, 𝑃), and a fluctuating part (𝐮 ′ , 𝑝'):

𝐮 = 𝐔 + 𝐮′ (9 a) p=P+p' ( 9 b)
With Reynolds time averaging defined as,

𝐟 ̅ = lim T→∞ 1 𝑇 ∫ 𝐟(𝑠) 𝑡+𝑇 𝑡 𝑑𝑠 (10)
and by applying this type of averaging on Eq. ( 7) and Eq. ( 8) using Eqs. ( 9), the Reynolds averaged equation of mass conservation, and the Reynolds Averaged Navier-Stokes (RANS), are derived as:

𝛁 • 𝐔 = 0 (11) 𝜌(𝐔 • 𝛁)𝐔 = -𝛁𝑃 + 𝛁 [ 2𝜇𝐕 -𝜌𝝉 ] ( 12
)
where 𝜌𝝉 is called Reynolds Stress tensor. The components of this tensor are equal to 𝜌𝐮 ′ ⨂𝐮′(⨂ stands for the dyadic product of two vectors) and the way one defines it determines the type of the flow model with which the simulation is to be done (Wilcox(1993) [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF]). For laminar flow, this term would be zero, because the velocity field is assumed to have no fluctuation in any direction. In turbulence modeling, there are two major ways for defining this term. In the first way, the Eddy Viscosity model assumes that Reynolds stresses are proportional to mean velocity gradients and are defined as:

𝜌𝐮 ′ ⨂𝐮′ = 𝜇 𝑡 𝐕 - 2 3 𝐈 ( 𝜌𝑘 ) (13)
where 𝑘 is the turbulence kinetic energy defined as 𝑘 =

𝐮′•𝐮′ ̅̅̅̅̅̅̅ 𝟐

, and 𝜇 𝑡 is the eddy viscosity or turbulent viscosity; in the 𝑘 -𝜀 model and the Shear Stress Transport Model 𝜇 𝑡 is respectively defined as:

𝜇 𝑡 = 𝐶 𝜇 𝜌 𝑘 2 𝜀 (14)
𝜇 𝑡 = 𝜌 𝛼𝑘 max ( 𝛼𝜔, 𝑆𝐹 ) [START_REF] Chouly | Simulation of the retroglossal fluid-structure interaction during obstructive sleep apnea[END_REF] where 𝐶 𝜇 and 𝛼 are constant parameters, 𝜀 is the turbulence eddy dissipation rate, 𝜔 is the specific dissipation rate, 𝑆 is an invariant of the strain rate matrix, and 𝐹 is a blending function. These variables are defined trough the following equations:

𝜀 = 𝜈𝛁𝐮 ′ : 𝛁𝐮′ ̅̅̅̅̅̅̅̅̅̅̅ ( 16
) 𝜔 = 𝜀 𝑘 (17) 𝐹 = tanh ( 𝑓 2 ) ( 18
)
where 𝑓 is defined as,

𝑓 = max ( 2√𝑘 𝛽𝜔𝑦 , 500𝜈 𝑦 2 𝜔 ) (19)
where ":" shows double contraction of two second-order tensors, 𝛽 is a constant, 𝑦 is the distance to the nearest wall, and 𝜈 is the kinematic viscosity.

The second common way for modeling Reynlods Stress Tensor term is called Reynolds Stress model. It is based on solving transport equations for all Reynolds stress tensor components and the dissipation rate. Two of the most commonplace methods of this type are SSG Reynolds Stress model and BSL Reynolds Stress model whose formulations could be found in [START_REF] Eisfeld | Reynolds stress modelling for complex aerodynamic flows[END_REF] [START_REF] Eisfeld | Reynolds stress modelling for complex aerodynamic flows[END_REF].

The simulation is conducted through three levels of iteration, (a) the field loop, used to converge toward a solution for each field with the associated solver, (b) the coupling loop, used for exchanging loads or displacements between the two fields, and, (c) the time loop which intends to make each step move forward in time. Within each field loop, the governing equations of the flow, Eq. ( 7) and Eq. ( 8), are solved for finding the value of pressure and velocity vector at nodes. During the coupling (stagger) loop, the results of the fluid field solver are applied to the relevant nodes of the solid domain using Eq. ( 2). The solid field solver computes the deformation of the structure and the results are applied on the fluid interface with solid domain as a new set of boundary conditions during the coupling loop using Eq. ( 1). This iterative coupling procedure continues until the convergence tolerance or the maximum number of iterations is reached. For a steady-state FSI simulation the solution for fluid domain is quasistatic and for solid domain is transient within each field loop (Gianopapa, 2004 [26]).

Material properties

We have seen above that the solid domain is considered to follow the neo-Hookean material model. Solid domain is of latex material and its elastic parameters are 𝐸 = 1.68MPa and ν = 0.499 (Chouly et al., 2008). Linearization of Eq. ( 4) at initial point gives,

𝐶 1 = 𝐸/3 (20)
And with the assumption of nearly incompressible material the material constants are given in Table 1. The fluid material is set as air at 25℃ with 𝜌 = 1.185 𝑘𝑔 𝑚 3 and 𝜇 = 1.831 × 10 -5 𝑘𝑔 𝑚𝑠 .

Finite Element Model

The geometry of the pharyngeal replica setup was created using the dimensions reported in Chouly et al. (2008) (Fig. 1). The model is build up in ANSYS release 15. Fig. 2 represents the finite element model of the latex surface. The mesh in Fig. 2(a) consists of a layer of 8-node quadrilateral SHELL181 elements, and the one in Fig. 2(b) consists of 20-node hexahedral SOLID186 elements along with 10-node SOLID 187 elements. The SOLID187 elements are distributed through the thickness and they count to 337 out of 1760 total elements. The number of elements and nodes of these two models are given in Table 2. The thickness of latex layer is 0.3𝑚𝑚. As shown, the latex surface is divided into two regions. The core region is in interaction with the airflow and the surrounding region is attached to the rigid pipe. Water flows in this pipe to apply a constant pressure against the core region. In order to simulate the conditions of the experimental setup an external pressure of 300Pa is applied on the upper surface of the core region and its inferior surface is supposed to act as the interface of the fluid and structure domains. The nodes of the surrounding region are fixed. The total time for simulation is considered to be 1s. The solution time step size is set as a function of inlet pressure to overcome convergence issues. It varies between 1e-5 to 1e-3 second. The finite element model of the setup airway is depicted in Fig. 3. The reason for using this relatively coarse mesh is to avoid using a mesh that captures the fine scales in the fluid in the turbulent models such as 𝑘 -𝜀 method. Indeed such fine scales generate turbulent vortices that unrealistically dissipate in the solution yielding an inaccurate simulation. The inlet total pressure, 𝑃 𝑖 , ranges from 100 Pa to 600 Pa and the outlet static pressure is considered to be 0 Pa. All pressure values are taken with respect to atmosphere. The walls of the airway are rigid except the section interacting with the solid domain.

Results

Influence of the mesh used for the solid domain

The pressure at the constriction, 𝑃 𝑐 , is used as a basis for comparison between the experimental setup and the model. Fig. 4 illustrates a comparison of 𝑃 𝑐 resulting from using either shell or volume elements to model the solid domain. These results are also presented in Table 3. The flow model was considered to be 𝑘 -𝜀 in this case. The coefficient of determination, 𝑅 2 , was computed as 0.92 for the model with volume elements, and 0.78 for that with shell elements. Comparing these values, we observed a trend for the volume element to work better for simulating such thin hyper-elastic surfaces interacting with the passage of fluid flow.

Influence of the flow model

In Fig. 5, the value of 𝑃 𝑐 is indicated vs. different values of inlet pressure using each one of the flow models described in section 2.1. The three-dimensional structure is chosen for the solid domain in this comparison. The results are also presented in Table 4. The 𝑅 2 parameter is calculated for these models and is given in Table 3. As can be seen in Fig. 5, and according to the value of 𝑅 2 for each case, the standard 𝑘 -𝜀 model seems to be slightly better than other approaches for modeling the fluid. However, it should be noticed that by using the laminar model, the simulations would be less complicated due to the elimination of time-dependent fluctuating turbulent parameters. In addition, the coefficient of determination in Chouly et al. (2008) [START_REF] Chouly | Modelling human pharyngeal airway: validation of numerical simulations using in vitro experiments[END_REF] is computed as 0.58 for an applied external pressure of 300Pa. Comparing the results obtained from solving a 2D flow model in Chouly et al. (2008) [START_REF] Chouly | Modelling human pharyngeal airway: validation of numerical simulations using in vitro experiments[END_REF] with those reported in Table 5, the necessity of 3D simulation for modeling the fluid-structure interaction problems is justified. The contour of deformation for the interface is shown in Fig. 6. This contour is the result of a 400 Pa inlet pressure set for the simulation with volume elements and 𝑘 -𝜀 flow model. The flow streamlines are also illustrated in Fig. 7. The streamline pattern before and after the constriction region justifies the asymmetric deformation of the pharyngeal part of the tongue. Maximum deformation has occurred after the constriction zone where the negative pressure over the surface is maximum. Therefore, the probable occlusion has the highest chance to occur after the initial stenosis region.

Conclusion

The goal of this study was to assess the accuracy of various 3D flow models in predicting the airflow properties of the pharyngeal airway in the case of major obstruction using a fully coupled fluid-structure interaction method. To this purpose, we relied on precise and quantitative experimental data coming from an in-vitro replica of the pharyngeal airway [START_REF] Chouly | Numerical and experimental study of expiratory flow in the case of major obstructions with fluid-structure interaction[END_REF].The solid material was assumed to be an isotropic, incompressible neo-Hookean hyperelastic one. The performances of shell and 3D finite element were compared with each other in the solid domain. Also, for the fluid domain, a comparison was made between laminar flow model and four Reynolds-Averaged Navier-Stokes (RANS) models, namely the standard 𝑘 -𝜀, the 𝑘 -𝜔 based Shear Stress Transport, the BSL Reynolds Stress, and the SSG Reynolds Stress model. The pressure at the constriction point was obtained from all numerical models and was compared with the experimental data. Nearly all of the models provided correct results. Among these, the results of standard 𝑘 -𝜀 flow model fitted best; however, the accuracy of laminar model results was reasonable and this model could be used in order for avoiding difficulties associated with tuning parameters in turbulence modeling. Both volume and shell elements worked well; though it could be seen that volume elements give slightly better results for modeling the solid. This finding is consistent with the results published in [START_REF] Landajuela | Coupling schemes for FSI forward prediction challenge: comparative study and validation[END_REF] [START_REF] Landajuela | Coupling schemes for FSI forward prediction challenge: comparative study and validation[END_REF]).

The result of this work could be further applied on the realistic human oropharyngeal model for the study and analysis of the Syndrome of Obstructive Sleep Apnea (OSA). A transient study could be done for simulating the respiratory cycles for both normal and patient cases. Also, possible treatments could be simulated before they are examined on a patient in order to reduce possible side effects and risk factors.

Figure 1 Figure 2

 12 Figure 1 Geometry of the pharyngeal replica setup

Figure 3

 3 Figure 3 Finite element model of the simplified pharyngeal airway Fluid mesh consists of 2082 nodes in 2216 finite volume elements from which 1222 are hexahedral.The reason for using this relatively coarse mesh is to avoid using a mesh that captures the fine scales in the fluid in the turbulent models such as 𝑘 -𝜀 method. Indeed such fine scales generate turbulent vortices that unrealistically dissipate in the solution yielding an inaccurate simulation. The inlet total

Figure 4 .

 4 Figure 4. Comparison of the results of shell and volume elements

Figure 5

 5 Figure 5 Comparison of the results of different fluid models. The solid domain was modeled with volume elements.

Figure 6 Figure 7

 67 Figure 6 Deformation of latex surface (in meter) for 𝑃 𝑖 = 400Pa (displacements scaled for better visualization)

Table 1

 1 Constants of the neo-Hookean 𝑾 function

	𝐶 1 (𝑃𝑎)	𝐾/2 (1/𝑃𝑎)
	560000	7.1429𝑒 -09

Table 2

 2 Number of nodes of the solid models

	2D Shell Model	3D Model
	3239	7819

Table 3

 3 Comparison of the results of shell and volume elements

		350
		300
		250
		200
	(Pa)	150
	Pc	100
		50
		0
		-50	0	100	200	300	400		500	600	700
		-100				Inlet Pressure (Pa)
					Experimental Data	SHELL-k-eps	SOLID-k-eps
							𝑃 𝑐 (Pa)
		Inlet pressure (Pa)	Experimental Data	Shell elements, 𝑘 -𝜀	Volume elements, 𝑘 -𝜀
		100			3.2		-33.0			-31.2
		200			-9.1		-39.3			-38.9
		300			-11.2		-20.8			-17.6
		400			13.6		32.5			42.7
		500			133		102.7			143.0
		600			290		177.2			237.7

Table 4

 4 Comparison of the results of different fluid models. The solid domain was modeled with volume elements.

		350
		300
		250
		200
	(Pa)	150
	Pc	100
		50
		0
		0	100	200	300	400	500	600	700
		-50
		-100				Pi (Pa)
		Experimental Data k-eps laminar SSG Reynolds Stress Shear Stress Transport BSL Reynolds Stress
						𝑃 𝑐 (Pa)
	Inlet pressure	Experimental	𝑘 -𝜀	Laminar	SSG Reynolds	Shear Stress	BSL Reynolds
		(Pa)	Data			Stress	Transport	Stress
		100	3.2	-31.2	-42.6	-47.0	-39.7	-38.7
		200	-9.1	-38.9	-57.6	-63.3	-54.3	-53.0
		300	-11.2	-17.6	-40.0	-45.3	-35.8	-34.6
		400	13.6	42.7	21.0	16.4	25.0	26.0
		500	133	143.0	123.5	118.3	126.8	127.7

Table 5

 5 𝑹 𝟐 parameters calculated for the different flow models interacting with the 3D solid mesh

	Laminar	𝑘 -𝜀	Shear Stress Transport	SSG Reynolds Stress	BSL Reynolds Stress
	0.86	0.92	0.88	0.84	0.88

Acknowledgement

This work has been supported by Center for International Scientific Studies and Collaboration (CISSC) and French Embassy in Tehran in the framework of the project Hubert Curien Gundishapur BIOSTOC.