Calcul numérique de l'intégrale de surface

Pour calculer la seconde intégrale dans (6),

$$\mathcal{I}_c(x) = \int\!\!\!\int_{\Sigma_c} \ensuremath{\m... ...Sigma_c} \Phi(x,y) \ \ensuremath{\mathbf{n}} (y) \mbox{ d}s(y)$$ (9)

un calcul analytique n'est pas envisageable.

La méthode de calcul de cette intégrale consiste à approcher la surface $\Sigma_c$ par une surface $\Sigma_h$ obtenue par interpolation quadratique par morceaux à partir d'un maillage de la surface $\Sigma_c$. On décompose l'intégrale sur $\Sigma_h$ en une somme d'intégrales sur les triangles $\widetilde{K}$ formant la surface $\Sigma_h$. La normale $\ensuremath{\mathbf{n}}$ est approchée par la normale au triangle $\widetilde{K}$et l'intégrale est évaluée à l'aide d'une formule de quadrature numérique.

Triangulation de la surface

On suppose que $\Sigma_c$ est l'image d'un domaine D polygonal fermé de ${\Bbb R}^2$ par un ${\cal C}^\infty$-difféomorphisme F, $\Sigma = F_j(D)$.

Pour chaque domaine D, on introduit une triangulation $\mit \widehat{T}_h$ de D au sens des triangulations usuelles en éléments finis,

$$D = \bigcup_{\widehat{K}\in\mit\widehat{T}_h} \widehat{K} .$$

Remarque: s'il n'existe pas de ${\cal C}^\infty$-difféomorphisme F tel que $\Sigma = F(D)$, il est possible de décomposer la surface $\Sigma_c$ en J morceaux $\Sigma_j ,\ j=1,\cdots, J$ pour lesquels $\Sigma_j$ est l'image d'un domaine D_jpolygonal fermé de ${\Bbb R}^2$ par un ${\cal C}^\infty$-difféomorphisme F_j.

On désigne par $\widehat K$ l'élément générique de la triangulation $\mit \widehat{T}_h$. Ses sommets sont notés $\widehat{S}_1,\widehat{S}_2$ et $\widehat{S}_3$. Sans restreindre la généralité, on supposera que D est un triangle rectangle isocèle.

On obtient une triangulation ${\mit T_h}$ de la surface $\Sigma_c$ en prenant l'image par F de la triangulation $\widehat{\mit{T}}_h$ de D. L'image du triangle $\widehat K$ par l'application F est noté K. On note S₁, S₂ et S₃ ses sommets qui sont définis par $S_i = F(\widehat{S}_i) , \ \ i=1,2,3$.

Soit $\sigma$ le triangle de référence

$$\sigma = \{(s,t)\in [0,1]^2 \ \backslash\ 0 \leq s+t \leq 1 \} .$$

On désigne par $\{\rho_i\}_{i=1...6}$ les trois sommets et les trois milieux d'arêtes de $\sigma$. On considère $f_{\widehat K}$ l'application linéaire de ${\Bbb R}^2$ dans ${\Bbb R}^2$ définie par les relations

$$f_{\widehat K} (\rho_i) = \widehat{S}_i,\hspace*{1cm} i=1,2,3 .$$

Cette application associe à tout point $(s,t) \in \sigma$, le point $(\hat{x}^1,\hat{x}^2) \in \widehat{K}$ tel que

$$(\hat{x}^1,\hat{x}^2) = (1-s-t)\ \widehat{S}_1 + t\ \widehat{S}_2 + s\ \widehat{S}_3 .$$ (10)

On note m_K l'application de ${\Bbb R}^2$ dans ${\Bbb R}^3$définie par

$$m_K = F \circ f_{\widehat K} .$$ (11)

À tout point $(s,t) \in \sigma$ l'application m_K associe le point $(x^1,x^2,x^3) \in K$ tel que

$$(x^1,x^2,x^3) = m_K (s,t) = F(f_{\widehat K}(s,t)) = F((1-s-t)\ \widehat{S}_1 + t\ \widehat{S}_2 + s\ \widehat{S}_3) .$$ (12)

Approximation de la surface

Pour approcher la surface $\Sigma_c$ on remplace l'application F, pour chaque triangle $\widehat{K}\in \widehat{T}_h$, par $\widetilde{F}_{\widehat K}$l'interpolant de Lagrange de degré 2 sur $\widehat K$. L'interpolant $\widetilde{F}_{\widehat K}$ vérifie

