Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 42, 1, pp. 3-20, Warsaw 2004 RANDOM VORTEX METHOD FOR THREE DIMENSIONAL FLOWS. PART I: MATHEMATICAL BACKGROUND Andrzej Styczek Piotr Duszyński Marta Poćwierz Jacek Szumbarski Institute of Aeronautics and Applied Mechanics, Warsaw University of Technology e-mail: jasz@meil.pw.edu.pl The paper presents a mathematical formulation of the Lagrangian me- thod suitable for numerical simulation of 3D viscous incompressible flows. The vorticity field is approximated by a large ensemble of vor- tex particles whichmovewith the fluid (advection) and perform random walks (diffusion). The charges of the particles change with time due to the stretching term in the governing equation. The construction of the vortex particles ensures that the approximated vorticity field is stric- tly divergence-free at any time instant. The boundary condition at the surface of an immersed body is satisfied by the creation of new vortex particles near the surface. Various properties of induced velocity and vorticity fields are also discussed. Key words: vortex methods, vortex stretching, Fokker-Planck- Kolmogorow equation, Itô stochastic differential equations 1. Introduction: motivation of the paper The problem of determination of the fluidmotion past an immersed body belongs to the most important and challenging problems in fluid mechanics. Fromthemathematical pointof view, the fundamental difficulty lies in the lack of rigorous results concerning the existence and regularity of solutions to the Navier-Stokes and continuity equations in the 3D case with ”sensibly” small viscosity or, which is essentially equivalent, a sufficiently large Reynolds num- ber Re. It is commonly known that the flow in the wake behind the immersed body becomes oscillatory with relatively small Reynolds numbers, although 4 A.Styczek et al. the incoming free stream is perfectly steady. Further increase of the Reynolds number Re leads to even more complicated spatio-temporal structure of the fluidmotion and, at sufficiently large Reynolds numbers, the flow attains the fully turbulent form. Investigations of initial stages of the laminar-turbulent transition are of fundamental importance as it is believed that there lies the key to control its development. In geometrically simple cases, remarkable progress has be- en achieved using linear and/or weakly nonlinear methods of the theory of hydrodynamic stability (see Ref. Schmid and Henningson (2001)). In more complicated (and technically interesting) cases, direct numerical simulations (DNS) are inevitable. During the last twenty years, the DNS methods achie- ved a high level of sophistication (see, for instance, Peyret, 2000; Deville et al., 2002; Cottet and Koumoutsakos, 2000). In general, the methods can be divided into two classes: the Eulerian (or field) methods and Lagrangian (or particle) methods. Typical methods from the first class are finite difference methods (FDM) or finite element methods (FEM). The Lagrangian methods influidmechanics aremostly representedbyvortexmethods,where the spatio- temporal evolution of the vorticity field is simulated by a large set of discrete vorticity carriers moving as individual, mutually interacting objects. Theme- thod described in this work belongs to this class. It should be remarked that the differentiation betweenEulerian andLagrangian is not as strict as itmight seem. There is a number of approaches which are in a sense ”hybrid” to name only the Large Eddy Simulation (LES) (Sagaut, 2002), methods with moving and/or deforming grids [15] or the arbitrary Lagrangian-Eulerian formulation applied to flows with moving and/or free boundaries (Deville et al., 2002). While the central idea in the Eulerian approach is a grid (or mesh), the Lagrangian approach is based rather on the usage of moving, spatially locali- zed individuals called particles. Theparticles are fictitious creatures – they are carriers of the vorticity or magnetization (see Buttke and Chorin, 1993; Sty- czek and Szumbarski, 2002). The method using vortex particles proved to be especially effective in 2D simulations – in this context it is usually referred to as themethod of vortex blobs (see, for instance, Błażewicz and Styczek, 1993; Szumbarski and Styczek, 1997; Styczek and Wald, 1995; Protas and Styczek, 2002). In such a case, the vorticity has only one nonzero component subject only to advection and diffusion (there is no stretching). The generalization for axisymmetric flows was formulated in Hedar and Styczek (1999). Since the governing equation can be cast into the form which formally resembles the 2D case, the axisymmetric vortex method does not differ much from the two dimensional one. For general 3D flows, the situation is much more complex – Random vortex method for 3D flows. Part I 5 there is an additional coupling between the vorticity and local variations of the velocity field. The effect of this coupling is, generically, the stretching of the vorticity due to local velocity gradients. In the region of the wake behind the immersed body, the velocity field is characterized by immense fine-scale non-uniformity. The velocity-vorticity coupling gives rise to specific synergism,which efficiently breaks upflow struc- tures into even smaller scales. Thus, there is a need for subtle discretization of the wake region, where the spatial scale of the velocity variations is small and time-dependent. The method of the vortex particles meets this demand in an apparently natural way: the vortex particles are born in vicinity of the body, shed from it and gathered in the wake area. This way, higher resolution in the regions with reach kinematics is attained. The other advantage of the vortex method is the elimination of the pressure. It should be remarked that other ways of the pressure elimination (not based on using the curl operator) lead to a complicated integro-differential problem. The paper is constructed as follows. In Section 2, we formulate the initial boundary problem and transform it to the velocity-vorticity formulation. We derive also some important results concerning the properties of the vorticity in the flow domain. In Section 3, we explain the construction of the 3D vortex particles. In Section 4, the velocity field is constructed and its properties are discussed in some details. In Section5, the splitting of the governing equation into advection/diffusion and the stretching parts are introduced. Next, we explain how the partial problems are solved during each time step of the method. In Section 6, the realization of the no-slip boundary conditions and creation of the vorticity at thematerial boundary in the form of new particles is considered. 2. The velocity-vorticity formulation of the flow problem Consider a nonstationary flowof a viscous fluid in the exterior of amotion- less body with the smooth surface Γ . The unknowns are the velocity field v and the pressure p. The specificmass is constant and equals to the unity. The fluidmotion is governed by the Navier-Stokes equation ∂tv+(v ·∇)v=−∇p+ν∆v (2.1) The velocity field is divergent-free ∇·v=0 (2.2) 6 A.Styczek et al. and satisfies the following boundary conditions v ∣ ∣ ∣ r∈Γ =0 v ∣ ∣ ∣ r→∞ =U∞ (2.3) At the initial time, the velocity is given as v ∣ ∣ ∣ t=0 =v0(r) (2.4) The vector field v0 satisfies conditions (2.2) and (2.3)1, and asymptotic condi- tion (2.3)2 at t= t0. We also assume that at infinity the pressure is spatially homogeneous. As it was already mentioned, the proof of existence and uniqueness of a regular (smooth) solution to the above problem is not yet available, unless some additional restrictions (symmetry, sufficiently small Reynolds number) are assumed. We will show that the pressure p is a nonlocal quantity. Indeed, applying the divergence operator to equation (2.1) one gets −∆p= ∂vi ∂xj ∂vj ∂xi (2.5) In the above, the summation convention has been applied. The following Neumann condition can be formulated at the surface Γ n ·∇p ∣ ∣ ∣ Γ = νn∆v≡−νn(∇×ω) ∣ ∣ ∣ Γ (2.6) The vorticity ω is defined as the curl of the velocity ω=∇×v (2.7) Bydefinition, the vorticity has zerodivergence.Poisson equation (2.5)with boundary condition (2.6) can be formally invertedwith the use of the generali- zed Green function. Then, inserting the pressure gradient into equation (2.1), one obtains the previouslymentioned nonlocal, integro-differential formof the governing equation. The other method of the pressure elimination consists in using the equation of vorticity transport (the Helmholtz equation) ∂tω+(v ·∇)ω= ν∆ω+(ω ·∇)v (2.8) The velocity field v can be then expressed as follows v=∇(rU∞(t)+ϕ)+ 1 4π ∫ suppω ω(t,ξ)× r−ξ |r−ξ|3 dξ (2.9) Random vortex method for 3D flows. Part I 7 where the potential ϕ is a harmonic function vanishing at infinity and the symbol suppω denotes the vorticity support (the closure of the subset in R3, where ω 6=0). The initial condition for the vorticity field is implied by condition (2.4) ω ∣ ∣ ∣ t=0 =ω0(r)≡∇×v0(r) (2.10) We also formulate the Neumann boundary condition − ∂ϕ ∂n ∣ ∣ ∣ Γ =nU∞+ n 4π ∫ suppω ω(t,ξ)× r−ξ |r−ξ|3 dξ (2.11) which – formally – allows for determination of the harmonic potential ϕ. Next, taking the component of the velocity v tangent to the surface Γ , we get n×(n×U∞(t)+n×∇ϕ) ∣ ∣ ∣ Γ +n× {n 4π × ∫ suppω ω(t,ξ)× r−ξ |r−ξ|3 dξ } ∣ ∣ ∣ Γ =0 (2.12) The nonlocality of formula (2.11) and (2.12) is evident. In a sense, it has been ”shifted” from the equation of motion. However, we will show later that con- ditions (2.11) and (2.12) are relatively easy and straightforward to implement in the method of vortex particles. Thevorticityω subjects toanadditional restriction.Letus integrateHelm- holtz equation (2.8) in the fixed region D bounded by thematerial surface Γ and the sphere S∞ of the very large radius R. d dt ∫ D ω dr+ ∫ Γ∪S∞ (nv)ω dS=−ν ∫ Γ∪S∞ n× (∇×ω) dS+ ∫ Γ∪S∞ (nω)v dS Since v|Γ = 0 and the vorticity vanishes at S∞ faster than r−2, one obtains for R→∞ d dt ∫ ω dr=−ν ∫ Γ n× (∇×ω) dS Using equation of motion (2.1) written for r→ rΓ ∈Γ , we get ∇p ∣ ∣ ∣ Γ =−ν∇×ω ∣ ∣ ∣ Γ Thus, the equality can be inferred −ν ∫ Γ n× (∇×ω) dS = ∫ ∇× (∇p) dr≡0 8 A.Styczek et al. which immediately leads to the following conclusion d dt ∫ suppω ω(t,r) dr=0 (2.13) One can see that the total charge of the vorticity is conserved in the presence of a motionless material boundary (the surface Γ) of the flow domain. 3. The particles of vorticity The vorticity field with a bounded support induces the velocity field with the asymptotic behavior v=C× r r3 +O(r−3) Such behavior can be concluded from the form of the Biot-Savart integral appearing in formula (2.9). The vector constant C is proportional to the charge of vorticity.Wewill generalize the above expression by postulating the velocity field induced by a single vortex particle placed at the origin in the following form v= Ω×r r3 F(r)=Ω×∇Φ=−∇× (ΩΦ) (3.1) The time-dependent vector Ω(t) is an individual characteristic of a particle. and the potential Φ is spherically symmetric and defined as Φ(r)= r ∫ 0 F(ξ) ξ2 dξ (3.2) The corresponding vorticity field is given as ω=Ω∆Φ−∇(Ω ·∇Φ) (3.3) The internal structure of the vortex particle is implied by the determination of the Laplacian of Φ ∆Φ≡ 1 r2 d dr ( r2 dΦ dr ) = f(r) (3.4) Random vortex method for 3D flows. Part I 9 If the following function is introduced F(r)= r ∫ 0 ξ2f(ξ) dξ (3.5) the vorticity field of the particle can be expressed in the following manner ω=Ω [ f(r)− 1 r3 F(r) ] − (Ωr)r r2 [ f(r)− 3 r3 F(r) ] (3.6) The total charge of vorticity in the particle can be calculated as Qω = ∫ ω dξ= ∞ ∫ 0 ∫ S(ξ) ω(ξ) dS(ξ)dξ The iterated integration is carried out over the sphere S(ξ) with the radius equal to ξ, and next along the variable ξ. The internal integration can be performed as follows ∫ S(ξ) (Ωξ)ξ ξ2 dS(ξ)= ∫ S(1) (Ωζ)ζ dS(1)= eβ ∫ S(1) Ωαζαζβ dS(1)= = 4π 3 eβδβαΩα = 4π 3 Ω In the above, the integrals for α 6= β vanish identically, and the remaining three integrals are equal. Since their sum is simply the surface of the unitary sphere, each integral is equal to 4π/3. Using this fact, the calculation of the integral with respect to ζ finally gives Qω = 8π 3 Ω ∞ ∫ 0 ξ2f(ξ) dξ= 8π 3 ΩF(∞) (3.7) Note that if F(r) vanishes for r > r0, the vortex particle carries zero total vorticity charge, and there is no induction in this area.Wewill postulate that the function f(r) vanishes for r > σ. The quantity σ will be referred to (conventionally) as the particle radius. The function f(r) can be unbounded for r → 0. However, to avoid singularity in the velocity field at the particle center, the function f(r) should satisfy the inequality |f(r)|< const ·r−αwith α < 1. In such a case v(0) = 0, and, additionally, the value of the potential Φ(r) exists in the particle center. If the function f vanishes for r>σ, then F(∞)= σ ∫ 0 ξ2f(ξ) dξ=F(σ) (3.8) 10 A.Styczek et al. We will normalize the function f so that F(∞) = 1. For r > σ, the poten- tial Φ can be expressed as follows Φ(r)= σ ∫ 0 F(ξ) ξ2 dξ+ r ∫ σ 1 ξ2 dξ=Φ(σ)+ 1 σ − 1 r =Φ(∞)− 1 r (3.9) Note also that the induced vorticity field outside the particle ”core” (the inte- rior of the sphere with the radius σ) is the potential vector field. Indeed, one can write ω=∇× (Ω×r r3 ) =−∇Ωr r3 =∇ ( Ω ·∇1 r ) (3.10) For a particularly simple function f(r), namely f(r)=      3 σ3 for r <σ 0 for r >σ (3.11) we get v=Ω×      r σ3 for r¬σ r r3 for r >σ (3.12) and ω=        2Ω σ3 ≡∇ ( 2 Ωr σ3 ) for r <σ −Ω r3 +3 (Ωr)r r5 ≡−∇ (Ωr r3 ) for r >σ (3.13) The inducedvelocity field is continuous in space.Thevorticity field is constant inside the particle core anddiscontinuous across the core boundary r=σ. The jump across the discontinuity [ω] is tangent to the core boundary (the normal component is continuous) and can be evaluated as follows [ω] =ω ∣ ∣ ∣ + −ω ∣ ∣ ∣ − = 3 σ3 [(Ωr)r σ2 −Ω ] (3.14) If the function f(r) is continuous then the tangent component of the vorticity is continuous aswell. If we ”mollify” the jumpattained by step function (3.11) at r = σ, we obtain the vorticity field with a rapid variation in vicinity of r = σ. Note that zero-divergence constraint imposed on the field ω implies that any vortex particle which induces in thewhole spacemust have an unbo- unded vorticity support. In such a case, the total charge of the particle Qω Random vortex method for 3D flows. Part I 11 definedby (3.7) is always different fromzero. If the total charge of an ensemble of the vortex particles is zero then, for a sufficiently large radius r (reaching beyond all particle cores in the ensemble), the induced velocity and vorticity fields will decay at a rate proportional to r−3 and r−4, respectively. 4. Properties of the velocity field The velocity field defined by expression (2.9) is divergence-free. Similarly, the divergence of the vorticity induced component determined by the Biot- Savart integral is zero. In can be concluded that the following equality holds for any closed surface immersed in the flow domain ∫ n ·∇ϕdS≡ ∫ ∂ϕ ∂n dS =0 (4.1) which ensures correctness of Neumann boundary problem (2.11). Using the Green formula, the harmonic function ϕ can be expressed as follows 4πϕ(r)= ∫ Γ ( ϕ(ξ) ∂ ∂nξ 1 |r−ξ| − ∂ϕ ∂nξ 1 |r−ξ| ) dSΓ (4.2) The expansion of the integral kernels and use of equality (4.1) yield the follo- wing asymptotic behavior with r→∞ ϕ(r)= const r2 +O(r−3) (4.3) Thevelocity field,which is the sumof theuniformstream,potential component ∇ϕ and the component induced by the set of vortex particles, can be written in the form v=U∞+∇ϕ+ N ∑ k=1 Ωk× (r−rk) |r−rk|3 F(|r−rk|) (4.4) In the above, rk denotes the location of the center of the kth particle (k=1, ..,N). The asymptotic form of the velocity field for r →∞ can be obtained by expanding the denominator |r−rk|−3 and using (4.3). The resulting form is v = U∞+ r r4 const + ( ∑ Ωk)×r r3 − ∑ (Ωk×rk) r3 + (4.5) + 3 ∑ [Ωk×rk(rrk)] r5 +O(r−4) 12 A.Styczek et al. The velocity field, (4.5), has to be compatible with the constraint imposed on the total vorticity charge, see (2.13). In view of (3.7), we conclude that the following equality must hold Ω0 = N ∑ k=1 Ωk = const If the velocity vanishes at the initial time instant, then obviously Ω0 = N ∑ k=1 Ωk =0 (4.6) The zero value of the total vorticity chargemeans that the total amount of the angularmomentum in the flow is finite. To show this, let us calculate the total angularmomentumof the fluid outside the ball containing all vortex particles. In the course of the calculations we will obtain the following integrals ∫ |r|>R r× [Ωk× (r−rk)] |r−rk|3 dr= ∫ |ξ|>R ξ× (Ωk×ξ) ξ3 dξ+rk× ∫ |ξ|>R Ωk×ξ ξ3 dξ The second integral in the right-hand side is zero. The first one can be trans- formed to the form 8π 3 Ωk ∫ |ξ|>R ξ dξ Summing up the contributions fromall vortex particles one concludes that the total amount of the angularmomentumwill remain finite only if condition (4.6) is satisfied. Note that the potential component can be transformed to the surface integral ∫ Γ n×rϕ dSΓ , which is also finite and therefore cannot prevent infinity in the momentumwhen condition (4.6) is violated. If the restriction was not imposed, an infinite amount of the angular vorti- citywouldbeproduced in theflowwithin afinite time.This singularity cannot be avoided in a way similar to how it is done in the case of the infinite mo- mentum.The latter difficulty can be removed simplyby choosing the reference frame related to themotionless fluid at infinity – in such a case the immersed bodymoves with respect to the ambient fluid. It is reasonable to use the ”corrected” form of the velocity v=U∞+∇ϕ+ N ∑ k=1 Ωk× (r−rk) |r−rk|3 F(|r−rk|)− Ω0×r r3 (4.7) Random vortex method for 3D flows. Part I 13 and the vorticity ω= N ∑ k=1 ωk(t,r−rk)−ω0(t,r) (4.8) where ωk(t,ξ)=Ωk(t) [ f(ξ)− 1 ξ3 F(ξ) ] − (Ωkξ)ξ ξ2 [ f(ξ)− 3 ξ3 F(ξ) ] In the above, the symbol ω0 denotes the vorticity field induced by the ad- ditional vortex particle located inside the immersed body (outside the flow domain). If (4.6) holds exactly, both velocity andvorticity remain unchanged.Howe- ver, in the course of numerical computation, inaccuracies are inevitable. Then, the correction terms in (4.7) and (4.8) ensure that condition (4.6) is verified. 5. Evolution of the vorticity field According to the Lie-Trotter Theorem (see Chorin and Marsden, 1997), we split Helmholtz equation (2.8) into the equation of advection/diffusion ∂tω+(v ·∇)ω= ν∆ω (5.1) and the ”stretching” equation ∂tω=(ω ·∇)v (5.2) The evolution operator corresponding to the Helmholtz equation can be expressed by a formal composition of the resolvents of equations (5.1) and (5.2) Ht = lim n→∞ (St/n ◦At/n)n (5.3) The emergence of the vorticity at material boundaries of the flow domain introduces an additional element in composition (5.3) – the operator of the vorticity creation Ψt (Chorin andMarsden, 1997) Ht = lim n→∞ (St/n ◦Ψt/n ◦At/n)n (5.4) 14 A.Styczek et al. Wewill determine the resolvent operators used in (5.4). Equation (5.1) for- mally resembles theFokker-Planck-Kolmogorov equation (seeGardiner, 1990), governing the spatio-temporal evolution of the density of the conditional pro- bability p(t,r|t0,r0) for the stochastic process determined by the Itô differen- tial equations dri =vi dt+ √ 2ν dW i ri ∣ ∣ ∣ t=t0 = ri0 (5.5) In the above, vi denotes an instantaneous velocity at the location of the ith ”diffusive particle”, while dW i denotes the infinitesimal increment of the vector Wiener process modeling the diffusion of the particles. In view of this interpretation, the vorticity field at the time instant t > t0 can be calculated as a conditional expectation as follows ω(t,r)= ∫ suppω0 p(t,r|t0,r0)ω0(r0) dr0 (5.6) The vorticity field determined by (5.6) satisfies equation (5.1) and the initial condition ω ∣ ∣ ∣ t=t0 =ω0(r) Equation (5.2) can be re-written using directional differentiation as follows ∂tω=(ω ·∇)v≡ d dλ v(t,r+λω) ∣ ∣ ∣ λ=0 (5.7) For a given velocity field, formula (5.7) defines a parameterized family of or- dinary differential equations. The parameters of this family are spatial coor- dinates (or, simply, the vector r). The resolvent operator for equation (5.2) is defined by the transformation ω(t,r) → ω(t+∆t,r) (r plays the role of the parameter), which is simply the process of numerical integration of the ordinary differential equations. Finally, the operator of the vorticity creation for the time period ∆t de- scribes production of such a boundary distribution of the vorticity so that the boundary condition for the velocity field (no-slip condition) is satisfied. An instantaneous vorticity distribution at the boundary is related via certain in- tegral operation to the vorticity field in the whole flow domain. If the particle representation of the vorticity field is used the integral relation is, as we will show in the next Section, approximated by an algebraic one. Random vortex method for 3D flows. Part I 15 Summarizing, the vorticity field at the time instant t+∆t is determined due to: • Modification of the vorticity field at the instant t caused by the vortex particles displacement during the time interval (t,t+∆t), accordingly to Itô equations (5.5) • Transformation of the vorticity field by integrating ordinary differential equations (5.7) describing the effect of the ”stretching” • Creation of the vorticity at the material boundary so that condition (2.12) is verified. 6. Vortex particle approximation of the vorticity field Assume that we have N vortex particles in the flow domain at the time instant t. The vortex charges ΩOk and coordinate vectors r O k (k = 1, ...,N) of all particles are known. During a short time interval ∆t the particles move accordingly with the Itô equations ∆r0k =vk∆t+ √ 2ν ∆Wk (6.1) and their vorticity charges change due to the stretching effect. Since at the particle central point, i.e. for r= rOk , the vorticity carried by the particle is ω ∣ ∣ ∣ 0 = 2 3 f(0)Ω (6.2) then the stretching equation takes the form dΩ0k dt = 3 2f(0) d dλ v ( t,r+λω(t,rk) ) ∣ ∣ ∣ λ=0 (6.3) and a new vorticity charge of the particle can be evaluated Ω0k(t+∆t)=Ω 0 k(t)+∆Ω 0 k (6.4) Since the vortex particles change both their positions and vorticity charges, condition (2.3)1 imposed on the velocity field is violated at the time instant t+∆t. In order to ensure (2.3)1, it is necessary to modify the boundary (or near-boundary) value of the vorticity field and calculate a new potential ϕ. The former can be achieved by introducing new vortex particles in a closed 16 A.Styczek et al. vicinity of the surface of the immersed body. The potential ϕ is sought in the form of the sum ϕ=ϕU +ϕO+ϕ ∗ (6.5) Each term in (6.5) represents a harmonic function, specified by appropria- te boundary conditions of the Neumann type. These conditions for two first functions in (6.5) are ∂ϕU ∂n ∣ ∣ ∣ Γ =−nU∞ (6.6) ∂ϕO ∂n ∣ ∣ ∣ Γ =n ( N ∑ k=1 Ω 0 k ) × r r3 F(r)−n N ∑ k=1 Ω 0 k× r−rk |r−rk|3 F(|r−rk|) It can be verified that both (6.6) verify the solvability condition of the Neu- mann boundary problem, and hence the functions ϕU and ϕO can be found using standardmethods. The third harmonic potential ϕ∗ is connected to newly created vortex particles.Wewill assume that M such particles are born at the fixed positions rnewk near thematerial surface at each time step. In general, each new particle induces the velocitywithnonzero tangent andnormal components at the body surface. The role of the potential ϕ∗ is to compensate the induction normal to the surface incurred by the newvortex particles. To achieve this, the potential is defined in the form of the sum ϕ∗ = M ∑ k=1 (Ωnewk1 ϕk1+Ω new k2 ϕk2+Ω new k3 ϕk3) (6.7) where ϕkα, α= 1,2,3, are harmonic functions fulfilling the Neumann boun- dary conditions as follows ∂ϕkα ∂n ∣ ∣ ∣ Γ =n eα×r r3 F(r)−neα× (r−r new k ) |r−rnewk |3 F(|r−rnewk |) ∣ ∣ ∣ r∈Γ (6.8) Using the same standard procedure as in the case of ϕU and ϕO, we can cal- culate 3M potentials related individually to all newly created vortex particles. The velocity field induced by these particles is given as v new = M ∑ k=1 3 ∑ α=1 Ωnewkα [ ∇ϕkα− eα×r r3 F(r)+ eα× (r−rnewk ) |r−rnewk |3 F(|r−rnewk |) ] (6.9) Random vortex method for 3D flows. Part I 17 Note that if the places of creation of the new particles are fixed, all poten- tials ϕkα, α = 1,2,3 can be calculated once and forever. For brevity, let us write formula (6.9) as vnew = M ∑ k=1 K̂(r,rnewk )Ωk (6.10) Note that the following equality holds nK̂(r,rnewk ) ∣ ∣ ∣ r∈Γ =0 (6.11) Inorder tofindthe M vector quantities Ωnewk (or,which is equivalent, 3M sca- lar values Ωnewkα ), we write the expression for the velocity as follows v=U∞+∇ϕU +∇ϕO+ N ∑ k=1 Ωoldk × (r−rk) |r−rk|3 F(|r−rk|)− (6.12) − N ∑ k=1 Ωoldk × r r3 F(r)+ M ∑ k=1 K̂(r,rnewk )Ω new k On the surface Γ the normal component of the velocity v is zero identically. One has to eliminate also the tangent component. To this end, we multiple equation (6.12), written for r ∈ Γ , by the orthogonal tangent versors ηβj, β = 1,2, and enforce the no-slip condition at the collocation points rΓj ∈ Γ , j = 1, ...,M. Moving all terms with the known quantities to the right-hand sides, we obtain an algebraic linear system of 2M equations M ∑ k=1 ηβjK̂(r Γ j ,r new k )Ω new k =−ηβj [ U∞+∇ϕU +∇ϕO+ (6.13) + N ∑ k=1 Ωoldk × (rΓj −rk) |rΓj −rk|3 F(|rΓj −rk|)− N ∑ k=1 Ω old k × rΓj r3 F(rΓj ) ] In order to get a solvable system of equations, we postulate M additional conditions in the form of M ∑ k=1 nΩnewk ∣ ∣ ∣ r=rΓ k =0 (6.14) It can be shown that thematrix of the full system is not singular.Moreover, if the vectors rnewk and r Γ k , k=1, ...,M, are fixed, than the coefficientmatrix of the obtained linear systemcanbe evaluated andLU-factored once and forever. 18 A.Styczek et al. 7. Final remarks The three dimensional generalization of the Vortex Blobs Method, for- mulated originally for 2D fluid motion, has been presented. The dynamics of vorticity in three dimensions is essentially more complex than in the 2D case, which makes the construction of the vortex method a nontrivial task. Three dimensional vorticity fields must obey restriction (2.13). The governing equ- ation contains the stretching term giving rise to additional variations of the vorticity due to local nonuniformity of the velocity field, (5.2). Therefore, the 3D vortex particles (analogues of two dimensional vortex blobs) evolve indi- vidually in time, while the total charge of the vorticity is constant in time. In two dimensional motion, there is a restriction imposed on the pressure field (corresponding to the conservation of the total vorticity), but the stretching effect is not involved. The proposed method of accounting for the stretching term is not the only possible. Nevertheless, it is algorithmically simple and seems to be ap- propriate for the vortex particles with small cores with highly concentrated vorticity. One of the main problems with the vortex methods is how to maintain computational efficiency during long-time simulations. Note that the new vor- tex particles are continually generated at the material boundaries of the flow domain, and thus their number gets larger from one step to another. In effect, the CPU time required to evaluate the velocity field at each time step grows in the course of computations. Known algorithms for fast summation can re- duce the numerical cost, however, not as spectacularly as in the 2D case. In addition, they are too complex to implement. The other method of cost re- duction is to put bounds on the number of particles in the flowdomain, which is achieved by selective removal and/or gluing. To this end various strategies can be proposed, but their reasonable application is limited by the lack of ri- gorous relation between the number of particles and the approximation order. Since these issues are still open, the evaluation of the results consists mostly in comparison with the experimental evidence (if available) and with other methods. The implementation details of the proposed vortexmethod and discussion of the obtained results will be presented in Part II of the paper. Acknowledgements This work has been supported by the State Committee for Scientific Research (KBN), grant No. 8TO7A02520. Random vortex method for 3D flows. Part I 19 References 1. Błażewicz J., Styczek A., 1993, The stochastic simulation of a viscous flow past an airfoil, J. Theor. Appl. Mech., Part 1: 31, 1; Part 2: 31, 4 2. Buttke T.F., Chorin A.J., 1993, Turbulence calculations in magnetization variables,Appl. Numer. Methods, 47, 12 3. Chorin A.J., Marsden J.E., 1997, A Mathematical Introduction to Fluid Mechanics, 3rd Ed. Vol.4 of Texts in Applied Mathematics, Springer Verlag, NewYork 4. Cottet G.H., Koumoutsakos P., 2000,Vortex Methods: Theory and Prac- tice, Cambridge University Press 5. Deville M.O., Fisher P.F., Mund E.H., High-Order Methods for Incom- pressible Fluid Flow, Cambridge University Press 6. Gardiner C.W., 1990, Handbook of Stochastic Methods, 3rd Ed. Springer Verlag 7. Hedar I.M., Styczek A., 1999, Random Vortex Method approach to axi- symmetric jet in a large tank, J. Theor. Appl. Mech., 37, 4, 863-872 8. PeyretR. (Ed.), 2000,Handbook of theComputational FluidDynamics,Aca- demic Press 9. Protas B., Styczek A., 2002, Optimal rotary control of the cylinder wake in the laminar regime,Physics of Fluids, 14, 7, 2073-2087 10. Sagaut P., 2002,Large Eddy Simulation for Incompressible Flows, An Intro- duction, Springer Verlag, 2nd Ed. 11. Schmid P.J., Henningson D.S., 2001, Stability and Transition in Shear Flows, AMS, 142, Springer Verlag, NewYork 12. Styczek A., Szumbarski J., 2002, On the magnetization-based Lagrangian methods for 2Dand3Dviscousflows.Part I:TheoreticalBackground,J.Theor. Appl. Mech., 40, 2, 339-356 13. StyczekA.,WaldP., 1995,Fast and efficientVortex-Blobs simulation of the flow past a circular cylinder,Arch. Mech. Eng., XLII, 3-4 14. Szumbarski J., Styczek A., 1997, The stochastic vortexmethod for viscous incompressible flow in a spatially periodic domain,Arch. Mech., 49, 1, 209-232 15. TheWWWpage of Guojun Liao at www.uta.edu/math/faculty/liao/pages/index.htm 20 A.Styczek et al. Stochastyczna metoda wirowa dla przepływów trójwymiarowych. Część I: Podstawy matematyczne Streszczenie W pracy sformułowano lagranżowskąmetodę wirową dla trójwymiarowych prze- pływów cieczy lepkiej. Pole wirowści jest w niej aproksymowane zbiorem cząstek po- ruszających się wraz z płynem oraz wykonujących ruch losowy (dyfyzja). Ładunki wirowości niesione przez cząstki zależą od czasu (efekt ”strechingu”). Konstrukcja cząstek zapewnia, że pole wirowości pozostaje bezźródłowe w trakcie symulacji. Wa- runek dla pola prędkości stawiany na powierzchni opływanego ciała jest realizowany drogą generacji nowych cząstek w sąsiedztwie tej powierzchni.Wpracy przedyskuto- wano również własności indukowanego pola prędkości i wirowości. Manuscript received September 15, 2003; accepted for print November 11, 2003