Grid Convergence Property of Three-Dimensional Measurement-Integrated Simulation for Unsteady Flow behind a Square Cylinder with Karman Vortex Street

The purpose of this study was to clarify grid convergence property of three-dimensional measurement-integrated (3D-MI) simulation for a flow behind a square cylinder with Karman vortex street. Measurement-integrated (MI) simulation is a kind of the observer in the dynamical system theory by using CFD scheme as a mathematical model of the system. In a former study, two-dimensional MI (2D-MI) simulation with a coarse grid system showed a fairly good result in comparison with a 2D ordinary (2D-O) simulation, but the results were degraded with grid refinement. In this study, 3D-MI simulation and three-dimensional ordinary (3D-O) simulation were performed with three grid systems of different grid resolutions, and their grid convergence properties were compared. As a result, all 3D-MI simulations reproduced the vortex shedding frequency identical to that of the experiment, and the flow fields obtained were very close, within 5% difference between the results, while the results of the 3D-O simulations showed variation of the solution under convergence. It is shown that the grid convergence property of 3D-MI simulation is monotonic and better than that of 3D-O simulation, whereas those of 2D-O and 2D-MI simulations for streamwise velocity fluctuation are divergent. The solution of 3D-MI simulation with a relatively coarse grid system properly reproduces the basic three-dimensional structure of the wake flow as well as the drag and lift coefficients.


Introduction
Both numerical simulation and measurement are essential techniques in analyses of flow problems.Flow measurement provides the actual property of a flow such as velocity and pressure within the range of measurement precision, but they are usually limited to partial information of the flow.Numerical simulation, therefore, plays an important role in elucidating the detailed structure of a flow.The numerical simulation yields detailed velocity and pressure fields of the flow, but its reliability depends on boundary and initial conditions as well as unknown disturbances.Therefore, neither numerical simulation nor measurement is expected to reproduce "real flows" exactly.
In order to reproduce real flows correctly, many studies have been conducted by integrating numerical simulation and measurement [1].Combining particle tracking velocimetry (PTV) and direct numerical simulation (DNS) with a linear combination, Suzuki et al. have developed a method that obtains unmeasurable quantities of an unsteady flow such as pressure fields and vorticity distributions [2], and have evaluated the data-assimilation capabilities of the method [3].Buffoni et al. have proposed an estimation method of a velocity field based on a non-linear low-dimensional model of the flow with velocity or shear stress measurement, treating an unsteady flow around a square cylinder with a low Reynolds number [4].Recently, the methodology combining measurement and simulation has been used to obtain the velocity in regions where no velocity information is available in particle image velocimetry (PIV) [5], and to estimate time-resolved velocity fields from non-time-resolved PIV data [6].The present authors have also studied measurement-integrated (MI) simulation, which is a kind of the state observer in dynamical system theory employing the CFD scheme as a model instead of a linear differential equation [7] [8].The validity of the MI simulation has been shown for several flow related problems.Conducting a numerical experiment of a developed turbulent flow in a square duct, Hayase and Hayashi [7] have shown that the MI simulation converges to the model solution of the real flow by adding pressure difference into the numerical simulation as feedback signals.As a medical application, Funamoto et al. [9] conducted the MI simulation of a blood flow using Doppler velocity measurement of medical ultrasonography, showing the complex blood flow velocity field in an aneurism.
The flow treated herein is a wake flow behind a square cylinder in a confined wind tunnel.Flows around bluff bodies provide fundamental problems in many engineering applications such as basic flow for vehicles and large structures, including important phenomena such as vortex shedding.A number of studies on the flow around cylinders have been conducted numerically and experimentally.The critical Reynolds number has been investigated in flows around circular [10]- [12] and rectangular [13]- [15] cylinders in which vortex shedding from two-dimensional bluff bodies reveals three-dimensional behavior.Turbulent wake flow behind a bluff body is also the subject of intense scrutiny and has been studied experimentally and numerically [16].Recently, large-scale computations such as LES of a turbulent wake [17] [18] and simulation of three-dimensional vortical structures behind two square cylinders [19] [20] have been performed.Traditional studies gave priority to clarifying the relation between some flow parameters and flow phenomena, investigating the Strouhal-Reynolds number re-lation experimentally and numerically in a wide range of Reynolds number [21].Davis et al. [22] and Sohankar et al. [23] have reported that the Strouhal number is influenced by many parameters, i.e., the aspect ratio of a rectangular cylinder, blockage and inlet/outlet boundary conditions, as well as Reynolds number.As mentioned above, it is difficult to obtain exact boundary and initial conditions of a real flow.Therefore, predicting the shedding frequency of the Karman vortex, which is fundamental characteristic of the wake flow, is difficult if based only on numerical simulation under noncontrolled flow conditions.In contrast, the MI simulation reproduced the vortex shedding frequency and the velocity fluctuation of the experimental flow by using measured pressure under unknown boundary conditions in a previous study [24].In a two-dimensional analysis of a Karman vortex street behind a square cylinder, Nisugi et al. [25] developed a hybrid wind tunnel, an on-line flow analysis system, by integrating numerical simulation and wind tunnel experiment based on the MI simulation, and showed that the MI simulation reproduces the velocity fluctuation and the shedding frequency of the Karman vortex in the case of a relatively low Reynolds number of 1200, even with a coarse grid system.In addition, allowing us to use a coarse grid system, the MI simulation has been able to reproduce a real flow in real time.The reproduced pressure fields have shown good agreement with measured pressure on the sidewall of a wind tunnel [26].If the reproduction of the detailed information of real flows becomes possible in real time, it would be expected that the control of the flow behind bluff bodies [27] [28] can be performed with high accuracy.In a recent paper, we have reported that the increase of feedback data is effective to improve the reproducibility of the velocity field by integrating numerical simulation and particle image velocimetry [29].However, previous studies of the MI simulation in a hybrid wind tunnel were performed in two-dimensional analysis and the reproducibility of three-dimensional flow was not investigated.In addition, the MI simulation with a coarse grid system showed good agreement with the experiment, but the simulation with a fine grid system resulted in an increase of the error of the velocity fluctuation in a previous two-dimensional study [25].Therefore, the influence of grid resolutions on the reproducibility of a flow needs to be investigated in three-dimensional analysis of the MI simulation.
The purpose of this study was to clarify the grid convergence property of a threedimensional measurement-integrated (3D-MI) simulation for a flow behind a square cylinder with a Karman vortex street for a relatively low Reynolds number.3D-MI simulation and three-dimensional ordinary (3D-O) simulation were performed and their grid convergence properties were investigated by comparing with those of laser Doppler velocimetry (LDV) measurement and previous two-dimensional (2D) analysis.Some results were discussed from fluid dynamic point view for instantaneous 3D flow structure and variations of lift and drag coefficients.

