The aim of this note is to present a numerical method to solve the Stokes
problem in a bounded domain with a Dirac source term, which preserves
optimality for any approximation order by the finite-element method. It is
based on the knowledge of a fundamental solution to the associated operator
over the whole space. This method is motivated by the modeling of the movement
of active thin structures in a viscous fluid.Comment: Comptes Rendus M{\'e}canique, Elsevier, 2015,
  http://www.sciencedirect.com/science/article/pii/S163107211400181

Lacouture, Loïc

English

arXiv

Loïc Lacouture

Comptes Rendus Mécanique (CRM)

C. R. Mecanique 343 (2015) 187–191Contents lists available at ScienceDirectComptes Rendus Mecaniquewww.sciencedirect.comA numerical method to solve the Stokes problem with a punctual force in source termUne méthode numérique pour la résolution du problème de Stokes avec une force ponctuelle en terme sourceLoïc LacoutureUniversité Paris-Sud, Laboratoire de Mathématiques d’Orsay, CNRS UMR 8628, Faculté des sciences d’Orsay, bâtiment 425, 91405 Orsay cedex, Francea r t i c l e i n f o a b s t r a c tArticle history:Received 19 June 2014Accepted 17 September 2014Available online 12 November 2014Keywords:Error estimatesFinite element methodStokesletThin structuresMots-clés :Estimations d’erreurMéthode éléments finisStokesletStructures finesThe aim of this note is to present a numerical method to solve the Stokes problem in a bounded domain with a Dirac source term, which preserves optimality for any approximation order by the finite-element method. It is based on the knowledge of a fundamental solution to the associated operator over the whole space. This method is motivated by the modeling of the movement of active thin structures in a viscous fluid.© 2014 Académie des sciences. Published by Elsevier Masson SAS. All rights reserved.r é s u m éLe but de cette note est de présenter une méthode numérique pour la résolution du problème de Stokes avec une force ponctuelle en terme source, qui assure l’optimalité de l’erreur d’approximation éléments finis. Elle s’appuie sur la connaissance explicite d’une solution fondamentale de l’opérateur linéaire associé. Cette méthode est motivée par la modélisation du mouvement de structures fines actives dans un fluide visqueux.© 2014 Académie des sciences. Published by Elsevier Masson SAS. All rights reserved.Version française abrégéeL’étude du mouvement de structures fines actives dans un fluide visqueux, tels que les flagelles permettant la nage de bactéries ou les cils impliqués dans le transport mucociliaire, conduit à considérer le problème de Stokes avec un second membre singulier. Dans l’asymptotique d’un cil dont le diamètre tend vers 0 et la vitesse vers l’infini, le terme source est en fait une distribution linéique de forces. Dans le but de pouvoir faire des calculs, puisque intégrer numériquement le long d’une courbe quelconque est difficile, nous approchons la distribution linéique de forces δΓ par une somme de forces ponctuelles ∑ciδi . Une preuve basée sur celle du théorème des sommes de Riemann permet de montrer qu’il y a convergence, au sens faible dans H−3/2−s , pour tout s > 0, de ∑ciδi vers δΓ lorsque le nombre N de masses de Dirac tend vers l’infini. On peut aussi préciser la convergence dans des espaces plus faibles, voir (1). La convergence des solutions E-mail address: loic.lacouture@math.u-psud.fr.http://dx.doi.org/10.1016/j.crme.2014.09.0081631-0721/© 2014 Académie des sciences. Published by Elsevier Masson SAS. All rights reserved.188 L. Lacouture / C. R. Mecanique 343 (2015) 187–191associées se déduit de l’inégalité (2), tirée de [1]. On est alors ramené à l’étude du problème de Stokes avec une force ponctuelle en terme source.Lorsqu’on considère un problème elliptique avec une masse de Dirac en second membre, en dimension d  2, ce second membre n’étant pas dans H−1, le problème sort du cadre variationnel standard basé sur l’espace de Sobolev H1. Si la méthode des éléments finis peut être définie au niveau discret, les résultats de convergence classiques ne sont a priori plus valables. Dans le cas du problème de Poisson, qui peut être vu comme une version scalaire et simplifiée du problème de Stokes, Scott a démontré dans [2] que la méthode éléments finis P1 converge en norme L2 à l’ordre 1 en 2d et 1/2 en 3d. Des estimations similaires ont été obtenues dans [3] avec une méthode de Galerkin discrète. De plus, Apel et ses co-auteurs ont montré dans [4] qu’en raffinant le maillage autour de la singularité, on retrouvait l’ordre de convergence classique. La méthode présentée, basée sur la connaissance explicite d’une solution fondamentale de l’opérateur linéaire associé, fait partie d’une classe de méthodes dites de soustraction, introduites en électroencéphalographie [5]. Elle permet de retrouver les ordres de convergence classiques sans raffinement de maillage.Pour fixer les idées, nous allons nous intéresser au problème de Stokes avec des conditions aux limites de type Dirichlet homogènes, voir le problème (4). La particularité de ce problème réside dans la singularité du second membre : un Dirac de force appliqué en un point x0 du domaine Ω . Pour cet opérateur, on connaît une solution fondamentale définie en domaine infini, appelée Stokeslet, que l’on note (uδ, pδ), voir (5). On obtient la solution (u, p) du problème (4) en ajoutant à (uδ, pδ) un relèvement régulier prenant ainsi en compte les conditions aux bords. La singularité de la solution (u, p) est contenue dans la solution fondamentale (uδ, pδ), et elle est localisée au point x0. Le principe de la méthode qui suit consiste à capturer cette singularité pour se ramener à la résolution d’un problème auxiliaire régulier.On commence par définir une fonction plateau χ , régulière, valant 1 sur un voisinage de x0 et 0 loin de ce point, voir Définition 3.1. On note ensuite u0 = χuδ et p0 = χ pδ , et g et h les fonctions définies en (6). D’après ces définitions, on remarque que les supports de g et h sont contenus dans une couronne centrée en x0, voir Fig. 1. De plus, les fonctions uδet pδ étant analytiques en dehors de x0, la régularité des fonctions g et h dépend directement de celle de la fonction χ . Finalement, pour obtenir la solution (u, p) de (4), il suffit de corriger les termes d’erreur g et h introduits en (6) en résolvant le problème elliptique régulier (7), dont on note v la solution. En effet, la fonction u := u0 +v est la solution du problème (4).Cette méthode permet de passer de la résolution d’un problème singulier à celle d’un problème auxiliaire régulier. Alors que le premier converge à un ordre faible [2], le second converge à l’ordre optimal, quel que soit l’ordre des éléments utilisés. En notant uh := u0 + vh , où vh est la solution numérique du problème (7) obtenue par une méthode éléments finis, on déduit de (8) que l’erreur commise sur u est la même que celle commise sur v, et on montre ainsi que la vitesse de convergence est optimale.Par exemple, si on utilise une méthode éléments finis P1, k = 0 suffit. On définit alors χ comme en (9), et on explicite g et h, valant respectivement (10) et (11) en dimension 2, et (12) et (13) en dimension 3. Après résolution numérique du problème (7), on obtient finalement une solution approchée uh dont l’erreur ‖u − uh‖L2 est en O (h2), quelle que soit la dimension, contre une erreur, avec une méthode directe, en O (h) en dimension 2 et en O (√h) en dimension 3.Cette méthode, présentée dans le cas du problème de Stokes, peut se généraliser à d’autres problèmes elliptiques li-néaires, comme le problème de Poisson avec une masse de Dirac en second membre. Les conditions aux limites de type Dirichlet homogènes peuvent aussi être remplacées par des conditions de type Dirichlet non homogènes, Neumann ou Ro-bin. Enfin, la linéarité, qui joue un rôle essentiel, permet en outre de résoudre le cas où le second membre est la somme d’un nombre fini de forces ponctuelles et d’une fonction lisse, tout en ne résolvant qu’un seul problème numérique.1. IntroductionIn order to model active thin structures in a viscous fluid, such as flagella connected to bacteria or cilia involved in the mucociliary transport, we have studied the Stokes problem with a singular right-hand side. In the asymptotic of a zero diameter cilia with an infinite velocity, the source term is a lineic distribution of forces, which, in order to ease computations, will be approximated by a sum of punctual forces. After having justified this approximation, we will present a numerical method to solve the Stokes problem with a punctual force in source term, and illustrate the results by numerical simulations.2. Approximation of the lineic distribution of forces by a sum of punctual forcesSince calculating an integral on any curve is numerically very difficult, the source term, noted δΓ , the lineic distribution of forces on a curve Γ , is approached by a sum of N punctual forces ∑ciδi uniformly distributed along Γ . The theorem of Riemann sums ensures that ∑ciδi weakly converges to δΓ in H−3/2−s , for all s > 0. Working in weaker spaces, it is possible to adapt the proof of the theorem of Riemann sums and specify the convergence:∥∥∥∥∥δΓ −N∑i=1ciδi∥∥∥∥∥H−2−s C√Nand∥∥∥∥∥δΓ −N∑i=1ciδi∥∥∥∥∥H−5/2−s CN(1)Moreover, using a result proved by Lions and Magenes in [1], which can be written in this caseL. Lacouture / C. R. Mecanique 343 (2015) 187–191 189‖u‖H2−r  C‖ f ‖H−r , ∀r  0 (2)where u is the solution to a regular elliptic problem with a source term f ∈ H−r , we can conclude that the solution uN of the Stokes problem with ∑ciδi right-hand side converges to the solution uΓ to the Stokes problem with δΓ source term, when N goes to infinity. Actually, we have:‖uΓ − uN‖−s  C√Nand ‖uΓ − uN‖−1/2−s  CN(3)Finally, the solution to the Stokes problem with a lineic distribution of forces is approached by the solution to the Stokes problem with a finite sum of punctual forces in the source term. By linearity and without loss of generality, in the following we will deal with a single punctual force.3. Numerical method to solve the Stokes problem with a Dirac source termIn dimension d  2, the δ-distribution is not continuous on H1, and so the solution to an elliptic problem with Dirac source term is not regular. Consequently, classical results for the convergence of the finite element method are not valid. In the case of the Poisson problem, which can be seen as the scalar version of the Stokes problem, Scott has shown in [2] that the P1-finite element method converges for L2-norm at the order 1 in dimension 2 and at the order 1/2 in dimension 3. Similar estimates have been obtained in [3] with a discrete Galerkin method. Moreover, it has been shown by Apel and his co-authors [4] that using graded meshes, it is possible to get numerically the classical order of convergence. The aim of this section is to present a numerical method that preserves optimality for any approximation order, without using mesh grading. It is based on the knowledge of a fundamental solution to the considered linear elliptic problem. This approach fits on the class of subtraction methods, introduced in [5] in the context of electroencephalography.3.1. Principle of the methodLet us consider the following problem, defined on a bounded open domain Ω ⊂ Rd ,⎧⎨⎩−μu + ∇p = δx0 F in Ωdiv u = 0 in Ωu = 0 on ∂Ω(4)where x0 is fixed in Ω and F is a vector of Rd . Let us note that a fundamental solution to problem (4) is known in dimensions 2 and 3:• d = 2, uδ(x) = 14πμ(− ln(|x|)I2 + xtx|x|2)F and pδ(x) = 12πx · F|x|2• d = 3, uδ(x) = 18πμ(I3|x| +xtx|x|3)F, and pδ(x) = 14πx · F|x|3 (5)where Id is the identity matrix. The fundamental solution (uδ, pδ) does not satisfy the boundary conditions, and so it is not the solution to problem (4). But this solution can be retrieved by adding a regular lifting, therefore the whole information on the singularity of the solution (u, p) is contained in the fundamental solution (uδ, pδ) and is located at x0. In order to extract this singularity, let us fix 0 < a < b < d(x0, ∂Ω) and define χ by Definition 3.1.Definition 3.1. Assume that χ is a bump function satisfying for some k  0,• χ ∈ H2+k(Rd),• χ|B(x0,a) = 1,• χ|B(x0,b)c = 0.Then, with u0 := χuδ and p0 := χ pδ , we define g and h asg = −μu0 + ∇p0 − δx0 F and h = div u0 (6)By the definitions of uδ , pδ and χ , supp(g) ⊂ Rba(x0) and supp(h) ⊂ Rba(x0), where Rba(x0) is the ring centered around x0, of internal radius a and external radius b, see Fig. 1. Since uδ and pδ are analytic on Ω \ {x0}, the regularity of functions g and h directly depends on the regularity of function χ , namely g ∈ Hk(Ω) and h ∈ Hk+1(Ω). Finally, it only remains to correct the terms u0 and p0 by solving the regular elliptic problem⎧⎨⎩−μv + ∇q = −g in Ωdiv v = −h in Ωv = 0 on ∂Ω(7)190 L. Lacouture / C. R. Mecanique 343 (2015) 187–191Fig. 1. Domain Ω .Fig. 1. Domaine Ω .and the solution to problem (4) is given by (u, p) := (u0 + v, p0 + q) = (χuδ + v, χ pδ + q), where u0 and p0 are explicitly known functions and (v, q) is the solution to problem (7). Noting (vh, qh) the numerical solution to problem (7) and defining uh := vh + u0 and ph = qh + p0, we have:‖u − uh‖Hs(Ω) = ‖v − vh‖Hs(Ω), for 0  s  k + 1‖p − ph‖Hs(Ω) = ‖q − qh‖Hs(Ω), for 0  s  k (8)Actually, this method allows us to switch from the numerical computation of the solution of a singular problem with Dirac source term (with a poor convergence rate) to the numerical computation of the solution of a regular problem with an optimal convergence rate, at any required precision in terms of regularity.3.2. Practical aspectsFor the sake of simplicity, the location of the Dirac source term will be the origin. First, we need to choose a suitable function χ . Actually, to take advantage of using P	-finite elements, 	  1, χ has to be H	+1(Rd), in order to ensure that g ∈ H	−1(Ω) and h ∈ H	(Ω), and finally to get an optimal order of convergence. For instance, for 	 = 1, let us define χ , as a radial function, by:χ(r) =⎧⎪⎨⎪⎩1 for r ∈ [0,a]2r3−3(a+b)r2+6abr+b2(b−3a)(b−a)3 for r ∈ [a,b]0 for r > b(9)where the function r is defined on Rd byr(x) = ‖x‖2The function g and h can be made explicit. According to this definition of χ , g and h vanish outside the ring a < ‖x‖2 < b. For a < ‖x‖2 < b, the expressions of g and h depend on the dimension,• for d = 2,g(x) = 32π(b − a)3r[((3r2 − 2(a + b)r + ab) ln r + 2r2 − 2(a + b)r + 2ab)I2 + (ab − r2)xtxr2]F (10)h(x) = 3(1 − ln r)(r2 − (a + b)r + ab)2πμ(b − a)3r x · F (11)• for d = 3,g(x) = 34π(b − a)3r2[((a + b)r − 2r2)I3 + (2ab − (a + b)r)xtxr2]F (12)h(x) = 3(r2 − (a + b)r + ab)2πμ(b − a)3r2 x · F (13)4. Numerical illustrationsIn this section, we illustrate our theoretical results by a numerical example. We define Ω as the unit square and x0 =(0.5, 0.5). Table 1 presents the L2-error for a direct method and a subtraction method respectively, for a characteristic mesh L. Lacouture / C. R. Mecanique 343 (2015) 187–191 191Table 1L2-error and estimated order of convergence for the direct method and the substraction method.Tableau 1Erreur L2 et ordre de convergence pour la méthode directe et la méthode de soustraction.h 2−3 2−4 2−5 2−6 2−7 e.o.c.Direct method 1.02 × 10−2 4.87 × 10−3 2.36 × 10−3 1.21 × 10−3 5.89 × 10−4 1.02Subtraction method 4.12 × 10−3 1.33 × 10−3 2.92 × 10−4 6.86 × 10−5 2.71 × 10−5 1.88Fig. 2. Section {y = 0.5} of the error |u − uh | for the direct method and the subtraction method with h = 0.125.Fig. 2. Coupe {y = 0.5} de l’erreur |u − uh | pour la méthode directe et la méthode de soustraction avec h = 0.125.size h, and the estimated order of convergence (e.o.c.). Fig. 2 illustrates the section {y = 0.5} of the error |u −uh| in the both cases. Numerical simulations evidence the fact that solving the auxiliary problem associated with the subtraction procedure of the singularity is more efficient than solving directly the problem with the Dirac source term.5. ConclusionTo model active thin structures in a viscous fluid, such as flagella connected to bacteria or cilia involved in the mu-cociliary transport, we have studied Stokes’ problem with a singular right-hand side: a punctual force. However, when this problem is solved numerically, the singularity causes a poor convergence of the approximate solution to the exact solution. The method presented in this note preserves optimality for any approximation order, without using mesh grading. If the examples are treated with homogeneous Dirichlet conditions, the same method is still valid in the case of non-homogenous Dirichlet or any affine boundary conditions (Neumann, Robin...), up to suitable adaptations. Similarly, the method can be generalized to the problem with a sum of a finite number of Dirac masses and a smooth term right-hand side.References[1] J.-L. Lions, E. Magenes, Problèmes aux limites non homogènes et applications, vol. 1, Dunod, Paris, 1968.[2] R. Scott, Finite element convergence for singular data, Numer. Math. 21 (1973) 317–327.[3] P. Houston, T.P. Wihler, Discontinuous Galerkin methods for problems with Dirac delta source, ESAIM Math. Model. Numer. Anal. 46 (2012) 1467–1483.[4] T. Apel, O. Benedix, D. Sirch, B. Vexler, A priori mesh grading for an elliptic problem with Dirac right-hand side, SIAM J. Numer. Anal. 49 (2011) 992–1005.[5] C.H. Wolters, H. Köstler, C. Möller, J. Härdtlein, L. Grasedyck, W. Hackbusch, Numerical mathematics of the subtraction method for the modeling of a current dipole in EEG source reconstruction using finite element head models, SIAM J. Sci. Comput. 30 (2007) 24–45.

A numerical method to solve the Stokes problem with a punctual force in source term

Abstract

Similar works

Full text

Available Versions

Comptes Rendus Mécanique (CRM)