CFD-FEM 2D MODELLING OF A LOCAL WATER FLOW. SOME NUMERICAL RESULTS Antonio Pasculli Department of Geotechnology, University of Chieti Campus Universitario Madonna delle Piane, 66013 Chieti (Italy) e-mail: a.pasculli@unich.it Abstract: Pasculli A.: CFD-FEM 2D Modelling of a local water flow. Some numerical results (Chieti, Italy) (IT ISSN 0394-3356, 2008). The paper describes a first step towards the elaboration of flexible numerical research tools to study fluvial erosion, sediment tran- sport, Fluid-Structures Interaction (FSI), through the application of the Computational Fluid Dynamics (CFD) and the Finite Element Method (FEM) techniques. A preliminary but necessary item to be explored is the analyses of the robustness of both mathematical and numerical features of the selected equations and the implemented algorithms. The unsteady Navier-Stokes equations, aimed to model a local 2D water flow, are introduced. To include turbulence, the URANS (Unsteady Reynolds Averaged Navier Stokes) Two Equations κ - ε model, with the Wall-Functions proposed by S. FAN et alii (1993), have been employed. At this step, important pheno- mena, such as free-surface influence, bed and banks morphologies changes, coupling between particles and fluid flow, have not con- sidered. Then the adopted numerical approach is discussed. Spatial discretization is carried out by the linear Finite Element Method (FEM). A structured meshing with h like adaptability algorithm was developed. Thus to avoid velocities and pressure instabilities rela- ted, respectively, to convective mode predominance respect to diffusive mode and incompressibility constraint (divergence free flow), the Characteristic-based split (CBS) algorithm and the method of Artificial Compressibility (AC) have been applied. The Dual Time Stepping method has been adopted with Internal External Local Time Stepping. The three steps of the CBS based on AC or the CBS- AC method (NITHIARASU & LIU 2006) are then applied and discussed. As first step, only some steady state parametric numerical experi- ments have been performed, considering a semi-implicit, approach. Poiseuille flow test has been carried out also to clarify and to investigate in depth CBS performance. Both laminar and turbulent regimes have been considered. The numerical results of the selec- ted turbulence modelling are compared with some simple analytical expressions, valid for Poiseuille flow. Useful and important sugge- stions regarding the numerical tuning of some intrinsic parameters, specific to the selected algorithms, are acquired through the discussion of the general performance of the implemented algorithms. Riassunto: Pasculli A.: Modellazione CFD-FEM 2D di un flusso locale acquoso. Discussione dei primi risultati numerici (Chieti, Italia) (IT ISSN 0394-3356, 2008). Il lavoro presentato in questo articolo descrive la prima fase di un'attività di ricerca il cui scopo è di elaborare degli strumenti numerici atti a studiare in maniera flessibile processi quali l'erosione delle sponde fluviali, il trasporto di sedimenti, l'interazione Fluido Struttura (FSI), mediante l'applicazione della Fluidodinamica Computazionale (CFD) e del Metodo agli Elementi Finiti (FEM). Innanzitutto è neces- sario analizzare la robustezza sia degli aspetti matematici sia degli aspetti numerici delle equazioni e degli algoritmi implementati. Sono state considerate le equazioni non stazionarie di Navier-Stokes per la simulazione di un flusso locale d'acqua in 2D. La turbolenza è stata inclusa mediante il modello a due equazioni κ - ε, ottenute tramite le equazioni mediate di Navier Stokes (URANS, Unsteady Reynolds Averaged Navier Stokes) con le “Funzioni di muro” (Wall-Functions), proposte da S. FAN et alii (1993). In questa fase iniziale, non sono state considerate importanti caratteristiche dei fenomeni in considerazione, quali l'influenza della “superficie libera”, il cam- biamento della morfologia al contorno, ossia sponde e letto fluviale, l'influenza reciproca tra i sedimenti e il flusso stesso. Quindi si è discusso l'approccio numerico selezionato. Per la discretizzazione spaziale è stato scelto il metodo agli elementi finiti (FEM) lineari di tipo “strutturato”, con infittimento di elementi triangolari (h like adaptability). Per evitare le instabilità numeriche connesse con le velo- cità e le pressioni, in particolare in regime a prevalenza convettiva piuttosto che diffusiva e per il soddisfacimento del vincolo di incom- pressibilità (con divergenza nulla delle velocità stesse), sono stati scelti, rispettivamente, il Metodo delle Caratteristiche applicato all'e- quazione della conservazione della quantità di moto suddivisa nella sua parte convettiva e diffusiva (Characteristic-based split algo- rithm CBS) e il metodo della compressibilità artificiale applicata, in particolare, all'equazione di continuità (Artificial Compressibility AC). Inoltre è stato applicato il metodo del doppio incremento temporale locale, “interno” ed “esterno” (Dual Time Stepping; Internal External Local Time Stepping). L'algoritmo CBS-AC (NITHIARASU & LIU 2006) si risolve mediante un approccio iterativo in tre passi succes- sivi. In questo lavoro è stato analizzato solo il moto stazionario. Sono stati elaborati degli esperimenti numerici mediante un approccio semi-implicito. Per un confronto con soluzioni analitiche si è scelto il classico moto alla Poiseuille. Lo scopo principale è stato l'ap- profondimento numerico delle modalità con le quali il metodo CBS risolve le equazioni. Sono stati considerati sia il regime laminare, sia il regime turbolento. In particolare, i risultati numerici relativi ai modelli di turbolenza sono stati confrontati con espressioni analitiche semplificate, valide per flussi alla Poiseuille. Mediante la discussione dei risultati, sono state ricavate delle utili e fondamentali indicazio- ni sui valori numerici di alcuni parametri, caratteristici degli algoritmi introdotti, che regolano la convergenza degli stessi. Key words: CBS stabilization, RANS turbulence, numerical modelling. Parole chiavi: stabilizzazione CBS, turbolenza RANS, modellistica numerica Il Quaternario Italian Journal of Quaternary Sciences 21(1B), 2008 - 215-228 INTRODUCTION The experimental studies of important processes and mechanisms, such as fluvial erosion, mass failure, sediment transport, are severely affected by intrinsic difficulties to acquire field data collections. A support for the numerical evaluation of important mechanical parameters, namely shear stresses near fluvial banks and bed rivers, the exchanged energy between the flow and the eroded material, particles lift and depositions, has been searching within the framework of hydraulic models and CFD approaches. Some models have been 216 A. Pasculli proposed using empirical data sets obtained from labo- ratory channels (SIMONS D.B.& SENTÜRK F. 1977, KNIGHT et al., 1984). Recently some progress has been made in using analytical models (KE A N J.W., SM I T H J.D. 2004, 2006). Although these approaches are promising, it is not yet clear whether such methodologies are entirely appropriate (RI N A L D I M. & DA R B Y S.E. 2008). Thus the practise of using CFD modelling techniques as useful support to evaluate field data in river flows, that are dif- ficult or impossible to measure, has now become esta- blished for a range of open-channel flow contests (LANE S.N., et al., 2000 ; DARBY S.E. et al., 2004). Actually, ero- sion and sediment transport phenomena involve multi- ple interactions among fluid flow, particles, moving boundaries, geotechnical characteristics of banks and bed. Thus the main issues are the numerical solution of the equations related to the selected models. Namely, Navier-Stokes and Newton's laws equations have to be solved with their coupling. In order to pursue this objec- tive, the CFD and the Computational Granular Dynamics (CGD) (POSCHEL, T.& SCHWAGER T. 2005) tech- niques could be applied. Within the framework of these methodologies several approaches are available. The phenomena complexity implies, necessarily, a splitting of the modelling in two or more different spatial scales: macro scales of some tenth meters wide or more, and smaller scales of few meters wide or less. Regarding the latter, a more accurate and, for this reason, more sophisticated and more CPU consuming time modelling is necessary. Very accurate local calculation could be employed to evaluate stochastically (ILIOPOULOS I., MITO Y., HA N R A T T Y T., et al., 2003; PA S C U L L I A. & SC I A R R A N . 2006, 2007) the main mechanical parameters to be introduced in cheaper models like the Shallow Water or others, even if less accurate. The most accurate approach is the D i r e c t Numerical Simulation (DNS) of fluid and particles flows (LI C., MOSYAK A., HETSRONI G. 1999, SCHMEECKLE M.W. & NELSON J.M. 2003, ZOHDI T.I. 2007). On the other hand, its application on large or medium scales is not very practical due to huge CPU time and memory, reque- sted for high and either medium Reynolds number. In order to get an unsteady high-frequency representation of the solution by a less expansive approach, L a r g e - Eddy Simulation (LES) has been investigated by many authors (RO D I W., et al.1997; SA G A U T P. 2002;). This technique, which is based on a low-pass filtering of the exact solution of the Navier Stokes equations, makes it possible to obtain a significant reduction in the com- plexity of the simulation by reducing the number of degrees of freedom. But LES is still subject to severe constraints when wall bounded flows are considered, because (at least theoretically) the internal region of the boundary layer needs to be quasi-directly resolved, yielding large computational costs. Because an accura- te unsteady description of the solution is not needed everywhere when dealing with practical engineering problems, the idea of using zonal approaches has emerged (QU É M É R E P. & SA G A U T P. 2002). DNS and/or LES could be used in small localized subdomains where an accurate description of the flow is requested, while computing the rest of the configuration with a low-accuracy method. Commercial CFD computer codes, in spite of the very large and qualified international research activities, several times lack of flexibility and internal intelligibility, so they do not always assure the solutions of those problems occurr ing within the Computational Geosciences framework and which could be solved (few) in enough satisfactory fashion with the now-days mathematical and numerical tools (AB A N T O J. et al. 2005; NITHIARASU P., ZIENKIEWICZ O.C. 2006, pag. 5544). Thus the aim of this work is to describe a first step of the process aimed to develop a research nume- rical tool, suitable for the employment of experimental correlations, theoretical models and so on, with the necessary flexibility and intelligibility. The starting point has been the 2D numerical implementation of the Reynolds Averaged Numerical Simulation (RANS) with its U n s t e a d y version (URANS), belonging to the low- accuracy method set, but less expensive than the afo- rementioned ones. Among many others (WARNER J.C. et a l. 2005) the k - ε two-equation turbulence closure models has been selected (LAUNDER B.E. & JONES W.P. 1972; JO N E S W.P. & LA U N D E R B.E. 1972). The spatial discretization have been realized by FEM approach. Linear triangular elements have been selected. The adoption of the Finite Element Standard Galerkin Methods implies instabilities in particular for convective mode predominance respect to diffusive mode . Furthermore incompressibility causes pressure instabili- t i e s. Several methods have been proposed, including Streamline upwind Petrov-Galerkin (SUPG), G a l e r k i n Least Squares (GLS), Finite Calculus (FIC) and more recently Subgrid Scale (SGS) approach (CH U N G T . J . 2006). In the present work the C h a r a c t e r i s t i c - b a s e d split algorithm (CBS) and the Method of Artificial com- p r e s s i b i l i t y (AC) have been applied. For the transient dual time stepping method has been adopted with Internal External Local Time Stepping (NITHIARASU & LIU 2006; NITHIARASU P. & ZIENKIEWICZ O.C. 2006). One pro- blems of the interaction between the fluid flow and both banks and bed surface is computing turbulent flows that are influenced by an adjacent soil (wall) boundary. Thus the 2D Wall-Functions, proposed by S. FAN et alii (1993), have been implemented. Since the main purpose of the present work is to verify the robustness of the implementation of the selected fluid flow and turbulence models and their mean numerical features, many important issues regar- ding fundamental characteristics of the phenomena under study have not been considered. Among others Free - Surface and Bed- Fluid Interaction effects (SOULI M. & ZOLESIO J.P. 2001; BAGHLANI A. & TALEBBEYDOKHTI N. 2007), Meshing - Updating (XU X. & HARADA K. 2003). Then some parametric studies regarding h adapti- ve structured meshing have been carried out. Poiseuille test has been utilized to explore how CBS methodology can affect the stabilization of the fluid velocities. Some important parameters had to be correctly tuned to reach convergence of the iterative procedure. The outcomes of the assumed turbulence model- ling obtained by turbulent Poiseuille test are compared with the Prandtl's Mixing Length Theory and Deissler's sublayer and buffer zones models. The integral equation related to the latter model has been solved by an itera- tive approach. Future developing will concern with further RANS analyses and the link among different-accuracy approa- ches applied at different scales: Shallow Water, URANS, LES, DNS, CGD. 1. A GENERAL 2D FLUID DYNAMIC MODELLING In this paragraph a general 2D fluid dynamic modelling, developed in both its mathematical and numerical features will be described and discussed. The algorithms, translated in a research computer code performed by the author by Visual Fortran 2000 Compilator through the implementation of about twenty thousands instructions, have been tested and some numerical results have been discussed in this paper. The physics conservation laws are the mathemati- cal framework of fluid dynamics. Thus the usual Navier Stokes equations have been considered (in S.I. units). 1.1 Mass conservation In Einstein summation convention: ∂ρ ∂ (ρ ui) –––– + ––––––– = 0 i = 1,2 (1) ∂t ∂xi where ρ is the water-sediments volumetric mass den- sity, ui the i th velocity component. 1.2. Momentum conservation equation For the momentum equation it is necessary to introduce the dissipative term due to water-sediment viscosity. In order to build constitutive laws related to viscous fluids the Strain Temporal Rates Tensor has to be introduced, which, for small deformation assumes the following expression: 1 ∂ui ∂uj ε. ij = –– ( –––– + ––––) (2) 2 ∂xi ∂xj It is expressed also by means of the following vector: ε. = [ε.11, ε . 22, 2ε.12] T (3) The relation between the Strain Tensor and the Stress Tensor, for Linear Newtonian and Isentropic Fluids is provided by two terms: the Deviatoric Stress- Strain Term and the Isotropic Term. The first one is: 1 1 τij = σij – –– δij σkk = 2µ (ε. ij – –– δij ε . kk) (4) 2 2 Where µ is the dynamic viscosity. The complete form is: ∂ux ∂uy ∂ux ∂uy(––– – –––) (––– + –––)∂x ∂y ∂y ∂x τ = µ ( ) (5) ∂uy ∂ux ∂uy ∂ux(––– + –––) (––– – –––)∂x ∂y ∂y ∂x The second term, Isotropic Term, is defined as pressure: 1 σ11 + σ22 p = – –– σkk = – –––––––– (6) 2 2 Volume viscosity has been neglected. In the Turbulent Regime flow, a further stress term, due to Turbulent Kinetic Energy k, should be introduced 217CFD-FEM 2D Modelling of a local ... ∂ui τt i j = µt ––– – ρk where ρ is the volumetric water sedi- ∂xj ments mixture density, µt is the turbulent dynamic viscosity which will be described in following para- graphs, while: k = half of the turbulent square velocity = –1 2 (u’2 x + u’2 y ) where u’2 x , u’2 y are velocities fluctuation. µ is the laminar dynamic viscosity. Thus: ∂ρui ∂[(ρuj) . ui ] ∂(τij + τt ij ) ∂ρ –––– + ––––––––– – –––––––– + ––– – ρgi = 0 i=1,2; j=1,2 (7) ∂t ∂xj ∂xj ∂xi 1.3. Energy conservation equation To complete the review of all the base equations necessary to elaborate a numerical model of a 2D local water flow including phenomena like fluid-structure interaction, free-surface, bed and river-bank erosion with their topography change, sediment deposition and flow mass accretion, it is mandatory to show, at least, the equation related to the advection of a scalar quan- tity, affected by a turbulent regime as well. Thus the form of the mechanical energy is briefly reported and described as example, notwithstanding its contribution is not discussed in this paper. In a 2D mechanical flow model, excluding thermo-chemical contribution and considering gravitational potential energy as a volume- tric source, the Specific Energy for Mass Unity is related only to Kinetic component, defined as: 1 E = – (u2 x + u2 y ) (8) 2 Applying the Energy Conservation Principle t h e following equation is deduced (in a compact form): ∂ρE ∂[(ρui) E ] ∂(ρui) ∂(τij + τt ij ) uj –––– + –––––––– + ––––– – –––––––––– – ρgi ui – Q = 0 (9) ∂t ∂xi ∂xi ∂xi ∂ui Where, in particular τt ij = µt ––– – ρk as previously intro- ∂xj duced, while ρgi ui is the variation during the transitory of the potential energy due to gravity. The Q [W] [m]- 3 term is the Volumetric Energy Source which could include, for example, buoyancy forces effects due to spatial variability of water-sediments mixture. 2. TURBULENT MODELLING Among others, the RANS (Reynolds Averaged Navier Stokes) Two Equations k - ε modelling with the 2D W a l l - F u n c t i o n s proposed by S.FAN et alii (1993), has been implemented. It consists, fundamentally in two equations which describe the generation, the diffu- sion, the convection and the dissipation of the energy due to turbulence. 2.1. Production and transport equation of the turbu- lent kinetic energy The turbulent kinetic energy is considered as a scalar quantity. Thus by means of the General Scalar Conservation Equation, the following relation yields: ∂ρk ∂ [( ρux) . k] ∂ [( ρuy) . k] ∂ µt ∂k ––– + –––––––––– + –––––––––– – –– [(µ + –––) ––– ] –∂t ∂x ∂y ∂x σk ∂y (10) ∂ µt ∂k – –– [(µ + –––) ––– ] + (P + G – ρε) = 0 ∂y σk ∂y where [k] = [m]2[s]-2 = semi square of total Turbulent k2 Kinetic Energy, µt = ρCµ fµ ––– turbulent viscosity; ε ∂ux 2 ∂uy 2 ∂ux ∂uy 2 P = µt { 2 . [ ( –––) + ( –––) ] + ( ––– + –––) } is the source ∂x ∂y ∂x ∂y of the turbulent kinetic energy due to Shearing Stress; G is the turbulence source due to Buoyancy Forces. 2.2. Production and transport equation of the dissi- pation rate of the turbulent kinetic energy By means of the same approaches adopted in the previous paragraph: ∂ρε ∂ [( ρux) . ε] ∂ [( ρuy) . ε] ∂ µt ∂ε ––– + ––––––––– + ––––––––– – –– [(µ + –––) ––– ] –∂t ∂x ∂y ∂x σk ∂x (11) ∂ µt ∂ε ε – –– [(µ + –––) ––– ] + –– (Cε1fε1 P + Cε1G –Cε2fε2ρε) = 0 ∂y σk ∂y k Where: [ε]=[m]2[s]-3 is the volumetric density rate dissi- pation of the turbulent kinetic energy. 2.3. Experimental parameters For these kind of turbulence models it is neces- sary to calibrate the mathematical equations with expe- rimental parameters. In this paper the following values have been selected: Cµ = 0.09; Cε1 = 1.44; Cε2 = 1.44; σt = 0.9; σk = 1; σε = 1.3 2.4. Wall Functions The following are the functions which introduce the reasonable lowering of the turbulence near rigid boundaries (Wall Functions) (S. FAN, et alii, 1993): k2 √ −−− Rk √ −−− Rk Rk Rk Rt = ––; fw = 1-exp {– ––––– +(––––– – ––––)[1-exp(– –––)] 3 };νε 2.3 2.3 8.89 20 fw fw Rk fµ = 0.4 ––––– (1 - 0.4 ––––– ) [1-exp ( – ––––––) ] 3 ; √ −−−−− Rt √ −−−−− Rt 42.63 0.4 Rt fε 2 = {1– ––– exp [ – ( –––) 2 ] 3 } f 2 w; 1.8 6 fε 1 = 1 dw Where: Rk = ––– √ − k, dw is the minimum distance of the ν considered point from the closest rigid boundary; µtν = –– is the kinematic viscosity [ν]=[m]2[s]-1 ρ 2.5. Boundaries conditions related with turbulence modelling Along rigid boundaries, the turbulent kinetic energy is assumed to be zero, while on boundary elements ∂k 0,5 2 εw = 2ν ( –––– ) where n is the normal direction. At the ∂n eventual inlet flow, it is assumed: k = 0.05 x √ −−−−−−−−− u2 x + −− u2 y . 3. NUMERICAL MODELLING The previous equations are characterized by non self adjoint differential operators, due to convective terms. As it is well known, in this case, the adoption of the Finite Element Standard Galerkin Methods i m p l i e s instabilities in particular for convective mode predomi- n a n c e respect to diffusive mode. Furthermore i n c o m- p r e s s i b i l i t y causes p r e s s u r e i n s t a b i l i t i e s. Thus to avoid velocities and pressure instabilities the C h a r a c t e r i s t i c - based split algorithm (CBS) and the Method of Artificial c o m p r e s s i b i l i t y (AC) have been, respectively, applied. For the transient dual time stepping method has been adopted with Internal External Local Time Stepping. The CBS based on Artificial compressibility AC or the CBS-AC (NITHIARASU & LIU 2006) method, despite other methods proposed in literature, consists in only three steps. The most important advantage of this procedure is the possibility to implement arbitrary order for the velocity and the pressure interpolation without an artifi- cial damping. In this way the satisfaction of the heavy LBB (BABUSKA 1973; BREZZI 1974) condition is avoided. In this paper the Standard Galerkin FEM (weighting functions are the shape functions) and Eulerian triangu- lar meshing with linear interpolation for both velocity and pressure have been selected (CH U N G T.J., 2006, pag. 273-283). Important consequences are: a m o r e simplicity and affidability of the implementation pro- cess; very fast algorithm for non high Mach n u m b e r flow, with a very good precision. The three steps are summarized as (NI T H I A R A S U P., CO D I N A R., ZI E N K I E W I C Z O.C. 2006): Step 1 ∂ (uj Ui) ∂(τij + τt ij) Δt ∂ ∂ (uj Ui) ΔU*i =Δt [– –––––– + ––––––– + (ρgj) + –– uk ––– (–––––– – ρgi)] n (12) ∂xj ∂xj 2 ∂xk ∂xj where Ui = ρui is the flux vector, while all the other sym- bols meanings are supposed to be evident. The application of the Standard Galerkin FEM, requests the employment of unknowns related to each node. The global solution is carried out by means of an assemblage of each single solution related to connec- ted spatial elements. Thus each unknowns, labelled at a global level as Ui where i=1,N with N = Total nodes number, must be labelled again at a local level as Umi, where m is a local index associated to each node belonging to the selected element. For the linear trian- gle element, m = 1 , 2 , 3 corresponding to triangle vertex (many other options are available). In particular regar- ding eq. (12), the following 3 equations, whose unk- nowns are the flux vector Umi with i=1,2 (for each spatial 218 A. Pasculli component x and y) and with m=1,2,3, should be sol- ved, following the Standard Galerkin FEM approach: ∫Ω Nm (∑3 k=1 Nk ΔUki) dΩ = ∂u–j (∑3 k=1Nk Uki) ∂Nm –Δt [ ∫Ω Nm ––––––––––––– dΩ+ ∫Ω –––– (∑3 k=1Nk (τij+τt ij)k)dΩ – ∂xj ∂xj Δt2 ∂ (u–k Nm) ∂u–j (∑3 k=1NkUki) – ∫Ω Nm∑ 3 k=1Nk (ρgj)kdΩ ] n + –– [ ∫Ω ––––––– (– –––––––––––––+ 2 ∂xk ∂xj + ∑3 k=1Nk (ρgj)k ) dΩ ] n +Δt [ ∫Γ N* (τij+τt ij) nj dΓ ] n where Nm, with m=1 , 2 , 3, are the Shape Functions which, in this case, are a linear x, y, expressions: Nm = am + bm x = cm y where am, bm, cm are simple functions of the global coordinate of the related mt h node (CH U N G T.J. 2006, pag. 273-283; ZIENKIEWICZ O.C., TAYLOR R.L. & NITHIARA- SU P. 2006, pp. 88-91) Step 2 1 ∂Un i ∂ΔU*i ∂ 2 pn ∂2Δp Δρ= (––) n Δp=– Δt [––– +θ1 ––––– –Δtθ 1(–––––+θ2–––––)] (13) c2 ∂xi ∂xi ∂xi∂xi ∂xi∂xi Step 3 ∂pn+θ2 Δt 2 ∂ ∂pn U**i =– Δt –––––– + ––– un k ––– ( –––– ) (14) ∂xi 2 ∂xk ∂xi for which the θ Wilson procedure (also called β Method, CH U N G T.J. 2006, pag. 72) has been adopted as well. The flux value calculated at time step 1 is corrected: Un+1 i = U n i + ΔU*i + ΔU**i , so at the end of the (n+1)th time step: un+1 i = Un+1 i /ρ. By these three steps, fluid-dynamics field (velocity and pressure) is evaluated. The momentum equations are non linear due to convective terms. Thus to linearize the system, the velo- city product of the convective terms are assumed to be formed by the value of the velocity evaluated at the pre- vious step, calculated as an average within the local ele- ment, corresponding to the value assumed by the velo- city at the Gauss point of the element (barycentre), ux1 + ux2 + ux3 uy1 + uy2 + uy3 and the value: u _ x = ––––––––––––– u _ y = –––––––––––––, 3 3 where ux1, uy1 are the element nodal velocity values. For example for the x stability term: b1b1 b1b2 b1b3Δt 2 Δt 2KuUx = – ––– A[u_xu _ x( b2b1 b2b2 b2b3 ) + 2 b3b1 b3b2 b3b3 c1c1 c1c2 c1c3 Ux1 + u _ xu _ y ( c2c1 c2c2 c2c3 ) ] ( Ux2 ) c3c1 c3c2 c3c3 Ux3 Where Δt is the time step, A is the triangular element area, while bi and ci are the parameters characterizing the Shape Functions previously introduced. The com- plete numerical solution of the selected modelling is summarized in the following list: 1. momentum equation solution without pressure terms (splitting procedure); 2. pressure calculation by Poisson equation; 3. correction of the velocities evaluated at Step 1; 4. energy and scalar quantity (like sediment concentra- tion) calculation; 5. calculation of the parameters to be introduced in the k - ε equations; 6. calculation of the production and transport (convecti- ve and diffusive) of the turbulent kinetic energy; 7. calculation of the turbulent kinetic energy dissipation. 4. TIME-STEPS It is well known that in convective predominance flux, in order to avoid instabilities, the time step should assume optimal values (Petrov Galerkin), function of u . h Peclet number Pe = ––––– where u is the velocity, h the 2k mesh characteristic size and k is the diffusion parame- ter (thermal or mass diffusion, kinematic viscosity and so on). Peclet number depends on local conditions. Thus a nodal time step has been introduced by the eva- luation of the minimum time step for each node. 4.1. Artificial compressibility In Step 2 (13), previously described, the Artificial c o m p r e s s i b i l i t y (AC), by means of the adoption of an enough low value of the sound velocity, has been intro- duced. (NITHIARASU P.& LIU C.B. 2006). This aproach is correct in the Steady State situation if any. So in the 1 1 Step 2: ( –– ) n Δp ≈ ( –– ) n Δp where β is an artificial para- c2 β 2 meter with a velocity dimension. It could be calculated by (ZIENKIEWICZ et al., 2006): β = max (ε, uconv , udiff) (15) where ε is a small constant in order to avoid zero values ν of velocities, while uconv = √  u2 x + u2 y  and udiff = – with ν = h kinematical viscosity, h = characteristic dimension of the triangular element. The final local time step is: h Δt = ––––––––––––– (16) √  u2 x + u2 y  + β It is possibile to apply the same AC procedure for the transient as well through the Dual Stepping approach. 4.2. Dual Time Stepping In order to recover the true transient for the artifi- cial compressibility AC, any of each of the three steps (12), (13) and (14) should be modified. Within each real time step Δt an internal iterative loop is carried out to go from time nth to time (n+1)th. Then, let's select the third step (14) in its matrix form and consider the following modified expression (for symbols meanings see NITHIA- RASU P.& ZIENKIEWICZ O.C. 2006): ΔUm Δtint Um – Un ––– = – M-1 u [ GT (pn + θ 2Δpm) + –––– Ppn ] – ––––––– (17) Δτ 2 Δt 219CFD-FEM 2D Modelling of a local ... 220 Fig. 1a - Derivate of the x velocity vs y direction, 20 nodes. Derivata rispetto ad y della velocità lungo x, 20 nodi. Fig. 1b - Derivate of the x velocity vs y direction, 6 nodes. Derivata rispetto ad y della velocità lungo x, 6 nodi. 3Um +1 – 4Um + 4Um-1 Where ΔUm = ––––––––––––––––– which allows the tran- 2 sient to be second order in time (CHUNG T.J. 2006, pp. 46-49), Δpm = pm – pn. When iterating within each inter- nal loop, Δtext and Δtin assume the meaning of pseudo time steps and the meaning of stabilization parameters increasing the acceleration towards the steady state convergence. Therefore for each iteration, a local values h of Δtext = –––––––––––– and Δti n = γΔte x t, with γ>1, obtained √  u2 x + u2 y  + β by attempts aimed to reach convergence, are numeri- cally provided. The iterations stop when the function defined as: √  ΣNodes i = 1 (|U|n+1  – |U|n)2  epsilon = –––––––––––––––––– (18) √  ΣNodes i = 1 (|U|n+1  )2  assumes a value lower than a prescribed tolerance, whose range is usually 10-4 ÷ 10-6. 5. MATRIX INVERSION AND DERIVATIVE In Step 1 (12) and Step 2 (13) of the split procedu- re, the solution only require the inversion of mass matrix related to velocities and pressure. Such steps are called explicit and generally are accomplished using approximation by a diagonal (lumped) form (ZIENKIEWICZ O.C. et alii 2006, pag.91-92). Such lumping for steady state problems make Steps 1 to Step 3 very fast and the related errors are of no consequences as terms involving time variation disappear at a converged (steady state) solution. However, for transient problems, quite serious errors can occur and in such cases an additional iteration is used to obtain a consistent solu- tion (ZIENKIEWICZ O.C. et alii 2006, pag.61-63). Different important parameters, for example the vertical gradient of the water-sediment mixture, depend on the first and on the second derivative. A FEM first derivative can be employed considering the derivative of the FEM approxi- mation of the velocity in each element, whose expression is: u(e) x = L1 . ux1 + L2 .ux2 + L3 . ux3, where L1 = ai + bi x + ci y x3 – x2 (in particular, for example, c1 = –––––– , A= element area, |2A| x2 and x3 abscissas of the selected triangular element, as previously discussed see: CH U N G T.J. 2006, pag. 273-283), uxi is the velocity values at the ith node of the element “e”. Thus: ∂u(e) x (x3 – x2) (x1 – x3) (x2 – x1) ––– = ––––––– ux1 + ––––––– ux2+ ––––––– ux3 (19) ∂y |2A(e)| |2A(e)| |2A(e)| The average is then obtained weighting the expression (19) by the area A( e ) of all the element to which the selected node belongs. It was assumed a parabolic ux (y) velocity profile, ux (m/s) = 1. + y . (1 – y), and a structured uniform meshing of a rectangular strip, similar to that displayed in Fig. 2a. Then in Fig.1a and Fig. 1b, comparison among different calculated value of the first derivative, has been reported. The comparison has been carried out considering the FEM expression, the velocity ratio along a selected vertical line (Difference Element Method) and the exact analytical expression (1–2y), for 20 and 6 nodes uniform grid. The figures show clearly that this kind of FEM approach is poor at the border respect to the DEM model. This is due to the particular selected averaging technique. Thus a particular care should be employed to adopt a different more suitable approach to provide a correct first and second derivatives especially, in par- ticular, close to solid boundary zones within which shear stress has to be evaluated. In general such pro- blems arose where discontinuities can occur. From the above discussion it follows the strong necessity of checking the robustness of the implemen- ted algorithms, before any attempts of doing numerical evaluation of actual phenomena. Thus this is the aim of the following paragraphs. 6. NUMERICAL TESTS Among available fluid dynamical analytical tests designed to check the correctness of numerical algo- rithm, the Poiseuille Problem has been selected. It is supposed that through a 0.5 m wide and 1. m. length A. Pasculli 221 horizontal longitudinal strip a water flow occurs. The Reynolds number is assumed equal to 103. The inlet velocity profile, within the whole transient, has been assumed parabolic along the longitudinal coordinate x, while zero along y. Thus, initially, the velocity is zero everywhere, except at the inlet of course. No slip condi- tion (zero velocities) was assumed on the two longitudi- nal boundaries. At the inlet the pressure is assumed to be p = 9.8.104 + 8.µL Pascal (according to the theory), where L is the length of the strip (equal to 1 m in this case) while, initially, in any other points the pressure is equal to 98000 Pa. The theoretical distribution of the pressure is uniform along direction transversal to flow vector (y direction) and linearly decreasing along the longitudinal coordinate. As the flowing surface is hori- zontal, the hydrostatic load is excluded. In all the calcu- lation an uniform or clustered structured meshing, obtained by 18 nodes along the y direction and 36 nodes along the longitudinal (x) direction, have been selected. Thus the velocity components along x, y and the pressure have been evaluated by the numerical algo- rithm everywhere, outlet section included, except at the inlet and on the two longitudinal boundaries. Furthermore both Laminar and turbulent Poiseuille cases have been studied and compared to each other. In all the “laminar” elaborations carried out in this paper, the following parameters have been selected: θ1 = 1., θ2 = 1. safety coefficient = 0.9; pressure stability c o e f f i c i e n t (rateo between internal and external time step) = 107; β reference velocity (AC model) = 1 m/s. The numerical problem, actually, is 2D. In fact the numerical result value of uy should be vanishing if the selected algorithms and its implementation are correct! As a first testing step, only steady state numerical elaborations have been carried out and discussed in this paper. 6.1. Laminar flow steady state In Fig. 2a and Fig. 2b, the total vector velocity profiles u = uxi + uyj, numerically evaluated for a struc- tured uniform meshing, after different iterative number of steps, have been reported. The y velocity component Fig. 2a - Uniform Structured meshing Discretizzazione Strutturata uniforme Fig. 2b - Laminar Poiseuille test; total velocity profiles at different iterative steps, uniform meshing Test di Poiseuille laminare; profili della velocità totale a differenti passi iterativi con distribuzione uniforme di elementi CFD-FEM 2D Modelling of a local ... 222 uyj, correctly, vanishes. In Fig. 3 the results of an elaboration, carried out without the CBS stabili- zation term, is shown. I t worths to note that without CBS stabili- zation term, the flow plots show an evident instabili- ties and the calculation does not converge. Then, the influence on the numerical results of the non uniform grid meshing has been explored. To this purpose the fol- lowing coordinates transforming algorithm has been introduced: β + 1 x-α –––– 1-α (2α + β ) . (–––––) + 2α – β β – 1 x’ = L . –––––––––––––––––––––––––––––– β + 1 x-α –––– 1-α (2α + 1) [ (–––––) + 1]β – 1 (20) β + 1 y-α –––– 1-α (2α + β ) . (–––––) + 2α – β β – 1 y’ = H . ––––––––––––––––––––––––––– β + 1 yα – ––– 1-α (2α + 1) [ (–––––) + 1]β – 1 where (x, y) e (x’, y’) , respectively, are nodes coordina- tes before and after the transformation. L is the length, while H is the wide or the height of the rectangle, 0 < α e 1 < β < ∞ are two parameters to be adjusted in input and from which the elements clustering density depends. Velocities numerical results, carried out on a cluste- red element grid, Fig. 4a, have been displayed in Fig. 4b. Fig. 3 - Laminar Poiseuille test; velocity profiles without CBS stabilization term Test di Poiseuille laminare; profili di velocità ottenuti senza il termine di stabilizzazione CBS Fig. 4b - Laminar Poiseuille test; total velocity profiles at different iterative steps, lateral boundaries meshing clustering Test di Poiseuille laminare; profili della velocità totale a differenti passi iterativi con addensamento laterale di elementi Fig .4a - Lateral boundaries meshing clustering Addensamento di elementi sui bordi laterali A. Pasculli 223 The plot of Fig. 5 have been obtained after 500 iterative steps, with just only lateral boundaries meshing clustering (Fig. 4a). The velocity profiles show a general good behaviour, but the numerical analysis reported in Fig. 5a and 5b, reveal a substantial differen- ce (about 30%) among numerical and analytical values of the velocities just along the boundaries where the results, for erosion phenomena elaboration, should be more accurate. Any way, it was reached an almost per- fect symmetry (within the 10th decimal digit) respect to the longitudinal axes, due to the symmetry of the grid, not obtainable with a different spatial covering strategy also if through the same type of element. The pressure profile is very satisfactory and uniform (0.45% maximum relative error) showing, any way, the same numerical qualitative trend of the velocities profiles. Figs. 6 set reports the results and the comparison obtained after 2000 iterative steps. The overall results for both velocities (0.9% as maximum relative error) and pressure (0.08% !) are very satisfactory. Furthermore other non reported plots show that the calculated numerical y component velocities, provided by the solution of the related 2D numerical algorithms (theore- tically zero), do not exceed 10-5 m/s. 6.2. Turbulent Poiseuille problem For the turbulent Poiseuille problem, the following parameters have been selected: θ1 = 1., θ2 = 1. s a f e t y c o e f f i c i e n t = 0.2; pressure stability coefficient ( r a t e o between internal and external time step) =2.; β referen- ce velocity (AC model), suggested for high Reynolds number, = 1 m/s. In particular, to obtain the convergence, it was necessary to drastically lower the safety coefficient: from 0.9 down to 0.2. Fig. 7a shows the averaged velo- Fig. 5a - Relative Error distribution of longitudinal velocity Errore relativo della velocità longitudinale Fig. 5b - Relative Error distribution of the pressure Errore relativo della pressione Fig. 6a - Laminar Poiseuille test; total velocity profiles after 2000 iterations, lateral boundaries meshing clustering Test di Poiseuille laminare; profili della velocità totale dopo 2000 iterazioni con addensamento laterale di elementi Fig. 6b - Relative Error distribution of longitudinal velocity Errore relativo della velocità longitudinale Fig. 6c - Relative Error distribution of the pressure Errore relativo della pressione CFD-FEM 2D Modelling of a local ... 224 cities profi le, obtained with an upper and lower meshing clustering after 3000 iterative steps. From Fig. 7a it is possible to outline the capability of these kind of algorithm to provide and to forecast the occurrence of a thin boundary layer, within which a strong velocity variation and for this reason a strong shear stresses occur. Fig. 7b and Fig c show, respectively, the cross distribution of the turbulent energy production and the cross distribution of the temporal rate of the dissipation of the turbulent energy. Their profiles are rather distor- ted especially close to the solid boundaries. Thus a further calculation with the total clustering of all the boundaries have been carried out: Fig. 8a. The following are the selected values of the clustering con- stants: α = 0.5 and β = 1.05. As Figs.8 display, turbulent kinetic energy and its dissipation rateo shapes are closer to a parabolic profi- le than the previous case. Also in this case a boundary layer occurs. 6.3. Numerical-Analytical Turbulent Poiseuille pro- blem comparison In order to make some comparisons with simpli- fied turbulence models, the Prandtl's Mixing Length Fig. 7a - Turbulent Poiseuille test; total velocity profiles after 3000 iterations, lateral boundaries meshing clustering Test di Poiseuille laminare; profili della velocità totale dopo 3000 iterazioni con addensamento laterale di elementi Fig. 7b - Distribution of the Turbulent Energy along y Distribuzione dell'Energia Cinetica turbolenta y Fig. 7c - Distribution of the Turbulent Kinetic Energy dissipation Distribuzione della dissipazione dell'Energia turbolenta lungo y Fig. 8a - Clustering of all boundaries Addensamento di elementi sui tutti i bordi Fig. 8b - Cross Turbulent Energy Energia Turbolenta lungo le sezioni Fig. 8c - Cross Turbulent Energy dissipation Dissipazione dell' Energia Turbolenta lungo le sezioni A. Pasculli 225 Theory and Deissler's sublayer and buffer zones models (BRODKEY R.S. 1967, pp.. 240-252; BIRD R.B. et al.1960, pp. 149-161) have been selected. The comparisons are necessarily semiquantitative since different assump- tions underlay each models. Dimensionless parameters are introduced: u u+ = –– is the dimensionless velocity, where ux is, u* in this case, the longitudinal ux velocity; u* = √  τ0 – is the reference velocity, while τ0 = µ – – – d ux is ρ dy the shear stress at the wall; u* ρ y+ = y –––– is the dimensionless distance from the wall. µ In Prandtl-Deissler approach the turbulence region is divided into three main areas: - viscous or laminar sub-layer 0 ≤ y+ < 5; - buffer layer 5 ≤ y+ < 26; - turbulent core 26 ≤ y+ For each area, separate equations yield. Within the turbulent core Prandtl introduced the logarithmic velocity distribution which can be rearranged in the fol- lowing dimensionless form: 1 u+ = –– ln y* + B (21) k By experimental considerations, Deissler suggested to assume k = 0.36. Further he found that at the border y+ = 26 between the buffer and the turbulent area, the dimensionless velocity assumes the value: u+ = 12.85. It follows that in equation 21 B = 3.8. Many other relations should be applicable in viscous and buffer layers, but this approach is undesirable. Thus efforts have been made to reduce the number of equations by combining the sublayer and buffer zones. Deissler, assuming an exponential form of the eddy viscosity, suggested the following integral equation: dy+ u+ = ∫ v+ 0 –––––––––––––––––––––––––––– (22) 1 + n2 u+ y+ [ 1 – exp ( – n2 u+ y+)] Equation (22) has been solved by a numerical ite- ration methods, whose algorithm has not been reported in this paper. Each integration has been realized by three Gauss points approach. The results are reported in Figs. 9. The analytical values are compared with the numerical results carried out with the nodalization reported in Fig. 8a Figs. 10 are related to analytical and numerical distribu- tion velocities calculated at the section 0.86482 m far from the inlet. Further, also the inlet velocity input has been considered in order to visualize how the parabolic input velocity distribu- tion is deformed by the selected turbulence model. As the Figs. 10 show, the match between numeri- cal and analytical results are satisfactory within the viscous layer, while in the other zones the comparison is not well enough. Other calculations, not reported in this paper, suggest that the problem is due to the ina- dequate number of elements to cover the geometrical section. CONCLUSION In this paper the effort and the care necessary to develop a CFD-FEM research computer code, aimed to perform numerical evaluations related to banks and bed river erosion, sediments transport and structures-fluid interactions have been just sketched. As a consequen- ce, it should be stressed the importance of all the essentially four related steps: the selection of the most suitable physical models by the undertaken study pur- pose point of view, the accuracy and consistence of the mathematical “translation” of the already selected phy- sical approach, the accuracy and consistence of the numerical tools employed to solve mathematical models and some experience and skil l on p r o g r a m m i n g. On the other hand several issues have been not considered, like the impact of the morphology and its change, shock front propagation, free surface flows, inclusion of the Computational Granular Dynamics and many other ones. The discussion of the results points out the sensi- tivity of the employed algorithms not only to the mesh size and to their distribution, but also to some intrinsic “experimental dials” (safe coefficients, explicit vs impli- cit ratio) and to a correct number of iterations, whose knowledge reveals to be fundamental. Any way the results discussed in the paper, regar- Fig.9a - Universal turbulent velocity Distribuzione universale delle velocità turbolente Fig.9b - Universal turbulent velocity in the viscous layer Distribuzione universale delle velocità turbolente nella zona viscosa CFD-FEM 2D Modelling of a local ... ding the comparison among numerical and available analytical values, carried out for P o i s e u i l l e p r o b l e m , related to both laminar and turbulent modelling, show good performance in regions close to the rigid boun- dary in which eventual erosion can occur. Suggestions come up in order to improve the turbulence numerical results by optimizing both distribution and number of the elements. Next efforts will be concerned with more checks on turbulence models, the exploration of other approaches like SPH (Smoothed Partic le H y d r o d y n a m i c s), Meshless method, PFEM (P a r t i c l e Finite Element Methods) which seem to be very promi- sing. ACKNOWLEDGMENTS The author would like to thank two anonymous referees for their effort and useful suggestions. BIBLIOGRAPHY ABANTO J., PELLETIER D., GADON A., TREPANIER Y., REGGIO M. (2005) -Verification of some commercial CFD codes on atypical CFD problems - AIAA Paper, Reno. BA B U S K A I. (1973) - The finite element method with Lagrangian multipliers - Numer. Math. 2 0, pp. 179-192. BA G H L A N I A. & TA L E B B E Y D O K H T I N. (2007) - A mapping technique for numerical computations of bed evo- lutions - Applied Mathematical Modelling 31, pp. 499-512. BI R D R.B., ST E W A R T W.E., LI G H T F O O T E.N. (1960) - Transport Phenomena - John Wiley & Sons, Inc. New York. BR E Z Z I F. (1974) - On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers - Rev. Francaise d'Automatique Inform. Rech. Oper., Ser. Rouge Anal. Numer. 8(R-2), pp. 129-151. BRODKEY R.S. (1967) - The phenomena of Fluid Motions- Addison-Wesley Series in Chemical Engineering. CH U N G T.J. (2006) - Computational Fluid Dynamics - Cambridge University Press. DA R B Y S.E., SP Y R O P O U L O S M., BR E S S L O F F N. A N D RI N A L D I M. (2004) - Fluvial bank erosion in meanders: A CFD modelling approach - In Garcia de Jalon Lastra, D. and Martinez P.V. (Eds.), Proceedings of the V Int. Ass. Of Hydraulic Engineering and Research, Madrid, pp. 268-273. FA N S., LA K S H M I N A R A Y A N A B. A N D BA R N E T T M. (1993) - Low-Reynolds number k – ε model for a steady turbulent boundary layers flows - AIAA Journal, 31, pp. 1777-1784. KEAN J.W., SMITH J.D. (2004) - Flow and boundary shear stress in channels with woody bank vegetation - in Bennet S.J. and Simon A. (Eds), Riparian Vegetation and Fluvial Geomorphology. American Geophysical Union, Washington, DC, pp. 237- 252. KEAN J.W., SMITH J.D. (2006) - Form drag in rivers due to small scale natural topographic features: 1) Regular sequences - J. Geophys. Res., 111, F04009, doi: 10.1029/2006JF000467. KN I G H T D., DE M E T R I O U J.D., HA M E D M.E., (1984) - Boundary shear in smooth rectangular channels - J. Hydraul. Eng. 110, pp. 405-422. JONES W.P., LAUNDER B.E. (1972) - The prediction of two 226 Fig.10a - Universal turbulent velocities comparison Confronto tra le velocità turb. adimensionali Fig.10b - Universal turb. velocities comparison in the visc. layer Confronto tra le velocità turb. adimensionali nella zona viscosa Fig.10.c - Turbulent and input velocities profile comparison Confronto tra la velocità analitica, numerica e di input A. Pasculli laminarization with a two equation model of turbu- l e n c e - Int. J. Heat Mass Transfer 1 5, pp. 301- 314. ILIOPOULOS I., MITO Y., HANRATTY T., et al., (2003) - A sto- chastic model for solid particle dispersion in a nonhomogeneous turbulent field - Int. J. of Multiphase Flow 29, 375-394 LAUNDER B.E., JONES W.P. (1972) - Mathematical Models of Turbulence - Academic Press, New York LANE S.N., BRADBROOK K.F., RICHARDS K.S., et al., (2000) - Secondary circulation cells in river channel con- fluences: Measurement artefacts or coherent flow structures - Hydrol. Process. 14, pp. 2047-2071. LI C., MOSYAK A., HETSRONI G. (1999) - Direct numerical simulation of particle-turbulence interaction - Int. J. of Multiphase Flow 25, pp. 187-200. NITHIARASU P., LIU C.B. (2006) - An artificial compressibi- lity based characteristic based split (CBS) scheme for steady and unsteady turbulent incompressible flows - Comput. Methods Appl. Mech. and Engrg. 195, pp. 2961-2962. NI T H I A R A S U P., ZI E N K I E W I C Z O.C. (2006) - Analysis of an explicit and matrix free fractional step method for incompressible flows - Comput. Methods Appl. Mech. and Engrg. 195, pp. 5537-5551. NITHIARASU P., CODINA R., ZIENKIEWICZ O.C. (2006) - The characteristic Based Split (CBS) scheme - a uni- fied approach to fluid dynamics - Int. Journal for Numerical Methods in Engrg. 66, pp. 1514-1546. PA S C U L L I A. & SC I A R R A N. (2007) - A 2D Numerical Modelling of watery debris flow including a stati- stical determination of the local erosion - 1s t American Landslide Conference, Vail, Colarado (USA) June 3-10, 2007. PASCULLI A. & SCIARRA N. (2006) - A probabilistic approa- ch to determine the local erosion of a watery debris flow - Int. Ass. For Mathematical Geology XIth International Congress. Liege September, 3th- 8th, 2006 (paper S08-08). PO S C H E L T., SC H W A G E R T. (2005) - C o m p u t a t i o n a l Granular Dynamics, Models and Algorithms - Springer Berlin Heidelberg New York. QU É M É R E P., & Sa g a u t P. (2002) - Zonal multi-domain RANS/LES simulations of turbulent flows - Int. J. Numer. Meth. Fluids 40; pp. 903-925. RINALDI M. & DARBY S.E. (2008) - Modelling river-bank- erosion processes and mass failure mechanisms: progress towards fully coupled simulations - Gravel Bed Rivers VI: From Process Understanding to River Restoration. H.Habersack, H.Piégay, M. Rinaldi Editors. Elsevier B.V. RO D I W., FE R Z I G E R J.H., BR E U R M., PO R Q U I É M. (1997) - Status of Large Eddy Simulation: Results of a W o r k s h o p - Transaction of the ASME 1 1 9, pp. 248-262. RO S A T T I G., FR A C C A R O L O L. (2006) - A well balanced approach for flows over mobile bed with high sediment transport - Journal of Computational Physics 220, pp. 312-338. SA G A U T P. (2002) - Large-Eddy Simulation for Incompressible Flows - Springer Berlin. SO U L I M., ZO L E S I O J.P. (2001) - Arbitrary Lagrangian- Eulerian and free surface methods in fluid mecha- nics - Comput. Methods Appl. Mech. Engrg. 191, pp. 451-466. SI M O N S D.B., SE N T Ü R K F. (1977) - Sediment Transport Technology - Water Resources Publications, Fort Collins, CO. SCHMEECKLE M.W. & NELSON J.M. (2003) - Direct numeri- cal simulation of bedload transport using a local, dynamic boundary condition - Sedimentology, 50, pp. 279-301. WARNER J.C., SHERWOOD C.R., ARANGO H.G., SIGNELL R.P. (2005) - Performance of four turbulence closure models implemented using a generic length scale method - Ocean Modelling 8, pp. 81-113. XU X. & HA R A D A K. (2003) - Automatic surface recon- struction with alpha-shape method - Visual Computer 19, pp. 431-443. ZI E N K I E W I C Z O.C., TA Y L O R R.L. & NI T H I A R A S U P. (2006) - Finite Element Method for Fluid Dynamics - Elsevier, 6thEdition. ZO H D I T.I. (2007) - Computational of strongly coupled multifield interaction in particle-fluid systems - Comput. Methods Appl. Mech Engrg. 1 9 6, pp. 3927-3950. 227 Ms. ricevuto il 15 marzo 2008 Testo definitivo ricevuto il 29 maggio 2008 Ms. received: March 28, 2008 Final text received: May 29, 2008 CFD-FEM 2D Modelling of a local ...