Experimental Apparatus
A hybrid wind tunnel is an on-line flow analysis system, which integrates numerical simulation and flow measurement.Figure 1(a) shows the configuration of the system, including a process flowchart of the MI simulation.The system consists of a wind tunnel with a square cylinder, pressure sensors (SSK, DP8A-2, 20 Pa, resonant frequency 260 Hz), a workstation (SGI, OCTAIN, 300 MHz, 2 CPU) and a supercomputer (SGI, Silicon Graphics Prism, 1.6 MHz, 128 CPU).The appearance of the experimental apparatus is shown in Figure 1  (5 mm in mesh size, 1 mm in thickness) are placed at 205 mm (upstream side) and 105 mm (downstream side) from both ends of the wind tunnel.The flow is made by suction with an air blower (Nishimuradenki, NK-200), which is connected to the outflow boundary of the wind tunnel through a flexible tube (inner diameter of 175 mm).The flow rate of the air blower is controlled by an inverter (Mitsubishidenki, FR-520-1.5K).
Pressure at the front and the both sides of the cylinder is measured to provide the information of the real flow for the numerical simulation (see Figure 3).In a previous study, the MI simulation was performed in two-dimensional analysis and the operation described above was performed in real time [8].In the present study, we dealt with three-dimensional analysis, and the MI simulation was performed off-line with pre-measured pressure because of increased computational time.

