|
|
(48 intermediate revisions by 3 users not shown) |
Line 1: |
Line 1: |
| = A fluid-structure interaction benchmark in turbulent flow (FSI-PfS-1a) = | | =Fluid-structure interaction in turbulent flow past cylinder/plate configuration I (First swiveling mode)= |
| {{UFRHeader | | {{UFRHeader |
| |area=2 | | |area=2 |
| |number=13 | | |number=13 |
| }} | | }} |
| | |
| | |
| __TOC__ | | __TOC__ |
| = Introduction = | | |
| | == Flows around bodies == |
| | === Underlying Flow Regime 2-13 === |
| | |
| | = Description = |
| | |
| | == Introduction == |
|
| |
|
| A flexible structure exposed to a fluid flow is deformed and deflected | | A flexible structure exposed to a fluid flow is deformed and deflected |
Line 22: |
Line 30: |
| simulation tools are required. | | simulation tools are required. |
|
| |
|
| The long-term objective of the present research project is the coupled | | The long-term objective of the research reported here is the coupled |
| simulation of big lightweight structures such as thin membranes | | simulation of big lightweight structures such as thin membranes |
| exposed to turbulent flows (outdoor tents, awnings...). To study these | | exposed to turbulent flows (outdoor tents, awnings...). To study these |
| complex FSI problems, a multi-physics code framework was recently | | complex FSI problems, a multi-physics code framework was recently |
| developed (Breuer et al., 2012). In order to assure reliable numerical | | developed (Breuer et al., 2012) combining Computational Fluid Dynamics (CFD) and Computational Structural Dynamics (CSD) solvers . In order to assure reliable numerical |
| simulations of complex configurations, the whole FSI code needs to be | | simulations of complex configurations, the whole FSI code needs to be |
| validated at first on simpler test cases with trusted reference | | validated at first on simpler test cases with trusted reference |
Line 32: |
Line 40: |
| developed is detailed. The CFD and CSD solvers were at first checked | | developed is detailed. The CFD and CSD solvers were at first checked |
| separately and then, the coupling algorithm was considered in detail | | separately and then, the coupling algorithm was considered in detail |
| based on a laminar benchmark. A 3D turbulent test case was also taken | | based on a laminar benchmark. A 3D turbulent test case was also calculated to prove the applicability of the newly developed |
| into account to prove the applicability of the newly developed
| |
| coupling scheme in the context of large-eddy simulations | | coupling scheme in the context of large-eddy simulations |
| (LES). However, owing to missing reference data a full validation was | | (LES). However, owing to missing reference data a full validation was |
Line 47: |
Line 54: |
| community for the technically relevant case of turbulent flows | | community for the technically relevant case of turbulent flows |
| interacting with flexible structures. | | interacting with flexible structures. |
| | |
| | == Review of previous work == |
|
| |
|
| The present study is mainly related to two former investigations | | The present study is mainly related to two former investigations |
| of Turek and Hron (2006,2010) and Gomes et al. (2006, 2012) on vortex-induced fluid-structure interactions. | | of Turek and Hron (2006, 2010) and Gomes et al. (2006, 2012) on vortex-induced fluid-structure interactions. |
| The well-known 2D purely numerical laminar benchmarks of Turek and Hron (2006, 2010) developed in a collaborative research effort on | | The well-known 2D purely numerical laminar benchmarks of Turek and Hron (2006, 2010) developed in a collaborative research effort on |
| FSI (DFG Forschergruppe 493) consists of an elastic cantilever | | FSI (DFG Forschergruppe 493) consists of an elastic cantilever |
Line 88: |
Line 97: |
| structural simulation and the grid adaptation of the flow | | structural simulation and the grid adaptation of the flow |
| prediction. | | prediction. |
| | |
| | == Choice of test case == |
|
| |
|
| Thus, in the present study a slightly different configuration is | | Thus, in the present study a slightly different configuration is |
Line 102: |
Line 113: |
| for the coupled problem) and the processing of the respective data to | | for the coupled problem) and the processing of the respective data to |
| guarantee a reliable reproduction of the proposed test case with | | guarantee a reliable reproduction of the proposed test case with |
| various suitable methods. A detailed description of the present test case is published in De Nayer et. al. (2013). | | various suitable methods. A detailed description of the present test case is published in De Nayer et al. (2014). |
|
| |
|
| The described test case FSI-PfS-1a is a part of a series of reference | | The described test case FSI-PfS-1a is a part of a series of reference |
Line 113: |
Line 124: |
| the structure deflections are completely two-dimensional and | | the structure deflections are completely two-dimensional and |
| larger. | | larger. |
| | | |
| = Test case description =
| |
| == Description of the geometrical model and the test section ==
| |
| | |
| FSI-PfS-1a consists of a flexible thin structure with a distinct
| |
| thickness clamped behind a fixed rigid non-rotating cylinder installed
| |
| in a water channel (see Fig. 1). The
| |
| cylinder has a diameter D=0.022m. It is positioned in the
| |
| middle of the experimental test section with <math>H_c \operatorname{=} H/2 \operatorname{=} 0.120\,m</math> <math>(H_c/D \approx 5.45)</math>, whereas the test section denotes a
| |
| central part of the entire water channel (see
| |
| Fig. 2). Its center is located at <math>L_c \operatorname{=}
| |
| 0.077\,m</math> <math>(L_c/D \operatorname{=} 3.5)</math> downstream of the inflow
| |
| section. The test section has a length of <math>L \operatorname{=} 0.338\,m</math>
| |
| <math>(L/D \approx 15.36)</math>, a height of <math>H \operatorname{=} 0.240\,m</math>
| |
| <math>(H/D \approx 10.91)</math> and a width <math>W \operatorname{=} 0.180\,m</math>
| |
| <math>(W/D \approx 8.18)</math>. The blocking ratio of the channel is
| |
| about <math>9.2\%</math>. The gravitational acceleration <math>g</math> points in
| |
| x-direction (see Fig.~\ref{fig:rubber_plate_geom}), i.e. in the
| |
| experimental setup this section of the water channel is turned 90
| |
| degrees. The deformable structure used in the experiment behind the
| |
| cylinder has a length <math>l \operatorname{=} 0.060\,m</math> <math>(l/D \approx 2.72)</math>
| |
| and a width <math>w \operatorname{=} 0.177\,m</math> <math>(w/D \approx 8.05)</math>.
| |
| Therefore, in the experiment there is a small gap of about <math>1.5
| |
| \times 10^{-3}\,m</math> between the side of the deformable structure and
| |
| both lateral channel walls.
| |
| The thickness of the plate is <math>h \operatorname{=} 0.0021\,m</math> <math>(h/D
| |
| \approx 0.09)</math>. This thickness is an averaged value. The material
| |
| is natural rubber and thus it is difficult to produce a perfectly
| |
| homogeneous 2 mm plate. The measurements show that the thickness of
| |
| the plate is between 0.002 and 0.0022 m. All parameters of the
| |
| geometrical configuration of the FSI-PfS-1a benchmark are summarized as follows:
| |
| | |
| [[File:table3.png]]
| |
| | |
| [[File:FSI-PfS-1a_Benchmark_Rubberplate_geometry0001_new.png]]
| |
| | |
| Fig. 1: Geometrical configuration of the FSI-PfS-1a Benchmark.
| |
| | |
| == Description of the water channel ==
| |
| | |
| | |
| In order to validate numerical FSI simulations based on reliable
| |
| experimental data, the special research unit on FSI (Bungartz et. al. (2006, 2010))
| |
| designed and constructed a water channel (Göttingen type, see
| |
| Fig. 2) for corresponding experiments with a
| |
| special concern regarding controllable and precise boundary and
| |
| working conditions Gomes et al. (2006, 2010, 2011, 2013). The
| |
| channel (2.8 m x 1.3 m x 0.5 m, fluid volume
| |
| of 0.9 m³) has a rectangular flow path and includes several
| |
| rectifiers and straighteners to guarantee a uniform inflow into the
| |
| test section. To allow optical flow measurement systems like
| |
| Particle-Image Velocimetry, the test section is optically accessible
| |
| on three sides. It possesses the same geometry as the test section
| |
| described in Fig. 1. The structure is
| |
| fixed on the backplate of the test section and additionally fixed in
| |
| the front glass plate. With a 24 kW axial pump a water inflow of up to
| |
| <math>u_{\text{max}}=6 m/s</math> is possible. To prevent asymmetries the
| |
| gravity force is aligned with the x-axis in main flow direction.
| |
| | |
| [[File:waterchannel_new.png]]
| |
| | |
| Fig. 2: Sketch of the flow channel (dimensions given in mm).
| |
| | |
| == Flow parameters ==
| |
| | |
| Several preliminary tests were performed to find the best working
| |
| conditions in terms of maximum structure displacement, good
| |
| reproducibility and measurable structure frequencies within the
| |
| turbulent flow regime.
| |
| | |
| [[File:inflow1.png]]
| |
| | |
| Fig. 3: Experimental displacements of the structure extremity versus the inlet velocity.
| |
| | |
| Fig. 3 depicts the measured extrema of the structure displacement versus the
| |
| inlet velocity and Fig. 4 gives the
| |
| frequency and Strouhal number as a function of the inlet
| |
| velocity. These data were achieved by measurements with the laser
| |
| distance sensor explained in Section [[UFR_2-13_Description#Laser distance sensor|Laser distance sensor]]. The
| |
| entire diagrams are the result of a measurement campaign in which the
| |
| inflow velocity was consecutively increased from 0 to <math>2.2 m/s</math>. At an inflow velocity of <math>u_{\text{inflow}}=1.385 m/s </math>
| |
| the displacement are symmetrical, reasonably large and well
| |
| reproducible. Based on the inflow velocity chosen and the cylinder
| |
| diameter the Reynolds number of the experiment is equal to
| |
| <math>\text{Re}=30,470</math>.
| |
| | |
| [[File:inflow2.png]]
| |
| | |
| Fig. 4: Experimental measurements of the frequency and the corresponding Strouhal number of the FSI phenomenon versus the inlet velocity.
| |
| | |
| Regarding the flow around the front
| |
| cylinder, at this inflow velocity the flow is in the sub-critical
| |
| regime. That means the boundary layers are still laminar, but
| |
| transition to turbulence takes place in the free shear layers evolving
| |
| from the separated boundary layers behind the apex of the
| |
| cylinder. Except the boundary layers at the section walls the inflow
| |
| was found to be nearly uniform (see
| |
| Fig. 5). The velocity components
| |
| <math>\overline{u} </math> and <math>\overline{v}</math> are measured with two-component
| |
| laser-Doppler velocimetry (LDV) along the y-axis in the middle of the
| |
| measuring section at <math>x/D=4.18 </math> and <math>z/D=0 </math>. It can be assumed that
| |
| the velocity component <math>\overline{w} </math> shows a similar velocity profile
| |
| as <math>\overline{v} </math>. Furthermore, a low inflow turbulence level of
| |
| <math>\text{Tu}_{\text{inflow}}=\sqrt{\frac{1}{3}~\left(\overline{u'^2}+\overline{v'^2}+\overline{w'^2} \right)}/u_{\text{inflow}}=0.02 </math> is measured. All
| |
| experiments were performed with water under standard conditions at
| |
| <math>T=20^{\circ}\,C</math>. The flow parameters are summarized in the following table:
| |
| | |
| {| class="wikitable"
| |
| |+ Flow parameters
| |
| |----
| |
| ! scope="row" | Inflow velocity
| |
| | <math>u_{\text{inflow}}=1.385\,m/s </math>
| |
| |----
| |
| ! scope="row" | Flow density
| |
| | <math> \rho_f=1000\,kg/m^3</math>
| |
| |----
| |
| ! scope="row" | Flow dynamic viscosity
| |
| | <math> \mu_f=1.0 \times 10^{-3}\,Pa\,s </math>
| |
| |}
| |
| | |
| [[File:flow_conditions.png]]
| |
| | |
| Fig. 5 Profiles of the mean streamwise and normal velocity as well as the turbulence level at the inflow section of the water channel.
| |
| | |
| == Material Parameters ==
| |
| | |
| Although the material shows a strong non-linear elastic behavior for
| |
| large strains, the application of a linear elastic constitutive law
| |
| would be favored, to enable the reproduction of this FSI benchmark by
| |
| a variety of different computational analysis codes without the need
| |
| of complex material laws. This assumption can be justified by the
| |
| observation that in the FSI test case, a formulation for large
| |
| deformations but small strains is applicable. Hence, the
| |
| identification of the material parameters is done on the basis of the
| |
| moderate strain expected and the St. Venant-Kirchhoff constitutive law
| |
| is chosen as the simplest hyper-elastic material.
| |
| | |
| The density of the rubber material can be determined to be
| |
| <math>\rho_{\text{rubber plate}}</math>=1360 kg/m<math>^3</math> for a thickness
| |
| of the plate h = 0.0021 m. This permits the accurate modeling
| |
| of inertia effects of the structure and thus dynamic test cases can be
| |
| used to calibrate the material constants. For the chosen material
| |
| model, there are only two parameters to be defined: the Young's
| |
| modulus E and the Poisson's ratio <math>\nu</math>. In order to avoid
| |
| complications in the needed element technology due to
| |
| incompressibility, the material was realized to have a Poisson's ratio
| |
| which reasonably differs from <math>0.5</math>. Material tests of the
| |
| manufacturer and complementary experimental/numerical structure test studies (static, dynamic and decay test scenarios) indicate that the Young's modulus is E=16 MPa and
| |
| the Poisson's ratio is <math>\nu</math>=0.48 (a detailed description of the structure tests is available in De Nayer et. al. (2013)).
| |
| | |
| {| class="wikitable"
| |
| |+ Structure parameters
| |
| |----
| |
| ! scope="row" | Density
| |
| | <math>\rho_{\text{rubber plate}}=1360\,kg/m^3</math>
| |
| |----
| |
| ! scope="row" | Young's modulus
| |
| | <math>E=16\,MPa</math>
| |
| |----
| |
| ! scope="row" | Poisson's ratio
| |
| | <math>\nu=0.48\,</math>
| |
| |}
| |
| | |
| = Measuring Techniques =
| |
| | |
| Experimental FSI investigations need to contain fluid and structure
| |
| measurements for a full description of the coupling process. Under
| |
| certain conditions, the same technique for both disciplines can be
| |
| used. The measurements performed by
| |
| Gomes et. al. (2006, 2010, 2013) used the same camera system for
| |
| the simultaneous acquisition of the velocity fields and the structural
| |
| deflections. This procedure works well for FSI cases involving laminar
| |
| flows and 2D structure deflections. In the present case the structure
| |
| deforms slightly three-dimensional with increased cycle-to-cycle
| |
| variations caused by turbulent variations in the flow. The applied
| |
| measuring techniques, especially the structural side, have to deal
| |
| with those changed conditions especially the formation of
| |
| shades. Furthermore, certain spatial and temporal resolutions as well
| |
| as low measurement errors are requested. Due to the different
| |
| deformation behavior a single camera setup for the structural
| |
| measurements like in Gomes et. al. (2006, 2010, 2013) used was not
| |
| practicable. Therefore, the velocity fields were captured by a 2D
| |
| Particle-Image Velocimetry (PIV) setup and the structural deflections
| |
| were measured with a laser triangulation technique. Both devices are
| |
| presented in the next sections.
| |
| | |
| === Particle-image velocimetry ===
| |
| | |
| A classic Particle-Image Velocimetry (Adrian, 1991) setup consists of a single camera obtaining two
| |
| components of the fluid velocity on a planar surface illuminated by a
| |
| laser light sheet. Particles introduced into the fluid are following the
| |
| flow and reflecting the light during the passage of the light sheet.
| |
| By taking two reflection fields in a short time interval <math>\Delta</math>t, the most-likely displacements of several particle groups on an
| |
| equidistant grid are estimated by a cross-correlation technique or a
| |
| particle-tracking algorithm. Based on a precise preliminary
| |
| calibration, the displacements obtained and the time interval
| |
| <math>\Delta</math>t chosen the velocity field can be calculated. To prevent
| |
| shadows behind the flexible structure a second light sheet was used to
| |
| illuminate the opposite side of the test section.
| |
| | |
| The phased-resolved PIV-measurements (PR-PIV) were carried out with a
| |
| 4 Mega-pixel camera (TSI Powerview 4MP, charge-coupled device (CCD)
| |
| chip) and a pulsed dual-head Neodym:YAG laser (Litron NanoPIV 200)
| |
| with an energy of 200 mJ per laser pulse. The high energy of
| |
| the laser allowed to use silver-coated hollow glass spheres (SHGS)
| |
| with an average diameter of <math>d_{\text{avg,SHGS}}</math>=10~µm and
| |
| a density of <math>\rho_{SHGS}</math> = 1400 kg m<math>^{-3}</math> as tracer
| |
| particles. To prove the following behavior of these particles a
| |
| Stokes number Sk=1.08 and a particle sedimentation velocity
| |
| <math>u_{\text{SHGS}}=2.18 \times 10^{-5}~\text{m/s}</math> is calculated
| |
| With this Stokes number and a particle sedimentation velocity which is
| |
| much lower than the expected velocities in the experiments, an eminent
| |
| following behavior is approved. The camera takes 12 bit pictures with
| |
| a frequency of about 7.0 Hz and a resolution of 1695 x 1211 px with respect to the rectangular size of the test
| |
| section. For one phase-resolved position (described in
| |
| Section~\ref{sec:Generation_of_phase-resolved_data}) 60 to 80
| |
| measurements are taken. Preliminary studies with more and fewer
| |
| measurements showed that this number of measurements represent a good
| |
| compromise between accuracy and effort. The grid has a size of 150 x 138 cells and was calibrated with an average
| |
| factor of 126 <math>\mu</math> m/px}, covering a planar flow field of
| |
| x/D = -2.36 to 7.26 and y/D = -3.47 to ~3.47 in the middle of the test section at
| |
| z/D = 0. The time between the frame-straddled laser
| |
| pulses was set to <math>\Delta</math> t=200 <math>\mu</math>s. Laser and camera were
| |
| controlled by a TSI synchronizer (TSI 610035) with 1 ns
| |
| resolution. The processing of the phase-resolved fluid velocity fields
| |
| involving the structure deflections is described in
| |
| Section "Generation of phase-resolved data".
| |
| | |
| === Laser distance sensor ===
| |
| | |
| Non-contact structural measurements are often based on laser distance
| |
| techniques. In the present benchmark case the flexible structure shows
| |
| an oscillating frequency of about 7.1 Hz. With the
| |
| requirement to perform more than 100 measurements per period, a
| |
| time-resolved system was needed. Therefore, a laser triangulation was
| |
| chosen because of the known geometric dependencies, the high data
| |
| rates, the small measurement range and the resulting higher accuracy
| |
| in comparison with other techniques such as laser phase-shifting or
| |
| laser interferometry. The laser triangulation uses a laser beam which
| |
| is focused onto the object. A CCD-chip located near the laser output
| |
| detects the reflected light on the object surface. If the distance of
| |
| the object from the sensor changes, also the angle changes and thus
| |
| the position of its image on the CCD-chip. From this change in
| |
| position the distance to the object is calculated by simple
| |
| trigonometric functions and an internal length calibration adjusted to
| |
| the applied measurement range. To study simultaneously more than one
| |
| point on the structure, a multiple-point triangulation sensor was
| |
| applied (Micro-Epsilon scanControl 2750, see
| |
| Fig. 6). This sensor uses a matrix of CCD
| |
| chips to detect the displacements on up to 640 points along a laser
| |
| line reflected on the surface of the structure with a data rate of 800 profiles per second. The laser line was positioned in a
| |
| horizontal (x/D = 3.2, see
| |
| Fig. 6(a)) and in a vertical
| |
| alignment(z/D = 0, see
| |
| Fig. 6(b)) and has an accuracy of
| |
| 40 µm}. Due to the different refraction indices of air,
| |
| glass and water a custom calibration was performed to take the
| |
| modified optical behavior of the emitted laser beams into account.
| |
| | |
| [[File:structure_sensors_scancontrolonly0001_new.png]]
| |
| | |
| Fig. 6: Setup and alignment of multiple-point laser sensor on the flexible structure in a) z-direction and b) x-direction.
| |
| | |
| = Numerical Simulation Methodology =
| |
| | |
| The applied numerical method relies on an efficient partitioned
| |
| coupling scheme developed for dynamic fluid-structure interaction
| |
| problems in turbulent flows (Breuer et al, 2012). The fluid flow is
| |
| predicted by an eddy-resolving scheme, i.e., the large-eddy simulation
| |
| technique. FSI problems very often encounter instantaneous
| |
| non-equilibrium flows with large-scale flow structures such as
| |
| separation, reattachment and vortex shedding. For this kind of flows
| |
| the LES technique is obviously the best
| |
| choice (Breuer, 2002). Based on a semi-implicit scheme the LES code
| |
| is coupled with a code especially suited for the prediction of shells
| |
| and membranes. Thus an appropriate tool for the time-resolved
| |
| prediction of instantaneous turbulent flows around light, thin-walled
| |
| structures results. Since all details of this methodology were
| |
| recently published in Breuer et al, 2012, in the following only a
| |
| brief description is provided.
| |
| | |
| == Computational fluid dynamics (CFD) ==
| |
| | |
| Within a FSI application the computational domain is no longer fixed
| |
| but changes in time due to the fluid forces acting on the structure.
| |
| This temporally varying domain is taken into account by the Arbitrary
| |
| Lagrangian-Eulerian (ALE) formulation expressing the conservation
| |
| equations for time-dependent volumes and surfaces. Here the filtered
| |
| Navier-Stokes equations for an incompressible fluid are solved. Owing
| |
| to the deformation of the grid, extra fluxes appear in the governing
| |
| equations which are consistently determined considering the
| |
| \emph{space conservation law (SCL)}(Demirdzic 1988 and 1990, Lesoinne, 1996). The SCL is expressed
| |
| by the swept volumes of the corresponding cell faces and assures that
| |
| no space is lost during the movement of the grid lines. For this
| |
| purpose the in-house code FASTEST-3D (Durst et al, 1996a, b)
| |
| relying on a three-dimensional finite-volume scheme is used. The
| |
| discretization is done on a curvilinear, block-structured body-fitted
| |
| grid with collocated variable arrangement. A midpoint rule
| |
| approximation of second-order accuracy is used for the discretization
| |
| of the surface and volume integrals. Furthermore, the flow variables
| |
| are linearly interpolated to the cell faces leading to a second-order
| |
| accurate central scheme. In order to ensure the coupling of pressure
| |
| and velocity fields on non-staggered grids, the momentum interpolation
| |
| technique of Rhie (1983) is used.
| |
| | |
| A predictor-corrector scheme (projection method) of second-order
| |
| accuracy forms the kernel of the fluid solver. In the predictor step
| |
| an explicit three substep low-storage Runge-Kutta scheme advances
| |
| the momentum equation in time leading to intermediate
| |
| velocities. These velocities do not satisfy mass conservation. Thus,
| |
| in the following corrector step the mass conservation equation has to
| |
| be fulfilled by solving a Poisson equation for the
| |
| pressure-correction based on the incomplete LU decomposition method
| |
| of Stone (1968). The corrector step is repeated (about 3 to 8
| |
| iterations) until a predefined convergence criterion (<math>\Delta{m} <
| |
| {O}(10^{-9})</math>) is reached and the final velocities and the
| |
| pressure of the new time step are obtained. In Breuer et al (2012) it
| |
| is explained that the original pressure-correction scheme applied
| |
| for fixed grids has not to be changed concerning the mass conservation
| |
| equation in the context of moving grids. Solely in the momentum
| |
| equation the grid fluxes have to be taken into account as described
| |
| above.
| |
| | |
| In LES the large scales in the turbulent flow field are resolved
| |
| directly, whereas the non-resolvable small scales have to be taken
| |
| into account by a subgrid-scale model. Here the well-known and most
| |
| often used eddy-viscosity model, i.e., the Ssmagorinsky (1963) model
| |
| is applied. The filter width is directly coupled to the volume of the
| |
| computational cell and a Van Driest damping function ensures a
| |
| reduction of the subgrid length near solid walls. Owing to minor
| |
| influences of the subgrid-scale model at the moderate Reynolds number
| |
| considered in this study, a dynamic procedure to determine the
| |
| Smagorinsky parameter as suggested by Germano et al (1991) was omitted and
| |
| instead a well established standard constant <math>C_s = 0.1</math> is used.
| |
| | |
| The CFD prediction determines the forces on the structure and delivers
| |
| them to the CSD calculation. In the other direction the CSD prediction
| |
| determines displacements at the moving boundaries of the computational
| |
| domain for the fluid flow. The task is to adapt the grid of the inner
| |
| computational domain based on these displacements at the
| |
| interface. For moderate deformations algebraic methods are found to be
| |
| a good compromise since they are extremely fast and deliver reasonable
| |
| grid point distributions maintaining the required high grid
| |
| quality. Thus, the grid adjustment is performed based on a transfinite
| |
| interpolation (Thompson et al., 1985). It consists of three shear
| |
| transformations plus a tensor-product transformation.
| |
| | |
| == Computational structural dynamics (CSD) ==
| |
| | |
| The dynamic equilibrium of the structure is described by the momentum
| |
| equation given in a Lagrangian frame of reference. Large deformations,
| |
| where geometrical non-linearities are not negligible, are allowed
| |
| (Hojjat et al, 2010). According to preliminary structure considerations, a total Lagrangian
| |
| formulation in terms of the second Piola-Kirchhoff stress tensor and
| |
| the Green-Lagrange strain tensor which are linked by the
| |
| St. Venant-Kirchhoff material law is used in the present study.
| |
| For the solution of the governing equation the finite-element solver
| |
| Carat++, which was developed with an emphasis on the prediction of
| |
| shell or membrane behavior, is applied. Carat++ is based on several
| |
| finite-element types and advanced solution strategies for form finding
| |
| and non-linear dynamic
| |
| Problems (Wüchner et al, 2005; Bletzinger et al, 2005; Dieringer et al, 2012). For the
| |
| dynamic analysis, different time-integration schemes are available,
| |
| e.g., the implicit generalized-<math>\alpha</math> method (Chung et al, 1993). In the
| |
| modeling of thin-walled structures, membrane or shell elements are
| |
| applied for the discretization within the finite-element model. In the
| |
| current case, the deformable solid is modeled with a 7-parameter shell
| |
| element. Furthermore, special care is given to prevent locking
| |
| phenomena by applying the well-known Assumed Natural Strain (ANS) (Hughes and Tezduyar, 1981; Park and Stanley, 1986) and Enhanced Assumed Strain (EAS) methods (Bischoff et al., 2004).
| |
| | |
| Both, shell and membrane elements reflect geometrically reduced
| |
| structural models with a two-dimensional representation of the
| |
| mid-surface which can describe the three-dimensional physical
| |
| properties by introducing mechanical assumptions for the thickness
| |
| direction. Due to this reduced model additional operations are
| |
| required to transfer information between the two-dimensional structure
| |
| and the three-dimensional fluid model. Thus in the case of shells, the
| |
| surface of the interface is found by moving the two-dimensional
| |
| surface of the structure half of the thickness normal to the surface
| |
| on both sides and the closing of the volume (Bletzinger et al., 2006).. On
| |
| these two moved surfaces the exchange of data is performed
| |
| consistently with respect to the shell theory (Hojjat et al., 2010).
| |
| | |
| == Coupling algorithm ==
| |
| | |
| To preserve the advantages of the highly adapted CSD and CFD codes and
| |
| to realize an effective coupling algorithm, a partitioned but
| |
| nevertheless strong coupling approach is chosen. Since LES
| |
| typically requires small time steps to resolve the turbulent flow
| |
| field, the coupling scheme relies on the explicit
| |
| predictor-corrector scheme forming the kernel of the fluid solver.
| |
| | |
| Based on the velocity and pressure fields from the corrector step, the
| |
| fluid forces resulting from the pressure and the viscous shear
| |
| stresses at the interface between the fluid and the structure are
| |
| computed. These forces are transferred by a grid-to-grid data
| |
| interpolation to the CSD code Carat++ using a conservative
| |
| interpolation scheme (Farhat et al, 1998) implemented in the coupling
| |
| interface CoMA (Gallinger et al, 2009). Using the fluid forces provided via
| |
| CoMA, the finite-element code Carat++ determines the stresses in the
| |
| structure and the resulting displacements of the structure. This
| |
| response of the structure is transferred back to the fluid solver via
| |
| CoMA applying a bilinear interpolation which is a consistent scheme
| |
| for four-node elements with bilinear shape functions.
| |
| | |
| For moderate and high density ratios between the fluid and the
| |
| structure, e.g., a flexible structure in water, the added-mass effect
| |
| by the surrounding fluid plays a dominant role. In this situation a
| |
| strong coupling scheme taking the tight interaction between the fluid
| |
| and the structure into account, is indispensable. In the coupling
| |
| scheme developed in Breuer et al (2012) this issue is taken into
| |
| account by a FSI-subiteration loop which works as follows:
| |
| | |
| A new time step begins with an estimation of the displacement of the
| |
| structure. For the estimation a linear extrapolation is applied taking
| |
| the displacement values of two former time steps into account.
| |
| According to these estimated boundary values, the entire computational
| |
| grid has to be adapted as it is done in each FSI-subiteration
| |
| loop. Then the predictor-corrector scheme of the next time step is
| |
| carried out and the cycle of the FSI-subiteration loop is
| |
| entered. After each FSI-subiteration first the FSI convergence is
| |
| checked. Convergence is reached if the <math>L_2</math> norm of the displacement
| |
| differences between two FSI-subiterations normalized by the <math>L_2</math> norm
| |
| of the changes in the displacements between the current and the last
| |
| time step drops below a predefined limit, e.g. <math>\varepsilon_{FSI} =
| |
| 10^{-4}</math> for the present study. Typically, convergence is not reached
| |
| within the first step but requires a few FSI-subiterations (5 to
| |
| 10). Therefore, the procedure has to be continued on the fluid
| |
| side. Based on the displacements on the fluid-structure interface,
| |
| which are underrelaxated by a constant factor $\omega$ during the
| |
| transfer from the CFD to the CSD solver, the inner computational CFD
| |
| grid is adjusted. The key point of the coupling procedure suggested in
| |
| Breuer et al (2012) is that subsequently only the corrector step of
| |
| the predictor-corrector scheme is carried out again to obtain a new
| |
| velocity and pressure field. Thus the clue is that the pressure is
| |
| determined in such a manner that the mass conservation is finally
| |
| satisfied. Furthermore, this extension of the predictor-corrector
| |
| scheme assures that the pressure forces as the most relevant
| |
| contribution to the added-mass effect, are successively updated until
| |
| dynamic equilibrium is achieved. In conclusion, instabilities due to
| |
| the added-mass effect known from loose coupling schemes are avoided
| |
| and the explicit character of the time-stepping scheme beneficial for
| |
| LES is still maintained.
| |
| | |
| The code coupling tool CoMA is based on the
| |
| Message-Passing-Interface (MPI) and thus runs in parallel to the
| |
| fluid and structure solver. The communication in-between the codes is
| |
| performed via standard MPI commands. Since the parallelization in
| |
| FASTEST-3D and Carat++ also relies on MPI, a hierarchical
| |
| parallelization strategy with different levels of parallelism is
| |
| achieved. According to the CPU-time requirements of the different
| |
| subtasks, an appropriate number of processors can be assigned to the
| |
| fluid and the structure part. Owing to the reduced structural models
| |
| on the one side and the fully three-dimensional highly resolved fluid
| |
| prediction on the other side, the predominant portion of the CPU-time
| |
| is presently required for the CFD part. Additionally, the communication
| |
| time between the codes via CoMA and within the CFD solver takes a
| |
| non-negligible part of the computational resources.
| |
| | |
| == Numerical CFD Setup ==
| |
| | |
| For the CFD prediction of the flow two different block-structured
| |
| grids either for a subset of the entire channel (w'/l=1) or for
| |
| the full channel but without the gap between the flexible structure
| |
| and the side walls (w/l=2.95) are used. In the first
| |
| case the entire grid consists of about 13.5 million control volumes
| |
| (CVs), whereas 72 equidistant CVs are applied in the spanwise
| |
| direction. For the full geometry the grid possesses about 22.5 million
| |
| CVs. In this case starting close to both channel walls the grid is
| |
| stretched geometrically with a stretching factor 1.1 applying in
| |
| total 120 CVs with the first cell center positioned at a distance of
| |
| <math>\Delta</math>y/D=1.7 x 10<math>^{-2}</math>.
| |
| | |
| | |
| [[File:Benchmark_FSI-PfS-1_full_and_subset_case.jpg]]
| |
| | |
| Fig. 7: X-Y cross-section of the grid used for the simulation (Note that only every fourth grid line in each direction is displayed here).
| |
| | |
| | |
| The gap between the elastic structure and the walls is not taken into
| |
| account in the numerical model and thus the width of the channel is
| |
| set to w instead of W. The stretching factors are kept below 1.1 with the first cell
| |
| center located at a distance of <math>\Delta</math>y/D=9 x 10<math>^{-4}</math> from
| |
| the flexible structure. Based on the wall shear stresses on the
| |
| flexible structure the average y<math>^+</math> values are predicted to be below
| |
| 0.8, mostly even below 0.5. Thus, the viscous sublayer on the
| |
| elastic structure and the cylinder is adequately resolved. Since the
| |
| boundary layers at the upper and lower channel walls are not
| |
| considered, no grid clustering is required here.
| |
| | |
| On the CFD side no-slip boundary conditions are applied at the rigid
| |
| front cylinder and at the flexible structure. Since the resolution of
| |
| the boundary layers at the channel walls would require the bulk of the
| |
| CPU-time, the upper and lower channel walls are assumed to be slip
| |
| walls. Thus the blocking effect of the walls is maintained without
| |
| taking the boundary layers into account. At the inlet a constant
| |
| streamwise velocity is set as inflow condition without adding any
| |
| perturbations. The choice of zero turbulence level is based on the
| |
| consideration that such small perturbations imposed at the inlet will
| |
| generally not reach the cylinder due to the coarseness of the grid at
| |
| the outer boundaries. Therefore, all inflow fluctuations will be
| |
| highly damped. However, since the flow is assumed to be sub-critical,
| |
| this disregard is insignificant. At the outlet a convective outflow
| |
| boundary condition is favored allowing vortices to leave the
| |
| integration domain without significant disturbances
| |
| (Breuer, 2002). The convection velocity is set to
| |
| <math>u_\text{inflow}</math>.
| |
| | |
| As mentioned above two different cases are considered. In order to
| |
| save CPU-time in the first case only a subset of the entire spanwise
| |
| extension of the channel is taken into account. Thus the computational
| |
| domain has a width of w'/l=1 in z-direction and the
| |
| flexible structure is a square in the x-z-plane. In this case a
| |
| reasonable approximation already applied in Breuer et al (2012) is to
| |
| apply periodic boundary conditions in spanwise direction for both
| |
| disciplines. For LES predictions periodic boundary conditions
| |
| represent an often used measure in order to avoid the formulation of
| |
| appropriate inflow and outflow boundary conditions. The approximation
| |
| is valid as long as the turbulent flow is homogeneous in the specific
| |
| direction and the width of the domain is sufficiently large. The
| |
| latter can be proven by predicting two-point correlations, which have
| |
| to drop towards zero within the half-width of the domain. The impact
| |
| of periodic boundary conditions on the CSD predictions are discussed
| |
| below.
| |
| | |
| For the full case with w/l=2.95 periodic boundary conditions
| |
| can no longer be used. Instead, for the fluid flow similar to the
| |
| upper and lower walls also for the lateral boundaries slip walls are
| |
| assumed since the full resolution of the boundary layers would be
| |
| again too costly. Furthermore, the assumption of the slip wall is
| |
| consistent with the disregard of the small gap between the flexible
| |
| structure and the side walls discussed above.
| |
| | |
| == Numerical CSD Setup ==
| |
| | |
| Motivated by the fact that in the case of LES frequently a domain
| |
| modeling based on periodic boundary conditions at the lateral walls is
| |
| used to reduce the CPU-time requirements, this special approach was
| |
| also investigated for the FSI test case. The detailed discussion of
| |
| this specific boundary modeling for the spanwise direction is given in
| |
| Section~\ref{section:bc}. As a consequence, there are two different
| |
| structure meshes used: For the CSD prediction of the case with a
| |
| subset of the full channel the elastic structure is resolved by the
| |
| use of 10 x 10 quadrilateral four-node 7-parameter shell elements. For the
| |
| case discretizing the entire channel, 10 quadrilateral four-node shell
| |
| elements are used in the main flow direction and 30 in the spanwise
| |
| direction.
| |
| | |
| On the CSD side, the flexible shell is loaded on the top and bottom
| |
| surface by the fluid forces, which are transferred from the fluid mesh
| |
| to the structure mesh. These Neumann boundary conditions for the
| |
| structure reflect the coupling conditions. Concerning the Dirichlet
| |
| boundary conditions, the four edges need appropriate support modeling:
| |
| on the upstream side at the rigid cylinder a clamped support is
| |
| realized and all degrees of freedom are equal to zero. On the opposite
| |
| downstream trailing-edge side, the rubber plate is free to move and
| |
| all nodes have the full set of six degrees of freedom. The edges which
| |
| are aligned to the main flow direction need different boundary
| |
| condition modeling, depending on whether the subset or the full case
| |
| is computed:
| |
| | |
| For the subset case due to the fluid-motivated periodic boundary
| |
| conditions, periodicity for the structure is correspondingly assumed
| |
| for consistency reasons. As it turns out, this assumption seems to
| |
| hold for this specific benchmark configuration and its deformation
| |
| pattern which has strong similarity with an oscillation in the first
| |
| eigenmode of the plate. Hence, this modeling approach may be used for
| |
| the efficient processing of parameter studies, e.g., to evaluate the
| |
| sensitivity of the FSI simulations with respect to slight variations
| |
| in model parameters shown in a sensitivity study. For
| |
| this special type of support modeling, there are always two structure
| |
| nodes on the lateral sides (one in a plane z=-w/2 and its twin in
| |
| the other plane z=+w/2) which have the same load. These two nodes
| |
| must have the same displacements in x- and y-direction and their
| |
| rotations have to be identical. Moreover, the periodic boundary
| |
| conditions imply that the z-displacement of the nodes on the sides are
| |
| forced to be zero.
| |
| | |
| For the full case the presence of the walls in connection with the
| |
| small gap implies that there is in fact no constraining effect on the
| |
| structure, as long as no contact between the plate and the wall takes
| |
| place. Out of precise observations in the lab, the possibility of
| |
| contact may be disregarded. In principle, this configuration would
| |
| lead to free-edge conditions like at the trailing edge. However, the
| |
| simulation of the fluid with a moving mesh needs a well-defined mesh
| |
| situation at the side walls which made it necessary to tightly connect
| |
| the structure mesh to the walls (the detailed representation of the
| |
| side edges within the fluid mesh is discarded due to computational
| |
| costs and the resulting deformation sensitivity of the mesh in these
| |
| regions). Also the displacement in z-direction of the structure nodes
| |
| at the lateral boundaries is forced to be zero.
| |
| | |
| == Coupling conditions ==
| |
| | |
| For the turbulent flow a time-step size of <math>\Delta t_{f} = 2 \times
| |
| 10^{-5}s~(\Delta t_{f}^{\ast} = 1.26 \times 10^{-3})</math> in
| |
| dimensionless form using <math>u_\text{inflow}</math> and D as reference
| |
| quantities) is chosen and the same time-step size is applied for the
| |
| structural solver based on the generalized-<math>\alpha</math> method with the
| |
| spectral radius <math>\varrho_\infty</math>=1.0, i.e, the Newmark standard
| |
| method. For the CFD part this time-step size corresponds to a CFL
| |
| number in the order of unity. Furthermore, a constant underrelaxation
| |
| factor of <math>\omega</math>=0.5 is considered for the displacements and the
| |
| loads are transferred without underrelaxation. In accordance with
| |
| previous laminar and turbulent cases in Breuer et. al. (2012) the FSI
| |
| convergence criterion is set to <math>\varepsilon_{FSI} = 10^{-4}</math> for the
| |
| <math>L_2</math> norm of the displacement differences. As estimated from previous
| |
| cases Breuer et. al. (2012) 5 to 10 FSI-subiterations are required to
| |
| reach the convergence criterion.
| |
| | |
| After an initial phase in which the coupled system reaches a
| |
| statistically steady state, each simulation is carried out for about
| |
| 4 s real-time corresponding to about 27 swiveling cycles of the
| |
| flexible structure.
| |
| | |
| For the coupled LES predictions the national supercomputer
| |
| SuperMIG/SuperMUC was used applying either 82 or 140 processors for
| |
| the CFD part of the reduced and full geometry,
| |
| respectively. Additionally, one processor is required for the coupling
| |
| code and one processor for the CSD code, respectively.
| |
| | |
| = Unsteady results =
| |
| | |
| In order to comprehend the real structure deformation and the
| |
| turbulent flow field found in the present test case, experimentally
| |
| and numerically obtained unsteady results are presented in this
| |
| section.
| |
| | |
| A high-speed camera movie of the structure deflection is available at: http://vimeo.com/59130974.
| |
| | |
| Figure 8 shows experimental raw signals of dimensionless displacements from a point located at a
| |
| distance of 9 mm from the shell extremity in the midplane of the test
| |
| section. In Figure 8a) the history
| |
| of the y-displacement <math>U_y^* = U_y / D</math> obtained in the
| |
| experiment is plotted. The signal shows significant variations in the
| |
| extrema: The maxima of <math>U_y^*</math> vary between 0.298 and 0.523 and the
| |
| minima between -0.234 and -0.542. The standard deviations on the
| |
| extrema are about <math>\pm 0.05~(\pm 12 \%)</math> of the mean value of the
| |
| extrema). Minor variations are observed regarding the period in
| |
| Figure 8a). Figure 8b)
| |
| and 8c) show the corresponding
| |
| experimental phase portrait and phase plane, respectively. The phase
| |
| portrait has a quasi-ellipsoidal form. The monitoring point trajectory
| |
| plotted in the phase plane describes an inversed 'C', which is
| |
| typical for the first swiveling mode. The cycle-to-cycle variations in
| |
| these plots are small. Therefore, the FSI phenomenon can be
| |
| characterized as quasi-periodic.
| |
| | |
| [[File:unsteady_LDT.png]]
| |
| | |
| Fig. 8: Experimental raw signals of dimensionless displacements from a point in the midplane of the test section located at a distance of 9 mm from the shell extremity.
| |
| | |
| | |
| Figure 9 is composed of eight
| |
| images of the instantaneous flow field (streamwise velocity component)
| |
| experimentally measured in the x-y plane located in the middle of the
| |
| rubber plate. These pictures constitute a full period T of the FSI
| |
| phenomenon arbitrarily chosen. As mentioned before, the shell deforms
| |
| in the first swiveling mode. Thus, there is only one wave node located
| |
| at the clamping of the flexible structure. At the beginning of the
| |
| period (t = 0) the structure is in its undeformed state. Then, it
| |
| starts to deform upwards and reaches a maximal deflection at t
| |
| = T / 4. Afterwards, the shell deflects downwards until its
| |
| maximal deformation at t =3T/4. Finally the plate deforms
| |
| back to its original undeformed state and the end of the period is
| |
| reached.
| |
| | |
| As visible in Fig. 9 the flow is
| |
| highly turbulent, particularly near the cylinder, the flexible
| |
| structure and in the wake. The strong shear layers originating from
| |
| the separated boundary layers are clearly visible. This is the region
| |
| where for the sub-critical flow the transition to turbulence takes
| |
| place as visible in
| |
| Fig. 9. Consequently, the flow in
| |
| the wake region behind the cylinder is obviously turbulent and shows
| |
| cycle-to-cycle variations. That means the flow field in the next
| |
| periods succeeding the interval depicted in
| |
| Fig. 9 will definitely look
| |
| slightly different due to the irregular chaotic character of
| |
| turbulence. Therefore, in order to be able to compare these results an
| |
| averaging method is needed leading to a statistically averaged
| |
| representation of the flow field. Since the FSI phenomenon is
| |
| quasi-periodic the phase-averaging procedure presented above is ideal
| |
| for this purpose and the results obtained are presented in the next
| |
| section.
| |
| | |
| | |
| [[File:unsteady_PIV.png]]
| |
| | |
| Fig. 9: Experimental unsteady flow field (x-y plane located in the middle of the rubber plate).
| |
| | |
| Prior to this, however, it should be pointed out that very similar
| |
| figures as depicted in Fig. 9 could
| |
| also be shown from the numerical predictions based on LES. Exemplary
| |
| and for the sake of brevity, Fig. 10 displays the
| |
| streamwise velocity component of the flow field in a x-y-plane solely
| |
| at t=3T/4. As expected the LES prediction is capable to
| |
| resolve small-scale flow structures in the wake region and in the
| |
| shear layers. Furthermore, the figure visualizes the deformed
| |
| structure showing nearly no variation in spanwise direction.
| |
| | |
| [[File:unsteady_LES.png]]
| |
| | |
| Fig. 10: Experimental unsteady flow field (x-y plane located in the middle of the rubber plate).
| |
| | |
| = Generation of Phase-resolved Data =
| |
| | |
| Each flow characteristic of a quasi-periodic FSI problem can be
| |
| written as a function <math>f=\bar{f}+\tilde{f}+f'</math>, where <math>\bar{f}</math>
| |
| describes the global mean part, <math>\tilde{f}</math> the quasi-periodic part
| |
| and <math>f'</math> a random turbulence-related
| |
| part (Reynolds et. al., 1972; Cantwell et. al., 1983). This splitting can also be written
| |
| in the form <math>f = <f> + f'</math>, where <math><f></math> is the
| |
| phase-averaged part, i.e., the mean at constant phase. In order to be
| |
| able to compare numerical results and experimental measurements, the
| |
| irregular turbulent part f' has to be averaged out. This measure is
| |
| indispensable owing to the nature of turbulence which solely allows
| |
| reasonable comparisons based on statistical data. Therefore, the
| |
| present data are phase-averaged to obtain only the phase-resolved
| |
| contribution <math><f></math> of the problem, which can be seen as a
| |
| representative and thus characteristic signal of the underlying FSI
| |
| phenomenon.
| |
| | |
| The procedure to generate phase-resolved results is the same for the
| |
| experiments and the simulations and is also similar to the one
| |
| presented in Gomes et. al. (2006). The technique can be split up into
| |
| three steps:
| |
| | |
| * '''Reduce the 3D-problem to a 2D-problem''' - Due to the facts that in the present benchmark the structure deformation in spanwise direction is negligible and that the delivered experimental PIV-results are solely available in one x-y-plane, first the 3D-problem is reduced to a 2D-problem. For this purpose the flow field and the shell position in the CFD predictions are averaged in spanwise direction.
| |
| | |
| * '''Determine n reference positions for the FSI Problem''' - A representative signal of the FSI phenomenon is the history of the y-displacements of the shell extremity. Therefore, it is used as the trigger signal for this averaging method leading to phase-resolved data. Note that the averaged period of this signal is denoted T. At first, it has to be defined in how many sub-parts the main period of the FSI problem will be divided and so, how many reference positions have to be calculated (for example in the present work n = 23). Then, the margins of each period of the y-displacement curve are determined. In order to do that the intersections between the y-displacement curve and the zero crossings (<math>U_y</math>=0) are looked for and used to limit the periods. Third, each period <math>T_i</math> found is divided into n equidistant sub-parts denoted j.
| |
| | |
| * '''Sort and average the data corresponding to each reference Position''' - The sub-part j of the period <math>T_i</math> corresponds to the sub-part j of the period <math>T_{i+1}</math> and so on. Each data set found in a sub-part j will be averaged with the other sets found in the sub-parts j of all other periods. Finally, data sets of n phase-averaged positions for the representative reference period are achieved.
| |
| | |
| The simulation data containing structure positions, pressure and
| |
| velocity fields, are generated every 150 time steps. According to the
| |
| frequency observed for the structure and the time-step size chosen
| |
| about 50 data sets are obtained per swiveling period. With respect to
| |
| the time interval predicted and the number of subparts chosen, the
| |
| data for each subpart are averaged from about 50 data sets. A
| |
| post-processing program is implemented based on the method described
| |
| above. It does not require any special treatment and thus the
| |
| aforementioned method to get the phase-resolved results is
| |
| straightforward.
| |
| | |
| The current experimental setup consists of the multiple-point triangulation
| |
| sensor described in Section "Laser distance sensor" and the
| |
| synchronizer of the PIV system. Each measurement pulse of the PIV
| |
| system is detected in the data acquisition of the laser distance
| |
| sensor, which measures the structure deflection continuously with 800
| |
| profiles per second. With this setup, contrary to Gomes et. al. (2006),
| |
| the periods are not detected during the acquisition but in the
| |
| post-processing phase. After the run a specific software based on the
| |
| described method mentioned above computes the reference structure
| |
| motion period and sorts the PIV data to get the phase-averaged
| |
| results. A more detailed description of the phase-averaging method is available in De Nayer et. al. (2013).
| |
| ---- | | ---- |
|
| |
|