$$\widetilde{F}_{\widehat K}(\widehat{S}_i) = F(\widehat{S}_i) = S_i \hspace*{1cm} i=1,\cdots, 6 .$$

On note L_i, $i=1,\cdots, 6$ les fonctions de base de l'interpolation de Lagrange de degré 2 sur $\sigma$. On note encore $\widetilde{m}_K$ l'interpolant d'ordre 2 de m_K donné par

$$\widetilde{m}_K(s,t) = \sum_{i=1}^{6} S_i\ L_i (s,t) = \sum_{i=1}^{6} m_K(\rho_i) \ L_i (s,t) .$$ (13)

La surface $\Sigma_c$ est approchée par la surface $\Sigma_h$ définie par

$$\Sigma_h = \bigcup_{\widehat{K} \in {\mit \widehat{T}_h}} \w... ...ehat K) = \bigcup_{\widetilde{K}\in{\cal T}_h}\widetilde{K} .$$ (14)

L'élément générique $\widetilde{K}$ de la triangulation ${\cal T}_h$est l'image du triangle $\widehat{K}\in{\mit \widehat{T}_h}$ par l'application  $\widetilde{F}_{\widehat K}$.

Approximation de l'intégrale

Nous avons $\mathcal{I}_c$ donnée par (10) qui s'écrit

$$\mathcal{I}_c(x) \sum_{\widehat{K}\in{\mit \widehat{T}_h}} \int\!\!\!\int_{F(\widehat{K})} \Phi(x,y)\ \ensuremath{\mathbf{n}} (y) \mbox{ d}s(y)$$ =

$$\sum_{\widehat{K}\in{\mit \widehat{T}_h}} \int\!\!\!\int_\sigma \Phi(x,m_K(s,t))\ \ensuremath{\mathbf{n}} (m_K(s,t))\ J(s,t) \mbox{ d}s \mbox{ d}t$$ = (15)

où

$$J(s,t) =\ \mid \partial_s m_K(s,t) \wedge \partial_t m_K(s,t) \mid .$$

Le vecteur normal unitaire à la surface $\Sigma_c$ est donné par

$$\ensuremath{\mathbf{n}} (m_k(s,t)) = n(s,t) = \frac{\partial_... ...} {\mid \partial_s m_K(s,t) \wedge \partial_t m_K(s,t) \mid} .$$ (16)

Il est en général mal aisé pour une surface $\Sigma_c$ quelconque d'expliciter l'application m_K et de calculer J. Pour pouvoir calculer $\mathcal{I}_c$, nous avons donc recours à l'interpolant $\widetilde{m}_K$ d'ordre 2 de m_K donné en (14). Ceci nous amène à considérer l'approximation suivante de $\mathcal{I}_c$,

 

$$\mathcal{I}_h(x) \sum_{\widetilde{K}\in{\cal T}_h} \int\!\!\!\int_\sigma \Phi(x,\w... ...ilde{m}_K(s,t)) \ \widetilde{n}(s,t)\ \widetilde{J} (s,t) \mbox{ d}s \mbox{ d}t$$ =

$$\sum_{\widetilde{K}\in{\cal T}_h} \int\!\!\!\int_{\widetilde K} \Phi(x,y) \ \widetilde{n}(y)\ \mbox{ d}s(y)$$ = (17)

$$\int\!\!\!\int_{\Sigma_h} \Phi(x,y) \ \widetilde{n}(y)\ \mbox{ d}s(y)$$ =

où

$$\tilde J(s,t) = \mid \partial_s \widetilde{m}_K(s,t) \wedge \partial_t \widetilde{m}_K(s,t) \mid$$ (18)

et ${\widetilde n}(s,t)$ est le vecteur

$${\widetilde n}(s,t) = \frac{\partial_s \widetilde{m}_K(s,t) \... ...ilde{m}_K(s,t) \wedge \partial_t \widetilde{m}_K(s,t) \mid} .$$ (19)

Il n'est pas possible de calculer l'intégrale dans (18) de manière exacte. Nous utilisons une méthode d'intégration numérique. Nous avons mis en évidence que des difficultés numériques liées à l'utilisation d'une formule de quadrature numérique apparaissaient quand le point x se situait à proximité de la surface $\Sigma_c$. Bien que les intégrales ne soient pas à proprement parler singulières, la précision de l'évaluation numérique se dégrade sérieusement en raison de la grande rapidité des variations de l'intégrant au voisinage de la surface (il se comporte en fonction de la distance à la surface r, comme r^-2).