Measurement-Integrated Simulation
Governing equations of MI simulation are the Navier-Stokes equation with a body force term (1) and the pressure Equation ( 2) for incompressible fluid.( ) The difference between MI simulation and ordinary simulation is the existence of the artificial body force f, which is a feedback signal derived from the difference between computed result and measurement data.Determination of the feedback signal is the key for the MI simulation design.
Since this study focused on the Karman vortex street, the measurement data are defined as the pressure on the both sides of the cylinder relative to the stagnation pressure (see Figure 3).Feedback forces F A and F B are added in the streamwise momentum equation at the control volumes in the upstream side of the pressure nodes.The forces are calculated proportional to the difference between the measured pressures AS P * , BS P * and computed ones P AS , P BS as follows [26]: where K* is the feedback gain (non-dimensional); ΔV is a control volume and L is a reference length (cylinder width D = 30 mm).Streamwise flows at the feedback points are accelerated or decelerated by the artificial forces F A , F B to reduce the difference between the measured pressure and the calculated one.It should be noted that the feedback forces are virtual forces to reduce the error of the numerical simulation.In previous two-dimensional analysis, the middle plane of the wind tunnel was treated, and pressures on the cylinder were measured in the plane for calculation of feedback forces.
In present study, the feedback forces calculated in the middle plane were applied to all feedback points along span-wise direction corresponding to the feedback law of the two-dimensional analysis (see Figure 3).

Computational Conditions
The governing equations were discretized with the finite volume method and solved with an algorithm similar to the SIMPLER method.The computational domain considered for analysis is shown in grid systems are summarized in Table 1.The relative CPU time in the table means the computational time relative to the one of the coarse grid A. The CPU time for the fine grid C is 55 times that for the coarse grid A with twice the number of CPUs.

