Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 2, 40, 2002 ON THE MAGNETIZATION-BASED LAGRANGIAN METHODS FOR 2D AND 3D VISCOUS FLOWS. PART 1 – THEORETICAL BACKGROUND Andrzej Styczek Jacek Szumbarski Institute of Aeronautics and Applied Mechanics, Warsaw University of Technology e-mail: jasz@meil.pw.edu.pl The paper presents the background of an alternative formulation of the Navier-Stokes equationusingavariable called themagnetization.Several variants of governing equations, based ondifferent choices of a particular gauge transform, are discussed.The remaining part of the paper is devo- ted to the formulation of a Lagrangian approach to 2D and 3D viscous flows. First, the carrier of the magnetization (the magneton) is defined and the corresponding induction law is derived. The instantaneous velo- city field is constructed as a superposition of contributions from a large set of magnetons and a uniform stream. An essential feature of the me- thod is a one-time-step operator splitting, consisting in the consecutive solution of three sub-problems: generation of themagnetization on solid boundaries, advection-diffusion of the magnetization and stretching. Key words: viscous flow, Navier-Stokes equations, magnetization, gauge transform 1. Introduction Viscous incompressible flows are usually describedby the set of theNavier- Stokes and continuity equations ∂tv+(v ·∇)v=−∇p+ν∆v (1.1) ∇·v=0 Theunknownsare the velocity field v and thepressure p.Thefluiddensity is assumedhere to beunary. It shouldbenoted that thepressure is determined modulo an arbitrary function of time. 340 A.Styczek, J.Szumbarski According to a well-known theorem, any vector field can be expressed as a sum of two components. One of these components is defined as a gradient of a scalar field, while the second one is a solenoid vector field. Thus, the divergence-free velocity v can be written as v=m−∇φ (1.2) The vector field m is called, after Buttke (Chorin, 1994), the ”magnetiza- tion”. The correspondence m → (v,φ) is unique, providing that the normal component of the velocity v at the boundary is equal to zero, and appro- priate regularity conditions are satisfied. On the other hand, the transform v→ (m,φ) obtained from (1.2) is not unique.This fact allows us to formulate the governing equations in alternative forms, employing variousmagnetization fields. In the paper we show a few of them, discuss their properties, formula- te a Lagrangian numerical method for the magnetization and finally present results of some numerical experiments. 2. The gauge transform Usingvelocity representation (1.2), theNavier-Stokes equation canbewrit- ten as ∂tm+(v ·∇)v= ν∆m+∇(∂tφ−ν∆φ−p)= 0 (2.1) Equation (2.1) can be viewed as a particular case of a more general equation ∂tm+(v ·∇)v= ν∆m+∇λ (2.2) where λ denotes some scalar field called a gauge field. We will refer to any particular selection of the gauge field λ as a gauge transform. Let us define the subspace of solenoidal vector fields V = {v∈L2(Ω) : ∇·v=0 in Ω, v ·n=0 at ∂Ω} Any square integrable vector field in Ω can be written as a sum of some element from the space V and the gradient of a certain scalar field. This decomposition is unique and the components are orthogonal with respect to the inner product in L2(Ω). Thus, the orthogonal projection operator can be defined m→Pm≡v∈V On the magnetization-based Lagrangian methods... 341 and equation (2.2) can be formulated as follows ∂tm+(Pm ·∇)Pm= ν∆m+∇λ (2.3) Assume that some gauge field λ has been chosen and the magnetization field obeys equation (2.3). Then the solenoidal vector field v =Pm satisfies the governing equations of a viscous liquid motion. Indeed, the continuity equation is fulfilled since v ∈ V . The magnetization can be expressed as m = v+∇φ. Insertion of this form into (2.3) yields after some algebra the following equation ∂tv+(v ·∇)v= ν∆v+∇(λ−∂tφ+ν∆φ) (2.4) Thus, the vector field v satisfies the Navier-Stokes equation, and the corre- sponding pressure is given as p= ∂tφ+ν∆φ−λ+f(t) (2.5) It can be shown that the gauge field λ can be chosen arbitrarily. Assume there are two differentmagnetization fields m1 and m2 corresponding to the gauge fields λ1 and λ2, respectively. The field m1 can be expressed as a sum m1 =v+∇φ1. It is immediate to verify that equation (2.3) for λ=λ2 admits the solution m2 in the following form m2 =v+∇φ2 Thus, the solution m2 defines the same velocity field and can be also verified that ∇(∂tφ1−ν∆φ1−λ1)=∇(∂tφ2−ν∆φ2−λ2) thus the resulting pressure fields differ only by a (time-dependent) constant. If an initial/boundary-value problem formulated for equation (2.3) permits a unique solution, then one concludes that themagnetization fields computed for different gauges will correspond to the equivalent velocity and pressure fields. Thus, the choice of a particular gauge is arbitrary. Note that the non- linear term in equation (2.3) is a bilinear one. Looking for convenient variants of equation (2.3), we assume the following form of the gauge field λ λ= 1 2 αv ·v+βv ·m+ 1 2 γm ·m (2.6) After elementary calculations, we find the gradient of the above expression ∇λ = α[(v ·∇)v+v× rotv]+γ[(m ·∇)m+m× rotm] (2.7) + β[(v ·∇)m+(m ·∇)v+v× rotm+m× rotv] 342 A.Styczek, J.Szumbarski In (2.7) α, β, γ are arbitrary constants. Any particular choice of these constantsdefines somegauge transform,which, in turn,determinesagoverning equation for the magnetization field. Having it solved one is able to find the potential φ, and then the pressure field can be determined with the use of (2.5). Here is a list of some particular, possibly interesting variants of the gauge transform and the corresponding form of the governing equation: Case 1: α= γ=1, β=−1 ∂tm+(v ·∇)m= ν∆m+(m ·∇)∇φ (2.8) p+ 1 2 v2 = f(t)+∂tφ−ν∆φ− 1 2 m ·m+m ·v Case 2: α=1, β=−1, γ =0 ∂tm+(v ·∇)m= ν∆m− (∇v)⊤m (2.9) p+ 1 2 v2 = f(t)+∂tφ−ν∆φ+m ·v Case 3: α=1, β=0, γ=−1 ∂tm+(m ·∇)m= ν∆m−∇φ× rotm (2.10) p+ 1 2 v2 = f(t)+∂tφ−ν∆φ+ 1 2 m ·m Case 4: α=1, β= γ =0 ∂tm−v× rotm= ν∆m (2.11) p+ 1 2 v2 = f(t)+∂tφ−ν∆φ Case 5: α=β= γ =0, ∂tm+(v ·∇)v= ν∆m (2.12) p= f(t)+∂tφ−ν∆φ On the magnetization-based Lagrangian methods... 343 Case 2 has been introduced by Buttke (Chorin, 1994), Russo and Smerka (1999) analysed cases 2, 4 and 5. Summers andChorin (1996)made use of the equations in case 2 and formulated the Lagrangian method for particles car- rying magnetization. Generally, the Buttke gauge seems to be most popular. Russo and Smerka remarked that the stretching term in the first equation in (2.9), i.e. (∇v)⊤m is similar to the vortex stretching term in the Helmholtz equation for the vorticity. This term describes themechanism of the deforma- tion of vortex structures in the flow, consisting in local stretching in one space direction and compression in the remaining two.This behavior is related to the existence of a large positive (or negative) eigenvalue of the tensor ∇v. The dominating eigenvalue of ∇v renders the magnetization (or vorticity) tend to be concentrated into the form of thin filaments oriented locally along the corresponding eigendirection. Changing the gauge onemay be able to reduce this phenomenon. It seems that the gauge defined in case 1 is a better choice. The stretching term in this case has the form (m ·∇)∇φ=mH(φ) where H(φ) is the Hessian, i.e. H(φ) = {∂xixjφ}. The trace of this matrix, equal to the sum of the eigenvalues, is not restricted by the incompressibili- ty condition. This fact can reduce the tendency for high concentration and instability. In this study, the formulation with inhomogeneous terms with respect to the magnetization will not be analysed. Any inhomogeneous term in the go- verning equation acts as a source of themagnetization, i.e. themagnetization is created at points where it was previously equal to zero. It seems to be unreasonable to formulate a Lagrangian (or particle) approach based on such formulation. In addition, the case 4 is not convenient for numerical compu- tations. According to Russo and Smerka (1999) this formulation suffers from essential instability which may be difficult to overcome. 3. The carrier of the magnetization In order to construct a numerical method based on the Lagrangian appro- ach one should define an elemental ”source” (or particle) of themagnetization field. It is convenient (and natural) to assume spherical symmetry of the spa- tial distribution of the magnetization insidesuch a particle. The radius of the particle can befinite – in such a case themagnetization differs fromzero inside 344 A.Styczek, J.Szumbarski the spherical ”core” and vanishes identically outside the core. If the radius is infinite, themagnetization fills thewhole space. The spatial distribution of the magnetization shouldbe, however, concentrated in a relatively small neighbor- hood of the particle center. The effect of the ”localization” can be achieved by applying, for instance, a ”slim” Gaussian distribution. In the remaining part of this work, we will refer to the magnetization particles (independentlyon thedetails of their construction) as themagnetons. Themagnetization field carried by a single magneton is defined as follows m=m0(t)g(r) (3.1) where r denotes the distance from the magneton center and g is a regular core function. The vector ”charge” of themagneton is, in general, time-dependent and it is characterized by the scaling factor m0(t). It should be clear that the potential term in decomposition (1.2) is defined as a solution to the following Poisson equation ∆φ= divm=(m0 ·∇)g(r) (3.2) The function φ can be sought in the form of (m0 ·∇)Φ(r), where Φ denotes a new unknown function. It is described by another Poisson’s equation. Indeed, insertion of the assumed form of φ in (3.2) gives 0=∆φ− (m0 ·∇)g=∆(m0 ·∇)Φ− (m0 ·∇)g=(m0 ·∇)(∆Φ−g) Since m0 is non-zero we conclude that: — 2D case, angular symmetry ∆Φ= g(r) ⇒ 1 r d dr ( r dΦ(r) dr ) = g(r) (3.3) — 3D case, spherical symmetry ∆Φ= g(r) ⇒ 1 r2 d dr ( r2 dΦ(r) dr ) = g(r) (3.4) Then, the derivative of Φ is given as dΦ dr =                1 r r ∫ 0 ξg(ξ) dξ for 2D case 1 r2 r ∫ 0 ξ2g(ξ) dξ for 3D case (3.5) On the magnetization-based Lagrangian methods... 345 The gradient of the function φ is expressed as follows ∇φ=∇[(m0 ·∇)Φ] =∇ (m0 ·r r dΦ dr ) (3.6) Finally, the velocity field induced by the single magneton with the center at the origin can be derived as u=                m0 [ g(r)− 1 r3 r ∫ 0 ξ2g(ξ)dξ ] − (m0 ·r)r r2 [ g(r)− 3 r3 r ∫ 0 ξ2g(ξ)dξ ] for 3D m0 [ g(r)− 1 r2 r ∫ 0 ξg(ξ)dξ ] − (m0 ·r)r r2 [ g(r)− 2 r2 r ∫ 0 ξg(ξ)dξ ] for 2D (3.7) In particular, the velocity induced at the center of the magneton is equal to u ∣ ∣ ∣ r=0 =        2 3 g(0)m0 for 3D case 1 2 g(0)m0 for 2D case (3.8) The reader may notice that, in the case of a finite core, the velocity induced outside the core is potential. Themagneton charge can be defined as Q=                a ∫ 0 ξ2g(ξ) dξ for 3D case a ∫ 0 ξg(ξ) dξ for 2D case (3.9) The symbol a denotes here the radius of the core. One can easily conclude from (3.7) that for r­ a the induced velocity is given by the formulas u=      −∇ ( Q m0 ·r r3 ) for 3D case −∇ ( Q m0 ·r r2 ) for 2D case (3.10) The boundary of the magneton core (the radius r = a) moves due to the self-induction with the velocity U =        2Q a3 m0 for 3D case Q a2 m0 for 2D case (3.11) 346 A.Styczek, J.Szumbarski The spherical shapeof the core is preservedduring the self-inducedmotion. Note that in the case when Q=0 there is no self-induction. In the case of an unbounded core, one must reasonably assume that the self-movement is equal to zero. This is equivalent to lim r→∞ 1 r3 r ∫ 0 ξ2g(ξ) dξ=0 for 3D case lim r→∞ 1 r2 r ∫ 0 ξg(ξ) dξ=0 for 2D case (3.12) Consider now the superposition of a number ofmagnetonswith characteri- stic vectors m0k andwith the centers located at points rk. The total velocity field is a sum of the contributions from each magneton v= ∑ k m0k(t)U(|r−rk|) (3.13) In the above U denotes the matrix operator defined by induction formulas (3.7). In the following sections of this part of the paper, we focuse on the 3D case. The two dimensional variant of the method can be derived analogously. 4. Lagrangian decomposition In the rest of this paper, we will focus on those variants of themagnetiza- tion equation, which can be cast in the following form ∂tm+(v ·∇)m= ν∆m+S(m) (4.1) Equation (4.1) can be re-written using thematerial derivative Dm Dt = ν∆m+S(m) (4.2) The above formdescribes the rate of change in time of themagnetization field seen by an observer moving with a fluid element. This form of the governing equation is an appropriate starting point for the formulation of a Lagrangian method. On the magnetization-based Lagrangian methods... 347 Amongdifferent variants of the governing equations discussed in Section 2, only first two cases lead to form (4.2). They differ with the shape of the stretching terms, namely S(m)=    (m ·∇)∇φ in case 1 (m ·∇)v in case 2 (4.3) In the remaining part of this Section, we explain the basic elements of the proposed numerical method. Wewill use theLagrangian decomposition of themagnetization field.First, we write the magnetization as the sum of magnetons currently present in the flow domain m= ∑ k m0kg(|r−rk|) (4.4) Themagnetonsmove accordingly to thefluidvelocity and self-inducedvelocity expressed by (3.11). The viscous effects can be simulated by randomwalks. Indeed, taking the set of the stochastic differential equations dxik = v i k dt+ √ 2ν dW ik x i k ∣ ∣ ∣ t=0 =xi0k (4.5) where dW ik denotes infinitesimal increments of the Wiener processes, we ob- tain the Fokker-Planck-Kolmogorov equation for the density of the transition probability (Gardiner, 1990) ∂tp+∂xi(v ip)= ν∆p (4.6) The initial condition formulated for the transition probability is following lim t→t0 p(t,r|t0,r0)= δ(r−r0) The solution to thehomogeneous advection-diffusionproblemcanbeexpressed as follows (Szumbarski and Styczek, 1997) m= ∫ Ω p(t,r|t0,r0)m0(r0) d3r0 ≈ ∑ k p(t,r|t0,r0k)4πm0k (4.7) Following Chorin and Marsden (1997), we introduce the following fractional- step approach. The time advancing the magnetization field at each interval [tn, tn+1] is carried out in three consecutive sub-steps: 348 A.Styczek, J.Szumbarski Step 1. Creation of new magnetons at the boundary At the beginning of a new time step the boundary distribution of the magnetization should be modified in order to satisfy the assumed boundary condition for the velocity field. Thus, we can write m ′ =m(n)+m (n+1) ∂Ω (4.8) where m (n+1) ∂Ω denotes the new magnetization created in close vicinity of the boundary at the beginning of the (n+1)th time step. This newmagnetization is such that the velocity ”induced” by the magnetization m′ defined as v′ = Pm′ satisfies the boundary condition. The magnetization field obtained in this sub-step can be expressed by introducing the operator B as follows m ′ =Bm(n) ≡m(n)+m(n+1)∂Ω (4.9) In the Lagrangian method, the magnetization m (n+1) ∂Ω is approximated by a certain number of new magnetons injected into the flow domain near the boundary of an immersed body. Detailed description of this procedure is given in the next Section. Step 2. Advection-diffusion In this sub-step, the magnetization is modified due to the transport by advection and diffusion. It means that all magnetons are moved to their new locations accordingly to the Itô equations dri′ =v ′ i′ dt+ √ 2ν dW i′ (4.10) The lower subscript i′ refers to all magnetons, including new magnetons in- jected near the boundary. Stochastic equations (4.10) are to be integrated over the time interval [tn, tn+1] using some numerical integration scheme. The length of the interval is τ = tn+1 − tn. Consequently, we obtain modified magnetization field. One can summarize this operation as follows m ′(r)= ∑ i′∈I′ m0i′g(|r−ri′|)⇒m′′(r)= ∑ i′∈I′′ m0i′g(|r−ri′ −∆ri′|) (4.11) In expression (4.11), the vector function mi′ describes the magnetization di- stribution inducedby the i′thmagneton. It should be remarked that the range I′′ of the second summation might be different than the range I′ of the first On the magnetization-based Lagrangian methods... 349 summation. This is so because a certain number of the magnetons always ”diffuses” across the boundary beyond the computational domain. The instantaneous magnetization field m′′ resulting from the advection- diffusion sub-step can be expressed by introducing another operator m′′ =Ξ1τm ′ (4.12) Step 3. Stretching In the last sub-step, one should account for the stretching term inmagneti- zation equation (4.1).Hence, the following (deterministic) differential equation ∂tm=S(m) (4.13) should be integrated over the time interval [tn, tn+1] with the initial condition m(tn)=m ′′. In effect, one obtains the final distribution of themagnetization at the time instant tn+1. Again, the operator notation can be used m(n+1) =Ξ2τm ′′ (4.14) The composition of the three operators introduced above defines completely the evolution of the magnetization field during one time step, namely m(n+1) =Ξ2τΞ 1 τBm (n) (4.15) The above formula can be used to obtain an approximate solution at the time instant t= t0+T (here t0 denotes the initial time). Indeed, fixing the number of time steps n, we can set τ =T/n and using recursively (4.15), we arrive at m ∣ ∣ ∣ t0+T ≈ (Ξ2τ Ξ1τB)nm ∣ ∣ ∣ t0 (4.16) The convergence of the splitting approach with is guaranteed by the Lie- Trotter formula (see Chorin and Marsden, 1997), which in our case can be written as follows m ∣ ∣ ∣ t0+T =(Ξ2T/nΞ 1 T/nB) nm ∣ ∣ ∣ t0 (4.17) being actually the limit form of (4.16). 350 A.Styczek, J.Szumbarski 5. Velocity decomposition and boundary conditions Consider a flow past a rigid body. The flow domain is unbounded and the stream at infinity is assumed homogeneous v ∣ ∣ ∣ ∞ =U∞ (5.1) On the surface of the body, the velocity field vanishes, i.e. v ∣ ∣ ∣ A =0 (5.2) We express the velocity field in the flow domain as a sum of the contributions induced by two ensembles of the magnetons. The first ensemble contains the magnetons created at earlier time steps and remaining in the flow domain at the considered time instant – these magnetons will be referred to as ”old” magnetons. The second ensemble is formedwith themagnetons created in the flow domain at the current instant of time, i.e. with ”new”magnetons. Thus, the complete velocity field can be written as follows v=U∞+ ∑ k m 0 0kU(|r−r0k|)+ M ∑ k=1 m n 0kU(|r−rnk|) (5.3) It is assumedthat all parameters m00k and r 0 k characterising thecurrent ”char- ge” and position of all old magnetons are known. On the other hand, the pa- rameters mn0k and r n k of the new magnetons are unknown and have to be determined at each time step of the flow simulation.We assume that the new magnetons will be generated near the body surface A. The surface A is divided into M small parts A1,A2, ..,AM, and the kth new magneton is placed over Ak, with its center located at a small distance above the surface.Thepositions rnk,k=1, ..,M of the ”generation points” are then defined and the only remaining unknowns are the characteristic vectors mn0k, k=1, ..,M. In can be concluded from (3.10) that m0U(r)→ 0 as r → ∞, and thus, asymptotic condition (5.1) is satisfied. Using boundary condition (5.2) one writes for r∈A M ∑ k=1 m n 0kU(|r−rnk|)=−U∞− ∑ k m 0 0kU(|r−r0k|) (5.4) Let eα ∣ ∣ ∣ Ai , α=1,2,3 denote the triple of the versors at a certain collocation point within the surface segment Ai. A possible choice is the triple consisting On the magnetization-based Lagrangian methods... 351 of the normal and two different tangent versors at each collocation point, or simply the three versors of the global Cartesian reference frame. Computing the scalar product of the last equation with three versors at each collocation point, we obtain (α=1,2,3) M ∑ k=1 m n 0kU(|rAi −r n k|)eα ∣ ∣ ∣ Ai =−U∞eα ∣ ∣ ∣ Ai − ∑ k m 0 0kU(|r−r0k|)eα ∣ ∣ ∣ Ai (5.5) where rAi denotes the position of a collocation point belonging to the seg- ment Ai. Consequently, we have a linear algebraic system of 3M equations. If the shape of the surface and other geometrical parameters are fixed in time, thematrix of the system is fixed as well and it can be evaluated (and possibly LU-factorized) once and forever. Solution to the linear system yields complete information about new gene- ration of the magnetons, and thus about the complete velocity and vorticity fields at a given instant of time. The velocity field satisfies the boundary con- dition at the chosen array of the collocation points. It should be remarked that alternative approaches to the boundary condi- tions are also possible. One can consider an integral-type rather than point- wise enforcement of the boundary condition on the surface. Such an approach would consist in integration of the tangent velocity along two different (possi- bly perpendicular) line segments located within the surface element Ai, and demanding these integrals to be equal to zero. The normal component of the velocity is set to zero at the collocation point, as previously. 6. Advection-diffusion and stretching The advection-diffusion operator working over one time step has been al- readymentioned. It describes a one-time-step advancing of themagnetization field. Let us summarize the essential steps defining this operator. Having com- plete knowledge of the magnetons created up to a given instant, we construct a ”magnetization layer”, which is adjacent to the boundary and consists of the newly born magnetons. These new magnetons are created at fixed loca- tions but their vector ”charges” are initially unknown. They are chosen so that the boundary condition for the velocity field is satisfied. At this point one has to solve the linear system of 3M equation, where M denotes the number of the newly generatedmagnetons. After that, themagnetization field is completely defined and all magnetons are moved to their new locations by 352 A.Styczek, J.Szumbarski performing one-time-step integration of stochastic differential equations (4.2). The instantaneous velocity of the jth magneton is expressed as vj =U∞+ 2Q a3j m0j + ∑ k 6=j m0kU(|rj −rk|) (6.1) i.e. as the sum of the free-stream velocity, self-induced velocity (3.11) and the velocity induced by other magnetons in the flow domain. The change of the positional vector rj of the magneton during one time step is determined by numerical integration of Itô differential equations (4.10). As an example, one can consider the Euler integration scheme leading to the following formula ∆rj =vj∆t+ √ 2ντ N(0,1)j (6.2) In the above, τ denotes the length of the time step and N(0,1)j denotes the random vector, the components of which are independent random variables with the standardGaussian distribution.More sophisticated, higher-order in- tegration schemes are also available. The reader should refer to Kloeden and Platen (1999) for more details. It should be noted that due to the presence of the random component in equation (6.2), the magnetons can penetrate the surface and jump into the body interior. This is a manifestation (on the ”kinetic” level) of the dif- fusion of the magnetization through the boundaries. The magnetons beyond the flow domain can be removed from the simulation or they can be reflected back.Kinetic equations (6.2), togetherwith the removal or reflectionapproach, complete the description of the advection-diffusion operator Ξ1τ. The last operator is defined by a stretching equation.We re-write it in the following form ∂tm α =−mβ∂xαvβ (6.3) in case 2, or as ∂tm α =mβ∂xα,xβφ (6.4) in case 1. Bothmatrix operators appearing in the right-hand sides of (6.3) and (6.4) are known. We substitute here the Lagrangian decomposition and solve equ- ation (6.3) or (6.4) for each magneton. This approach is justified by the follo- wing fact. If for a given matrix operator Bαβ and vectors mk, k=1,2, ..,M one has ∂tm α k =Bαβm β k (6.5) On the magnetization-based Lagrangian methods... 353 then, as a result of the formal linearity with respect to m, the basic equation holds. Inserting the expression for mαk we write g(|r−rk|) dmα0k dt =    −mβ0kg(|r−rk|)∂xαvβ m β 0kg(|r−rk|)∂xα,xβφ (6.6) Integrating the above equation over the support of the core function g we get 4π dmα0k dt =B0αβm β 0k (6.7) where the matrix operator depends on rk only and has the following form B0αβ = ∫ ρ¬a { vβ(rk+ρ) −∂xβφ(rk+ρ) } dg(ρ) dρ xα ρ dΩρ (6.8) In the derivation of (6.8) the equality g(a)=0 and the integration by parts have been employed. As we see, the stretching is described by the system of 3N ordinary dif- ferential equations with constant coefficients. Its dimension is equal to the tripled number of all magnetons present in the flow domain at a given instant of time. Thus, its dimension is large and can be excessively large for long flow simulations. 7. Final remarks The Lagrangian approach to the velocity-vorticity formulation known as the Vortex Blobs Method (see Protas (2000) for exhaustive list of references) allows solvingmany interesting problems. In the vortex method, the formula- tion involves only kinematic quantities: the velocity and vorticity. Once these fields are determined, the pressure can be recovered (protas et al., 2000). As a physical quantity, the pressuremust be a univalued function of time and spa- ce variables. This fact brings a fundamental condition constraining the total charge of the vorticity, being nontrivial for flows in multi-connected domains. It is essential that the implementation of this condition gives an effect of stabi- lization of a numerical process (see Błażewicz and Styczek, 1993; Szumbarski andStyczek, 1999; Protas, 2000), which otherwise tends to generate physically absurd ”solutions”. 354 A.Styczek, J.Szumbarski The existence of an additional condition of global nature is one of themain differences between the vortex and the magnetization methods. The pressure recovery does not require any further restriction of themagnetization because it is based on explicit formula (2.5) (in the vortex method it is based on a formula for the pressure gradient). There exists, however, a numerical evidence, which seems to indicate that some kind of (yet unknown) stabilizing condition is indeed necessary. Russo and Smerka (1999) and Summers and Chorin (1996) presented a sample of numerical computations. They encountered unstable behavior rendering im- possible sensible long-time simulations of flow evolution. These observations are in agreement with the results obtained by the authors, discussed in Sec- tion 2 of this paper. Acknowledgements This work has been supported by the State Committee for Scientific Research (KBN), grant No 7T07A02214. 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: No 1, Part 2: No 4 2. Chorin A.J., 1994,Vorticity and Turbulence, Springer Verlag 49, 1, 209-232 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. Gardiner C.W., 1990, Handbook of Stochastic Methods, 3rd Ed. Springer Verlag 5. Kloeden P.E., Platen E., 1999,Numerical Solution of Stochastic Differen- tial Equations, Springer Verlag, NewYork 6. Protas B., 2000, Analysis and control of aerodynamic forces in the plane flow past a moving obstacle – application of the vortexmethod, Ph.D. Thesis, WarsawUniversity of Technology 7. Protas B., Styczek A., Nowakowski A., 2000, An effective approach to computation of forces in viscous incompressible flow, J. Comp. Phys., 159, 231-245 8. Russo G., Smerka P., 1999, Impulse formulation of the Euler equations: general properties and numerical methods, Journal of Fluid Mechanics, 391, 189-209 On the magnetization-based Lagrangian methods... 355 9. Summers D.M., Chorin A.J., 1996, Numerical vorticity creation based on impulse conservation,Proc. Nath. Acad. Sci. USA, 93, 1881-1885 10. Szumbarski J., Styczek A., 1997, The stochastic vortexmethod for viscous incompressible flow in a spatially periodic domain,Arch. Mech., Lagrangeowska metoda magnetyzacji dla dwu i trójwymiarowych ruchów płynu lepkiego. Część I – Podstawowe sformułowania Streszczenie W artykule przedstawione jest sformułowanie problemu granicznego dla równań Naviera-Stokesazużyciemtzw.polamagnetyzacji. Sformułowanienie jest jednoznacz- ne, lecz wiąże się z przyjętą transformacją cechowania. Rozważane są różne postacie tej transformacji i dokonuje się wyboru odpowiednich wariantów. Pole magnetyzacji przedstawione jest w formie lagrangeowskiej.Wprowadza się cząstki będące źródłami tego pola i określa się związane z ich zbiorem pole prędkości. Cząstki magnetyzacji (zwane magnetonami) poruszają się w indukowanym polu prędkości, wykonują ruch losowy odpowiadający dyfuzji i podlegają przekształceniuw sposób opisany członem źródłowym (tzw. stretching). Warunek brzegowy formułowany na opływanym ciele jest realizowany przez tworzenie w każdej chwili nowych cząstek ulokowanych w bli- skim otoczeniu powierzchni ciała. Manuscript received March 2, 2001; accepted for print August 31, 2001