Ce type d'intégrales quasi-singulières se retrouve lors de la mise en uvre des méthodes d'équation intégrale de frontière. Deux approches sont possibles afin de calculer avec précision les intégrales singulières [#!Huang93!#]. La première consiste à utiliser des transformations géométriques (changements de variables) afin d'affaiblir la singularité. La seconde consiste à augmenter le nombre de nuds de quadrature, soit en augmentant l'ordre de la formule de quadrature, soit en divisant le domaine d'intégration. Nous avons adopté cette dernière solution, en raison de la simplicité de sa mise en uvre et de son efficacité. Nous reprenons une méthode introduite par K.E. Atkinson dans [#!Atkinson85b!#], [#!Atkinson89!#]. Cette méthode consiste à distinguer deux types de formules de quadrature en fonction de la distance du point P au triangle $\widetilde{K}$ sur lequel doit être évaluée l'intégrale.

Pour calculer

 

$$\mathcal{I}_h^k(x) \int\!\!\!\int_{\widetilde{K}} \Phi(x,y)\ \widetilde{n}(y) \mbox{ d}s(y)$$ =

$$\int\!\!\!\int_\sigma \Phi(x,\widetilde{m}_k(s,t)) \ (\partial_s ... ...tilde{m}_k(s,t) \wedge \partial_t \widetilde{m}_k(s,t)) \mbox{ d}s \mbox{ d}t ,$$ = (20)

nous nous donnons un paramètre $\delta$ et nous distinguons les deux cas suivants:

- Si Si $dist(x,\widetilde{K}) > \delta$, alors l'intégrant est considéré comme étant régulier (on veut dire par là que ni l'intégrant, ni ses dérivées n'ont de fortes variations ou ne deviennent très grands). Dans ce cas, on utilise la formule de quadrature à trois points, du second ordre,

$$\int_\sigma f(s,t) \ \textrm{d}\, s \textrm{d}\, t= \frac{1}{6} \sum_{j=4}^{6} f(\rho_j)$$ (21)

où les nuds de quadrature $\rho_j$ sont les milieux des arêtes du triangle de référence $\sigma$.

- Si $0 < dist(x,\widetilde{K}) < \delta$, alors l'intégrant est de type "quasi-singulier"; lorsque le point x est pris proche d'un triangle $\widetilde{K}$, l'intégrant devient très grand et varie dans de grandes proportions en fonction du point $y \in \widetilde{K}$. La méthode d'intégration doit être à même de prendre en compte ce comportement quasi-singulier de l'intégrant. Nous avons recours à la formule de quadrature à sept points, d'ordre 5, référencée (T2:5-1) dans [#!Stroud71!#] p. 314,

 

$$\int_\sigma f(s,t) \ \textrm{d}\, s \textrm{d}\, t \frac{9}{80} f(\frac{1}{3},\frac{1}{3}) + B \left\{ f(\alpha,\alpha) + f(\alpha,\beta) + f(\beta,\alpha)\right\}$$ =

$$+ C \left\{ f(\gamma,\gamma) + f(\gamma,\zeta) + f(\zeta,\zeta)\right\} ,$$ (22)

où

$$\left\{ \begin{array}{ccc} \displaystyle \alpha = \frac{6-\sq... ...C = \frac{155+\sqrt{15}}{2400} .\\ &&\\ \end{array}\right .$$

La stratégie d'intégration consiste à subdiviser $\sigma$ en sous-triangles à mesure que la distance de x à $\Sigma$diminue et à approcher l'intégrale sur chacun de ces triangles par la formule de quadrature (23). Pour cela on se fixe N_d un nombre maximum de subdivision de $\sigma$et un réel positif $\varrho$. On distingue alors les cas suivants:

- Si $dist(P,\widetilde{K}) \leq \varrho$ alors on prend un paramètre de subdivision l=N_d.

- Si $\varrho < dist(P,\widetilde{K}) \leq 2\ \varrho$ alors on prend un paramètre de subdivision $l= \max\{N_d - 1, 0\}$.

- On poursuit ainsi jusqu'à ce que la distance $dist(P,\widetilde{K})$soit suffisamment importante pour que l'on puisse utiliser la formule à trois nuds (22).

Les valeurs des paramètres N_d et $\varrho$ sont fixées de manière empirique.