Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
Numerical Solution of Viscous Flows about Submerged and Partly Submerged Bodies P. G. Esposito Istituto Nazionale di Studi ed Esperienze di Architettura Navale Roma, Italy G. Grazian~ and P. Orlandi University di "La Sapienza" Rome, Italy Abstract In this paper a finite difference method based on the Navier-Stokes equations in generalized curvilinear co- ordinates is applied to the simulation of unsteady vis- cous flows with free surface. An accurate analysis of the method is performed by comparing numerical and experimental results for several physical cases related to ship hydrodynamics. ~ Introduction The flow of a viscous fluid past a body shows a very complex behaviour when a free surface is involved. In fact the free surface strongly affects the flow field. The detailed analysis of the interactions between body and fluid is very important in the field of naval hydro- dynamics in order to quantify the viscous and wave resistance. Experimental investigation can provide global information at reasonable cost, while velocity and pressure measurement are difficult or even impos- sible and very expensive, instead numerical simulation gives the complete flow field quantities provided par- ticular care is devoted to implement the numerical scheme. The validation of the numerical solution by a comparison with the experimental findings is nec- essary to verify the accuracy of scheme. The best results, however, are obtained when the two fields of analysis compenetrate each other. In this paper a numerical method for the solution of the Navier-Stokes equations in general curvilin- ear coordinates based on an implicit finite difference scheme is presented. The physical domain bounded by the moving free surface is mapped onto a fixed com- putational domain. The time dependent coordinate transformation introduces extra terms that account for the grid velocity. For a complete test of the nu- merical method, the scheme has been applied to dif- ferent physical problems with increasing complexity. As a first case the 2-D flow field past a wall-mounted 607 square obstacle in a channel has been considered [1~. When the geometry is mapped onto a computational domain, the coordinate lines turn out to be highly distorted. The comparison with experimental results suggests that very fine meshes are required for an ac- curate resolution of flow reversal. To verify the treatment of the time dependent ge- ometry, the flow field inside a rectangular enclosure with a moving indentation has been studied and the temporal evolution has been compared with the ex- perimental findings. The 2-D free surface flow over a semicircular bump on a straight wall has been considered because it is closely related to ship design. This case is relevant for the presence of large recirculating regions and of the free surface. Experimental results are also available [2~. 2 Governing Equations The simulation of a two dimensional incompressible viscous flow has been considered as a first step to build a numerical scheme for 3-D flows. The gov- erning equations are: , . . ~ ovi Ski ovi o + tat {3x ~2 s . 't Ire t3xj[3xj ~ ~ (1) (2) where valve are the Cartesian velocity components and fi = (O. - 1/Fr2) is the gravitational force. The equations are nondimensionalized by the maximum inlet velocity U and a length L characteristic of each case examined. The Reynolds and the Froude n~,m- bers are Re = UL/zJ and Fr = U/~fi~, where ~ is the kinematic viscosity and 9 the acceleration of gravity. When the free surface is present and is adverted by the unknown velocity field, the physical domain varies in time. The numerical solution of the gov-
erning equations describing such complex phenomena requires an accurate computational scheme. Among the several possibilities to represent the behaviour of the free surface, such as the Marker and Cell method [3], in this paper it has been chosen to map the phys- ical domain (x`) into a Cartesian computational do- main (6j) at each time step. Interpolations to enforce boundary conditions are thus avoided. The choice of Cartesian velocity components as dependent variables allows a simple form of the equa- tions in the (I coordinate system. The adoption of other velocity components as unknowns require the in- troduction of higher order metric coefficients, as shown in Ref. A. Introducing the time dependent coordinate trans- formation xi = xi((j,t) into equations (1) and (2), it follows 1 8~jVi) _ 1 ski _ o J B6i J B(i ~ + ~ JO (qjoi) = __~j p + ~ ~ ~ ak! ~ Vi + fi t4y J (3) where qi = Levi are the fluxes and qi = qi_)j(~` ac- counts for the time dependence of the mapping. The notation is given in Appendix A. Eq. (3) is a very important constraint for incom- pressible flows, which is fulfilled in a simple way by locating the fluxes qi at the cell sides. The equations for fluxes qs are obtained by a lin- ear combination of (4~. The momentum equations become it, ~ ~` LIZ q ~` + ~]j ~5 (q q Ink) = _o~ij UP + ~ of ~ ~kZ ~ qa~j + pi f j (by where Eq. (5) can be rewritten as aq RA2 (qi) = _ pi _ R Ri = Eli _ Fs + G ~ (D~ +D2) (7) the convective term is the pressure gradient is Hi = .T1i,'~ (q q Ilk) pi _ - the body force term is - ails IMP (8) (9) F = ~ f (10) and the quantity G' = ~qT ,'`zip phi (11) accounts for the time dependence of the grid. The following definitions have been used: ^2`qi; = ~ ~ ~~ Nikko: -tip' no sum over i (12) and Ds = ~ Hi ~ C!kiquip S ~ i (13) D2 = ~ ~ ~~ 't ~~ q p~ (14) where /\2 contains only second derivatives of qi along Ok, while the remaining part of the diffusive term is accounted for in D' and D2 3 Numerical Mode! Eq. (6) are discretized on a staggered grid where the pressure p is defined at the node (i, j), A at (i+1/2, j) and q2 at (i, j + 1/2~. The discretized equations need the metric coefficients at (i,j), (i + 1/2,j), (i,j + 1/2) and (i + 1/2, j + 1/2~. These metric quantities are evaluated by second order centered finite differ- ences. Interpolation of the metric coefficients intro- duces truncation errors A. To avoid such problem, in the present work, the metric quantities are eval- uated on a grid twice finer than the computational grid. The governing equations are solved by a frac- tional step method [6,7,8], where a first step yields the non-solenoidal field ql. q (q ~ _ ~ ^2(qi _ Lint = _ (, Kiln L t3(Ri~n _ (Ri~n- (15) To have a second order time accurate method the Ri have been discretized by the Adams-Bashforth scheme requiring a restriction on the time step (Courant num- ber < 1 ). The direct inversion of Eq. (15) requires a large amount of memory and CPU time. The approx- imate factorization scheme t9,10] reduces the number of operations. It consists on the approximation of the LHS by the product of two three-diagonal matrices. Eq. (15) becomes: (1A1) (1A2) (qi _ (quiz) = i\t(Pi~n _ ~ [3(Ri~n(Ri~n-l] + + 2(A1 + A2~(qi~n Where Ak accounts for the discrete differential oper- 60g
ators of the laplacian along Eke The velocity field (qs~n+t is obtained by introduc- ing a scalar 4} at the second step, ~qi~n+l = qi _ Id §~ (17) By applying the discrete divergence operator at the grid point (i,j), the equation for 4} follows 1 ~ aims 2 {q- the kinematic condition J phi Ads ~/~` Phi (18) A significant amount of the CPU time of the en- tire computation is devoted to the solution of Eq. (18~. Direct methods, like LU decomposition, yield a solu- tion within round-off error, but their application to very fine grids (N = N' x N2) is impossible both for memory requirements and for number of operations CONS. Generallyiterative methods do not efficiently reduce low wave number errors. For such reason a lin- ear Correction Storage multigrid method [11,12] has been applied with a reduction of CPU time. The pressure field p appearing in Eq. (6) is evalu- ated from 4) 2 ( 2Re ~ k~ t~ Act) After the second step the overall accuray remains of O(~\t2~. 4 Boundary Conditions No-slip boundary conditions have been assumed on solid walls. An assigned profile of the Cartesian com- ponent vt together with v2 = 0 has been imposed at the inlet, and the fluxes are derived by Eq. (A.6~. At the outflow a developed Dow condition Levi Gina =0 (20) has been assumed when free surface problems are treated. The assumption of radiative conditions, that is done for the other cases, gives rise to instabilities in presence of free surface. Work is in progress to avoid this numerical instability. The Navier-Stokes equations require one kinematic and two dynamic conditions at the free boundary. The balance of the forces between the inner and the outer fluid gives the dynamic conditions. With zero external pressure and stresses, the normal and shear stresses must vanish. These conditions cannot be easily imple- mented by an implicit scheme. Therefore the simpler conditions of ref. t15] haste been used: Su' = Thai 0 (21) For slightly deformed grids the previous conditions are approximately equivalent to require the vanishing of normal and tangential stresses. The kinematic condition enforces the fluid parti- cles to remain on the free surface. If the free surface configuration is described by the function z2=~(=l,t) (22) D(x2rI) = 0 (23) yields `~77 = v2 _ v: 077 = v2vi C1 (24) At every time step the free surface is moved us- ing the previous equation. To eliminate the possible occurrence of instabilities, the new configuration is smoothed with the same filter used by [13,14,15~. 5 Results To test the complexities associated with the described numerical scheme, several cases have been chosen, each one enlighting a particular complexity. The flow in a steady domain and the flow in a domain where a boundary moves with a prescribed law have been con- sidered as test cases before solving the free surface flow past a semicylinder. 5.1 Flows in steady complex domains The use of a general curvilinear coordinate system is the main characteristic of the present method. There- fore to analyze the influence of the grid on the numer- ical solution the flow past an obstacle inside a channel has been considered. The sketch of the domain and the mesh obtained by an analytical function based on a conformal transformation is given in Fig.1. Fig. 2 shows, by a stream function plot, that the numerical method reproduces the main features revealed in the experiment of Ref. A. To capture the relevant de- tails at least a 64 x 32 grid with refinements in the regions of flow reversal is necessary. The reattach- ment length Xr of the main separation bubble, is the significative parameter of this flow. The numerical values of Xr, together with the experimental data, are given in table 5.1 for two Reynolds numbers {,Re = 144 and Re = 288`J. The numerical simulations have been Re Xr Num. Xr Exp. 144 5.5 6.7 288 9.5 12.1 Table 1: Values of Xr by present work and by Ref. [1] 609
done respectively by a 64 x 32 and a 96 x 32 grid. The coarseness of the grid is the main reason to explain the disagreement between the numerical and the experi- mental results. A second reason, to a lesser extent, is ascribed to the fact that the calculation was inter- rupted before a real stationary solution was reached. As it occurs for the backward facing step, numerical solutions usually present a very steep growth of X, within a short time. Later on Xr grows in time very slowly and the steady state is reached using a large amount of CPU time. However the obtained results can be acceptable as a presentation of the method. A more accurate computation to obtain a physically interesting solution will be presented in the future. 5.2 Flows in complex domains mov- ing with a prescribed law In order to verify the accuracy of the numerical scheme to describe flows driven by an unsteady boundary, the flow inside a channel with a moving indentation has been considered. For such a case experimental [16] and numerical [17] results allow to test the treatment of the grid time dependence. The numerical results of ref. [17] have been obtained by a vorticity-stream function method on a very fine grid (1440 x 48~. The same time dependent law of the indentation consid- ered in Ref. [16,17] has been used. Although this flow is not relevant to ship hydrodynamics, it has been con- sidered as a test case. A sequence of counter rotating eddies close to the upper and lower walls and sepa- rated by a wavy core flow is generate downstream the indentation. The eddies move downstream with a de- fined group velocity. This physical problem has been solved for Re = 760 arid St = 0.025 and the same flow time history reported by [17] has been obtained. The time sequence of istant~neous streamlines is re- ported in Fig. 3 for a sequence of time steps within one period and it shows the same position, shape and temporal development of the eddies of Ref. t17~. The values of the positions of the crests and troughs cor- responding to eddies B. C and D are shown in Fig. 4 in comparison with numerical t17] and experimental [16] values. In this case the better agreement derives from the much finer grid (384 x 32) used in this case with respect to that used in the previous case. These preliminary numerical results assess that the numerical model provides good flow simulation in presence of steady and unsteady irregular domains. 5.3 Free surface flows free surface flows involve the difficulty of irregular and time dependent domain, where boundary moves according to the velocity field. The correct enforce- ment of free-slip conditions gives the deformation of the free surface which strongly affects the inner veloc- ity field. To our knowledge few experimental results concerning free surface flows are available. One ex- periment [2] deals with the flow over a semicylindrical obstacle. We have chosen to verify the validity of the approximated form of the free surface boundary con- dition specified by Eq. (21~. A proper treatment of the free surface conditions is very difficult particularly when implicit schemes are used. In this case the ex- act free-slip conditions require an iterative procedure which increases the CPU time. Preliminary results were obtained at lower Re and higher Fr then those used in Ref. [2] because a coarse grid was employed. However also in this case the free surface is largely deformed. The velocity vector plot, shown in Fig. 5 for Re = 400 and Fr = 0.5, qual- itatively resembles the expected flow field. Also the dynamic pressure reported in Fig. 6 and the develop- ment of the recirculating regions qualitatively agree. However, it is to be noted that the flow field simu- lation is not possible up to the steady state because the free surface conditions are not correctly enforced when the deformation of the interface becomes rele- vant. A better treatment of such boundary conditions is required to achieve a complete and accurate flow simulation. 6 Conclusions The purpose of this paper is restricted only to the description of the numerical method, which has been tested by using coarse grids, certainly larger than those required by the physics of the flow. This is the main reason for the discrepancies between the numerical and experimental results in the two cases firstly con- sidered. More work should be done to define the best grid transformation for the wall mounted obstacle case. On the contrary, the case of the flow generated by a moving indentation shows that the grid time de- pendence is very accurately handled by our scheme. In fact our results computed on a relatively coarse mesh are in very good agreement with the numerical results obtained with a larger number of arid Nina.: in Ref. [17~. At the moment the lack of the method to accu- rately describe free surface flows is mainly related to two intimately connected aspects: the choice of the grid transformation and the approximations on the boundary conditions at the free surface. The present results show that the generation of a grid with co- ordinate lines orthogonal to the boundaries reduces the truncation errors. If such effect is large for solid boundaries, the effect is amplified at the free surface. The generation of coordinate lines nearly orthogonal 610 O ~ _ _ a_ _ _ ~ _ Am
to the free surface is more difficult and requires a high CPU time but avoids the approximations in at the free surface boundary conditions. In the future our work will be focused on these aspects and only afterwards real ship related Bows will be successfully simulated. Appendix A Metric Notation In this section a brief review of the tonsorial notations used throughout our paper is given. Let us consider the coordinate transformation be- tween a Cartesian reference system xi and a curvilin- ear system A, being the base vectors respectively e' and gj. If we introduce the coefficients c, defined as Ad (Aft) the Cartesian components vi of a vector v are related to the contravariant components uj [37 v = vied = ujgj (A.2) by the transformation . . . V$ = CS u] (A.3) If we consider the inverse transformation c~i then pi = ~c~~. vi (A.4) If the vector v introduced above represents a ve- locity field the fluxes qj are defined as qj = Juj = J(c-~'vi =~y~vi (A.5) and v = J 1c~qj = I3~.q' (A.6) where J = detect) is the Jacobian of the transforma- tion, eye = J(c-~' and,l'7 = J-~c'. Using the matrix c we can write the metric tensors as 9ij = gi gj = cici. (A.7) gij = g' . gj = (c~1)~(c~l)' (A.8) 9 = det(gij) = J2 (A.9) The gradient and the Laplacian operators in curvilin- ear coordinate are V¢)i = gi] 0¢ aci ( aci) (A.10) (A.ll) the metric quantities aid which appear in the equa- tions for fluxes are simply given by Hi} = Jgi] (A.12) References [1] Tropez C.D., Gackstatter R. UThe Flow Over Two-Dimensional Surface-Mounted Obstacles at Low Reynolds Numbers J. Fluids Egg. 107 (1985), 109 Miyata H., Matsukawa C., Kajitani H. "Shallow Water Flow with Separation and Breaking Wave" read at The Autumn Meeting of The Society of Naval Architects of Japan Nov. 1985 Welch J.E., Harlow F.H., Shannon J.P., Daly B.J. "The MAC method: A computing tech- nique for solving viscous, incompressible, tran- sient Duid-flow problems involving free surfaces", LA 3425, Los Alamos Scientific Labs. (1966) Levi Civita T., The absolute differential calcu- lus" Dover, 1977 Thompson J.F., Warsi Z.U.A., Mastin C.W. "Numerical Grid Generation" North-Holland, 1985 ~ Chorin A.J., "A Numerical Method for Solving Incompressible Viscous Flow Prblems" J. Comp. Phys. 2 (1967),12 t7] Kim J., Main P. Application of a Fractional- Step Method to uncompressible Navier-Stokes Equations" J. Comput. Phys. 59 (1985), 309 t8] Orlandi P., Kim J., Main P. Numerical Solution of 3-D Flows Periodic in One Direction and with Complex Geometries in 2-D" unpublished t9; Beam R.M., Warming R.F., "An Implicit Finite- Difference Algorithm for Hyperbolic Systems in Conservation-Law Form" J. Comp. Phys. 22 (1976), 87 [10] Briley W.R., McDonald H., "Solution of the Mul- tidimensional Compressible Navier-Equations by a Generalized Implicit Method" J. Comp. Phys. 24 (1977), 372 611
[11] Brandt A. "Multi-level Adaptive Grid Solutions to Boundary Value Problems" Math. Comp 31 (1977), 330 [12] P.G. Esposito, P. Orlandi "A Multigrid Solver for the Pressure Equation for Idcompressible Navier-Stokes Computations in Curvilinear Co- ordinates" unpublished [13) Longuet-Higgins M.S., Cokelet E.D. ~The Defor- mation of Steep Surface Waves. I. A Numerical Method of Computation" Proc. Roy. Soc. A 350 (1976), 1 [14] Haussling H.J., Coleman R.M. "Finite Differ- ence Computations Using Boundary-fitted Coor- dinates for Free-surface Potential Flows Gener- ated by Submerged Bodies" in Proceedings of the 2r~d Ir~tern. Conf. or~ Numerical Ship Hydrody- namics, Berkeley, 1977, 221 . H [15] Miyata H., Sato T., Baba N. "Difference Solution of a Viscous Flow with Free Surface Wave about an Advancing Ship" J. Comput. Phys. 72 (1987), 393 [16] Pedley T.J., Stephanoff K.D. ~Flow along a Ch~-nnel with a Time-Dependent Indentation in One Wall: the Generation of Vorticity Waves" J. Fluid Mech. 160 (1985), 337 [17~ Ralph M.E., Pedley T.J. "Flow in a Channel with a Moving Indentation~ J. Fluid Mech. 190 (1988), 337 I I I I I I I I t_1_~=~ J- l- ~L I I I I I I I I I I I I I I I I I I I l l l l 1 LL I I~ ~ 1 1 1 1 1 T 1 1 1 1 1 1 1 1 T I T 1 1 1 1 I I I I I ~! ~ ~ _ 1 1 1 1 1 1 1 T I I I 1 1 1 1 1 1 1 1 1 1 L I I I ~ Ih ~ I I I rl I T ~ T T I I T T rT 1 ~ I I I ,,,,,, ,~ ~ l · ,~ ,,,, I , I , , , , , I , , I , , , , ~ Fig. 1 A sample 32 x 16 grid for the solution of the Dow past a surface-mounted obstacle in a channel (h/H = 0.5, l/h = 4) xr Fig. 2 Stream function plot for Re = 288; dashed lines correspond to negative values. The insterval between contour lines is 0.02 for negative values and 0.1 for positive ones. 612
t =0.3 t = 0.4 t=O.S t=0.6 A t =0.7 E t =0.8 t =0.9 t = 1.0 Fig. 3 Instantaneous streamline plots within one period for the Dow in channel with a moving indentation 613
1.00- 0.80- 0.60 - 0.40- 0.20 - 0.00- 1.00- O.BO- 0.60- 0.40- 0.20 0.00 1.00- 0.80 - 0.60 0.40- 0.20 - 1 1 1 1 1 1 1 1 1 1 1 0.00 0 1 2 3 4 5 6 7 8 9 10 11 0 1 · Exper. ~ Numer. Am Pres. Work x (a) Eddy B A// I' 1 1 6 7 8 9 10 11 · Exper. Numer. Pres. Work 0 1 x (b) Eddy C 1 1 2 3 4 5 x (c) Eddy D Fig. 4 Time histories of the eddies B. C and D for the flow in a channel with a moving indentation 614 ~7 l 1 i 1 1 2 3 4 5 6 · Exper. ~ Numer. +} Pres. Work 1 1 7 8 9 l 10 11
- = - = / - ~ l ~ JO :?ig. 5 Velocity vector plot at t=0.8, Re = 400, Fr = 0.5 for the flow over a semicylindrical bump Fig. 6 Pressure contour map at t=0.8, Re = 400, Fr = 0.5 for the flow over a semicylindrical bump. The increment between two isolines is O.OS 615