Results and Discussion
To evaluate the grid convergence of the 3D-MI simulation, 3D-MI and 3D-O simulations were performed with the above-mentioned three grid systems.Figure 4(a) compares the time histories of pressure at the feedback point A (see Figure 3) between the 3D-O simulations and measurement.The measured pressure oscillates periodically corresponding to the vortex shedding from the square cylinder.The ordinary simulation with the coarse grid A and middle grid B fails to reproduce the pressure oscillation because the grid resolutions are insufficient for the ordinary simulation.The pressure oscillation of fine grid C shows nearly the same amplitude as that of the measurement, but is naturally out of phase with that of the measurement.The frequencies of the vortex shedding are compared among computations and experiment.Figure 5    the shedding frequency of the Karman vortex is affected by the upstream velocity profile.On the other hand, both 2D-and 3D-MI simulations reproduce the Strouhal number exactly identical to that of the experiment in all the grid systems.
Grid convergence was also investigated for distributions of the mean streamwise velocity and its fluctuation in comparison with the laser Doppler velocimetry (LDV) measurement in former work [24].LDV equipment (TSI, LDP 100) consists of an opti-   To evaluate the reproducibility of the velocity field quantitatively, we define the error norm as follows.
( ) ( ) where m u and m u * are time averaged values of the computation and LDV measurement for streamwise velocity or its fluctuation; i is grid index; and N is the number of the data in the domain.Figure 9(a) shows the error norm of the mean streamwise velocity u mean with grid spacing including the results of the former two-dimensional analysis   for the MI simulation with fine grid C. The contours of u and v velocities in the above figure represent the existence of the vortex shedding.In the figure of w, the flows heading to the middle of the wind tunnel appear periodically from connections between the cylinder and wind tunnel walls.Figure 11(b) shows the corresponding results of the 3D-MI simulation with coarse grid A. In the upper figures, large structures of u and v velocity components are well reproduced although the w velocity component is rather weak.In the lower figures, the results of u and v velocity components are symmetric on the centerline, while that of w is anti-symmetric.The distribution of u velocity is qualitatively the same as that of fine grid C although the recirculation area behind the cylinder is somewhat larger.Large structures of the v and w velocity distributions were also reproduced properly.Difference between the results of two grids seems more significant for w velocity in comparison with those of the other velocity components probably because the magnitude of the w velocity is smaller than those of the others.
Figure 12 shows time histories of drag and lift coefficients for the square cylinder obtained from the coarse gird 3D-MI and 3D-O simulations to be compared with that of the fine grid 3D-MI simulation.In Figure 12(a) the drag coefficient of the coarse grid 3D-O simulation is 1.6 and that of the coarse grid 3D-MI simulation is 2.2 in time average, which is close to that of the fine grid 3D-MI simulation of 2.4.Note that the drag coefficient of the fine grid 3D-MI simulation is slightly higher than the experimental value of 2.1 [22], probably due to the difference of the inlet boundary condition and blockage ratio.The coarse grid 3D-MI simulation, however, fails to reproduce the fluctuation of the drag.This is probably attributable to the fact that the coarse grid 3D-MI simulation cannot reproduce the detachment of the recirculation area from the cylinder (see Figure 11

Conclusion
The purpose of this study was to clarify grid convergence property of 3D-MI simulation  for a flow behind a square cylinder with Karman vortex street.3D-MI and 3D-O simulations were performed with three grid systems of different grid resolutions.All 3D-MI simulations reproduced the vortex shedding frequency exactly identical to that of the experiment, and the flow fields obtained were very close, within a difference of 5% between the results, while the results of the 3D-O simulations showed variation of the solution under grid convergence.It is shown that the grid convergence property of 3D-MI simulation is monotonic and better than that of 3D-O simulation, whereas those of 2D-O and 2D-MI simulations for streamwise velocity fluctuation are divergent.The solution of 3D-MI simulation with relatively coarse grid system properly reproduces the basic three-dimensional structure of the wake flow as well as the drag and lift coefficients.In this study, the relatively low Reynolds number of 1200 was treated in order to investigate fundamental grid convergence property of the 3D-MI simulation in comparison with the former 2D-MI simulation and experiment.Fully turbulent flow at higher Reynolds number is an important issue in practical applications, which is our future work.

Figure 2 .
Figure 2. Geometry of wind tunnel and computational domain (units in mm).
Pressure holes are located at the middle of the surfaces and are connected with the pressure sensors through a copper tube with a diameter of 3 mm and a length of 200 mm and a plastic tube with a diameter of 8 mm and a length of 200 mm.The resonant frequency of the tube is estimated to be 430 Hz under the assumption of a closed air column.Output signals of the pressure sensors are imported into the workstation with an analog-digital converter (VMIC, MVIPCI-3322, 1000 Hz, 16 bit) through a low-pass filter (NEC SANEI, 9B02, 10 Hz) to extract signals corresponding to vortex shedding.In the workstation, electric noises of the signals are suppressed by a digital low-pass filter (IIR-filter, Butterworth, 20 Hz) and the signals are transferred to a supercomputer with socket communication.The uncertainty in the pressure measurement was estimated to be 6%.Major error comes from sensor drift and electrical noise.Using the pressure signals, the MI simulation is performed on the supercomputer to obtain a solution corresponding to the real flow.Details of the MI simulation are described in the next section.The result of the MI simulation is transferred to the workstation and the flow field is visualized.

Figure 3 .
Figure 3. Details of the square cylinder.

Figure 2 .
Dimensions related to the domain are, respectively, 37.0D in the streamwise direction, 6.7D in the wall-normal and the span-wise directions.Three staggered grid systems in the Cartesian coordinate were employed.The numbers of grid points were, respectively, N x × N y × N z = 111 × 21 × 21 (grid A), 222 × 42 × 42 (grid B) and 444 × 84 × 84 (grid C).In a coarse grid system (grid A), the grid spacing was defined as h x = h y = h z = D/3, and the non-uniform meshes were partially employed near the sidewall of the wind tunnel in the wall-normal and span-wise directions to match the position of the sidewall.Although the present grid system is too coarse for ordinary simulation, the two-dimensional analysis in the previous study revealed that the velocity field of the Karman vortex street was correctly reproduced with this grid system[24].The uniform velocity profile with U = 0.605 m/s, corresponding to the Reynolds number Re = 1200 of previous studies[24]-[26], was applied at the inflow boundary.The convective boundary condition was used at the outflow boundary.The time step of computation was set as Δt = 0.01 s (non-dimensional time step 0.2).The frequency of the Karman vortex shedding was f s = 2.81 Hz at the Reynolds number of 1200 in wind tunnel experiment.The optimal feedback gain was so determined that the calculated pressure at the feedback point best agrees with the measured one.The gains determined by test computations were K* = 5.4 for grid A, K* = 10.8 for grid B and K* = 21.6 for grid C. The case of null feedback gain K* = 0.0 is equivalent to the ordinary numerical simulation without the feedback mechanism.Computational conditions for all

Figure 4 (
b) shows the corresponding results for the 3D-MI simulations.All MI simulations show in-phase oscillation with the measurement and good agreement in the amplitude.The effects of the feedback forces are especially seen in the results of the coarse grid A. Large feedback force appears in coarse grid A due to a large control volume, whereas the feedback gain is small.
K = 5.4 Grid B, K = 10.8Grid C, K = 21.6 cal head (diode laser power 50 mW, wavelength 690 nm, beam spacing 50 mm, focal length 250 mm) and a digital signal processor (Doppler frequency range 100 Hz to 75 MHz, downmix frequency 40 MHz).Oil mist with an average diameter of 1 μm was supplied by a fog generator (DANTEC, 2010).The number of valid data obtained in LDV processing was about 500 samples per second.In total, 413 measurement points were arranged in the middle plane of the wind tunnel, and the streamwise velocity was measured around the cylinder −3.5D ≤ x ≤ 4.5D at an interval of 0.33D and downstream 4.5D ≤ x ≤ 14.5D at an interval of 1.0D.Figure6(a)shows the distribution of the mean streamwise velocity with a recirculation area behind a square cylinder.Figure6(b)shows the distribution of the streamwise velocity fluctuation with two peaks due to the vortex shedding behind the cylinder.Both distributions show asymmetry due to non-uniformity of the inflow velocity distribution which was generated through the filter and the steel mesh mentioned in Section 2.1.

Figure 7
Figure 7 compares distributions of the mean streamwise velocity u m between 3D-O simulations (left) and 3D-MI simulations (right) with coarse grid A, middle grid B and fine grid C. The mean streamwise velocities u m are calculated for 20 s for grids A and B, and 10 s for grid C in steady oscillation, and are normalized by averaged inflow velocity U. Ordinary simulations with grid A and grid B show recirculation areas extending downstream to 6D and 4D, respectively, due to the poor grid resolution.The mean velocity field with grid C agrees with the experiment of Figure 6(a).In contrast, the MI simulations with grid A and grid B significantly improve the size of the recirculation area in comparison with the ordinary simulation.The result of MI simulation with grid C is almost the same as that of the ordinary simulation.

Figure 8 Figure 6 .Figure 7 .
Figure 8 compares distributions of the root-mean-square values of the streamwise velocity fluctuations u rms between 3D-O simulations (left) and 3D-MI simulations (right) with the three grid systems.These values are evaluated in the time intervals and
Figure 11.Distribution of the velocity components u (left), v (middle) and w (right) of the MI simulation on x-y plane (upper) and x-z plane (lower): (a) grid C; (b) grid A.
Figure12shows time histories of drag and lift coefficients for the square cylinder obtained from the coarse gird 3D-MI and 3D-O simulations to be compared with that of the fine grid 3D-MI simulation.In Figure12(a) the drag coefficient of the coarse grid 3D-O simulation is 1.6 and that of the coarse grid 3D-MI simulation is 2.2 in time average, which is close to that of the fine grid 3D-MI simulation of 2.4.Note that the drag coefficient of the fine grid 3D-MI simulation is slightly higher than the experimental value of 2.1[22], probably due to the difference of the inlet boundary condition and blockage ratio.The coarse grid 3D-MI simulation, however, fails to reproduce the fluctuation of the drag.This is probably attributable to the fact that the coarse grid 3D-MI simulation cannot reproduce the detachment of the recirculation area from the cylinder (see Figure11(a) and Figure11(b)).In Figure12(b), the coarse grid 3D-O simulation does not reproduce the oscillation of the lift coefficient, but the result of the coarse grid 3D-MI simulation agrees with that of the fine grid MI simulation.