S'abonner à un flux RSS
 

Conditions aux limites

De Wikhydro

Conditions aux limites

L’onde incidente $ \Psi^{inc}\, $ est une onde plane régulière (c’est-à-dire monochromatique et monodirectionnelle) de période T, de direction $ \alpha_{inc}\, $par rapport à l’axe O x et d’amplitude $ \Psi_0=\frac{gH^{inc}}{2\omega} $$ H^{inc}\, $est la hauteur de la houle incidente. L’onde incidente s’écrit donc :

$ \Psi^{inc}=\Psi_0 e^{jk(x\cos \alpha_{inc}+y\sin \alpha_{inc})} $              (2)

La frontière du domaine est formée de frontières de deux types principaux (voir figure 4) :

       -Les frontières réfléchissantes constituées essentiellement par les ouvrages portuaires ;

       -Les frontières ouvertes en entrée ou en sortie.

Fig. 4 : Port schématique.

Frontière réfléchissante

Les frontières réfléchissantes sont caractérisées d’une part par un angle d’incidence $ \theta\, $ par rapport à la normale à la frontière et d’autre part par un coefficient de réflexion K R qui est supposé en général indépendant de l’angle d’incidence. La condition aux limites s’écrit :

$ \frac{\partial\Psi}{\partial n}=jk\cos \theta\frac{1-K_R}{1+K_R}\Psi $              (3)

Calcul des angles d’incidence

Si la frontière se trouve en zone éclairée (voir figure 4), il est recommandé de prendre l’angle $ \alpha_inc\, $ de la houle incidente. Si la frontière se trouve en zone d’ombre (voir figure 4), il est recommandé de choisir le calcul de l’angle $ \theta\, $ à partir d’un point source qui est un point de diffraction. Le calcul des points de diffraction peut être effectué à partir d’un algorithme récursif schématisé par la figure 5.

Fig. 5 : Algorithme récursif pour la détermination des angles d’incidence.

Algorithme [A B] : Tous les points de la frontière entre les points A et B sont initialement rattachés au point de diffraction A. En parcourant cette frontière dans le sens trigonométrique, on détermine deux points de diffraction A1 et A2. Les contours [A A1], [B1 A2] et [B2 B] restent rattachés au point de diffraction A. Pour les contours [Ai Bi] avec i variant ici de 1 à 2, on applique récursivement l’Algorithme [Ai Bi].

L’application de l’algorithme récursif sur le Port de Brest est présentée Figure 18.

Méthodes itératives

La direction locale de propagation de l’onde de houle peut être exprimée à partir de la fonction de phase S(x,y) (Isaacson et Qu, 1990) :

$ \alpha (x,y)=\arctan \left[\frac{\frac{\partial S(x,y)}{\partial y}}{\frac{\partial S(x,y)}{\partial x}}\right] $               (4)

Cette expression est rigoureuse seulement pour une onde plane se propageant sur un fond plat et une frontière parfaitement absorbante (KR=0). La méthode itérative est la suivante (Beltrami et al., 2001) : A l’itération i, une valeur $ \theta\, $i est proposée sur chacun des éléments de la frontière. La résolution de l’équation de pente douce donne la valeur des phases Si+1 sur la frontière. L’équation (4) permet le calcul de la direction locale $ \alpha\, $i+1. On en déduit l’angle d’incidence à la frontière $ \theta\, $i+1. Pour une frontière réfléchissante Steward et Panchang (2001) présentent différentes méthodes pour estimer la direction des ondes incidentes.


       1. Prenant la partie imaginaire de l’équation (3), on obtient une relation qui permet d’obtenir l’angle d’incidence $ \theta\, $. Cette expression est satisfaite exactement par la solution analytique :

$ \theta=\arccos \left[\frac{1}{k}\frac{1+K_n}{1-K_n}\frac{\partial S(x,y)}{\partial n}\right] $              (5)

       2. En substituant $ \theta\, $, obtenu grâce à l’équation (4) dans (3), une condition aux limites non-linéaire est obtenue. Les simulations numériques montrent que cette méthode ne reproduit pas la solution analytique lorsque le coefficient de réflexion KR n’est pas nul.

       3. Isaacson et al. (1993) utilise la dérivée tangentielle de la phase :

$ \theta=\left[\arcsin \frac{1}{k}\frac{\partial S}{\partial \tau}\right] $              (6)

Cette formulation ne dépend pas du coefficient de réflexion KR . La méthode itérative ne converge pas lorsque $ \theta\ge 40^\circ $.

       4. Steward et Panchang (2001) indiquent que les méthodes 1 et 3 ne sont applicables que pour des incidences proches de la normale et la méthode 2 pour des frontières parfaitement absorbantes. C’est la raison pour laquelle ils ont présenté une nouvelle méthode. En combinant les équations (5) et (6), on obtient :

$ \cos \theta \frac{\partial S}{\partial n}-\sin \theta \frac{\partial S}{ \partial n}= k \cos \theta \sin \theta \frac {2K_R}{1+K_R} $              (7)

       Tous calculs faits, on obtient :

$ \theta = \arctan \left[\frac{\frac{\partial S}{\partial \tau}}{\frac{\partial S}{\partial n}+k \frac{2K_R}{1+K_R}\cos \theta}\right] $              (8)

Pour une condition parfaitement absorbante (KR=0), on retrouve l’équation (4). Des estimations de l’angle $ \theta\, $peuvent être obtenues de manière itérative à partir de l’expression implicite (8).

Coefficient de réflexion

Le coefficient de réflexion KRest exprimé à partir de différentes formules empiriques en fonction du type d’ouvrages.

Cas des pentes lisses : le coefficient de réflexion KR peut être calculé à l’aide de la formule de Battjes (1974) :

$ \begin{cases} K_R=0,1\xi_p^2 \quad si K_R<1 \\ K_R=1\qquad sinon\end{cases} $              (9)

$ \xi_p\, $, aussi appelé paramètre d’Iribarren $ I_r\, $est égal à $ \frac{\tan \alpha}{\sqrt{S_{op}}} $. Le coefficient $ \alpha\, $ est l’angle de la pente de l’ouvrage et $ S_{op} $est la cambrure au large. Seelig (1983) propose la formule suivante :

$ K_R=\frac{a\xi_p^2}{b+\xi_p^2} $              (10)

avec un choix recommandé a=1 et b=5,5.

Cas des pentes recouvertes d’un enrochement : le fait d’introduire une rugosité sur la carapace contribue à dissiper une partie de l’énergie incidente et diminue le coefficient de réflexion. L’influence de cette rugosité est modélisée par le paramètre a de la formule de Seelig (1983). Ce paramètre prend la valeur 1 pour les pentes lisses, mais devient, pour les pentes rugueuses et poreuses, une fonction des deux paramètres suivant$ \sqrt{\frac{\partial}{\lambda}}\cot \alpha $ : et de $ \frac{H_{inc}}{H_d} $ avec $ \delta\, $ le diamètre des enrochements, $ \lambda\, $ la longueur d’onde, $ H_{inc}\, $la hauteur de la houle incidente et$ H_d\, $ la hauteur de déferlement définie par Goda (1975) avec $ h_s\, $ la profondeur en pied d’ouvrage :

$ H_d=0,17\lambda \Big(1-\exp {\Big[-4,712\frac{h_s}{\lambda}\Big]}\Big) $              (11)

Fig. 6 : Valeurs du coefficient '''a''' pour des pentes recouvertes d’un enrochement.

Les abaques donnant a se trouvent dans Seelig et Ahrens (1981) ou Benoit (1994).

Cas des digues à talus et des enrochements : la modélisation des digues à talus est plus complexe car la dissipation d’énergie intervient à plusieurs niveaux (carapace, sous-couche, noyau). Néanmoins, il est proposé d’adopter encore l’expression (10) avec a=0,16 et b=6,6.

Le coefficient de réflexion KR est un nombre complexe même s’il est souvent pris réel dans les codes de calcul d’agitation de houle. A partir d’un grand nombre de mesures expérimentales 2D et 3D sur les structures côtières, Sutherland et O’Donoghue (1998) donnent une formule empirique pour la phase $ \lambda\, $ en fonction du paramètre 3D $ \chi_3\, $ :

$ \gamma=-8,84\pi\chi_3^{1,25} $              (12)

Le paramètre 3D$ \chi_3\, $ dépend de la profondeur $ h_3\, $ au pied de la structure, de la pente $ \alpha\, $ de la structure, de la période de la houle T et de l’angle d’incidence $ \theta\, $

$ \chi_3= \frac{1}{\tan \alpha}\sqrt{\frac{h_s\cos \theta}{gT^2}} $              (13)

Couche dissipative

Comme vu plus haut avec l’équation (4), l’approche par coefficient de réflexion demande la connaissance du coefficient complexe KR, de la période de la houle T et de la direction $ \theta\, $ d’incidence de la houle sur la structure qui est a priori inconnue. Il est donc intéressant de pouvoir proposer une formulation qui n’utilise pas cet angle. C’est le cas du modèle à dissipation qui sera présenté plus tard. Nous allons établir ici la correspondance entre un coefficient de réflexion KR et une couche dissipative caractérisée par son coefficient de dissipation f et son épaisseur $ \mathit{l}\, $ et placée devant une paroi parfaitement réfléchissante.

Fig. 7 : Couche dissipative à la frontière.

Houari (2002) montre que le coefficient de réflexion devant la couche dissipative est donné par l’expression suivante :

$ K_R=\frac{K_0\cos \theta+j\beta_d k_y\tan (k_yl)}{K_0\cos \theta-j\beta_d k_y\tan (k_yl)} $              (14)

avec $ k_y=\sqrt{k^2-k_0^2\sin ^2\theta} $ et $ \beta_d=\frac{gF_0}{cc_g} $. $ F_0\, $ est une expression complexe dépendant du coefficient de dissipation $ \mathit{f}\, $ et de l’épaisseur de la couche poreuse dissipative. Si l’on choisit une couche non dissipative, alors $ \beta_d=1\, $ et l’on vérifie que $ \left[K_R\right] $est égal à 1. Choisissant une digue poreuse occupant toute la profondeur, Houari trace les variations du coefficient $ K_R\, $ en fonction de l’épaisseur et de l’angle d’incidence (voir figure 8 et 9 avec le paramètre $ \alpha = \frac{\omega^2h_s}{g} $). Le coefficient de réflexion est stabilisé pour une épaisseur de trois longueurs d’onde. En ce qui concerne l’angle d’incidence, le coefficient de réflexion diminue jusqu’à une direction comprise en 60 et 80 degrés puis remonte jusqu’à 1 lorsque l’on approche 90 degrés, soit une incidence rasante. Ce résultat théorique est en accord avec les résultats expérimentaux (Benoit, 1994).

Fig. 8 : Variations du coefficient de réflexion $ K_R $en fonction de l’épaisseur relative l/$ \lambda $ de la digue ($ \alpha=0,1 $;$ \theta=0 $ et $ \mathit{f}=1 $)

Fig. 9 : Variations du coefficient de réflexion $ K_R $ en fonction de l’angle d’incidence ($ \alpha=1 $;$ \mathit{f}=1 $;l/$ \lambda=0,25 $)

Frontière ouverte

Les méthodes pour les solutions de problème d’onde dans des domaines infinis ont été développées depuis les années 1970 (Givoli, 1992). Elles sont utilisées dans de nombreux domaines d’application comme l’acoustique, l’électromagnétisme, l’hydrodynamique ou la géophysique. Trois grands types de méthode ont émergé : les méthodes exactes avec équation intégrale, les conditions de radiation et les couches dissipatives.

Equation intégrale : solution analytique avec frontière circulaire

Fig. 10 : Domaine de calcul avec frontière circulaire (gauche) et module de la hauteur de houle – houle d’Ouest avec une île circulaire (droite).

Certaines conditions aux limites ont pour handicap de demander un angle d’incidence que l’utilisateur ne connaît pas a priori. Il est possible de s’en affranchir si l’on considère que le domaine $ \Omega_2\, $ extérieur au domaine d’étude est de profondeur constante h0, soit de nombre d’onde constant k0. Le problème s’écrit alors sous la forme d’une équation intégrale. Si l’on choisit en plus une frontière ouverte circulaire (exemple d’une île), l’équation intégrale peut s’exprimer analytiquement sous la forme d’une série de Fourier.

$ \frac{\partial\Psi}{r}\Big(R,\theta\Big)=k_0\sum_{n=0}^\mathcal{1}\int_{0}^{2\Pi}m_n \cos \Big[n(\theta-\theta')\Big]\Psi(R,\theta')\,d\theta' +k_0\Psi_0 \sum_{n=0}^\mathcal{1}h_n\cos \big[n\theta\big] $

avec $ \delta_{ij}\, $ le symbole de Kronecker et les valeurs de coefficients suivants calculés à l’aide des fonctions de Bessel et de Hankel de première espèce :


$ \begin{cases}m_n=\frac{1}{\pi(1+\delta_{n0})}\frac{{H'}_n^1(k_0R)}{H_n^1(k_0R)} \\ h_n=j^n(2-\theta_{n0}\Big[{J'}_n(K_0R)-J_n(K_0R)\frac{{H'}_n^1(k_0R)}{H_n^1(k_0R)}\Big]\end{cases} $              (14)

Equation intégrale : solution analytique sur une secteur angulaire

Fig. 11 : Domaine de calcul avec frontière circulaire (gauche) et phase de la hauteur de houle – houle d’Est avec une digue semi-infinie (droite).

Nous considérons dans cette section un port schématique comprenant deux jetés rectilignes formant un angle $ \theta_1\, $ entre elles (voir figure 11), considérées semi-infinies et parfaitement réfléchissantes. Ces deux jetés définissent un secteur angulaire dont le sommet est noté O. Le domaine d’étude $ \Omega_1\, $ est fermé par un arc de cercle de centre O, de rayon R s’appuyant sur les deux jetés. La profondeur, dans le domaine $ \Omega_1\, $est variable et notée h(x). On considère encore que le domaine $ \Omega_2\, $extérieur au domaine d’étude, est de profondeur constante $ h_0\, $, soit de nombre d’onde constant $ k_0\, $. La sollicitation est une onde plane harmonique incidente $ \Psi^{inc}\, $ de module $ \Psi_0\, $ au point O et de direction $ \delta_0\, $.

$ \Psi^{inc}= \Psi_0 e^{jr\cos (\delta-\delta_0)} $              (15)

Sergent et al. (2002) ont exprimé l’équation intégrale qui régit ce type de conditions aux limites à partir d’une série de modes propres d’un secteur angulaire :

$ \frac{\partial\Psi}{\partial r}(R\theta)=k_0 \sum_{n=0}^\mathcal{1}\int_{0}^{\theta_1}m_n\cos [v_n\theta]\cos [v_n\theta']\Psi(R,\theta')\,d\theta' + 2\mu k_0 \Psi_0\sum_{n=0}^\mathcal{1}h_n \cos [v_n\theta]\cos [v_n\theta'] $               (16)

avec $ \mu=\frac{\pi}{\theta_1} $ et $ v_n = n \mu \, $ les valeurs de coefficients suivants calculés à l’aide des fonctions de Bessel et de Hankel de première espèce :

$ \begin{cases}m_n=\frac{2}{\theta_1(1+\delta_{n0})}\frac{{H'}_{v_n}^1(k_0R)}{H_{v_n}^1(k_0R)} \\ h_n=e^{jv_n}(2-\theta_{n0}\Big[{J'}_{v_n}(K_0R)-J_{v_n}(K_0R)\frac{{H'}_{v_n}^1(k_0R)}{H_{v_n}^1(k_0R)}\Big]\end{cases} $              (17)

Equation intégrale : fonction de Green avec frontière quelconque

Une autre méthode sur une frontière quelconque $ \partial\Gamma $consiste à écrire l’équation intégrale suivante :

$ \frac{1}{2}\Psi(x_0)=\int_{\partial\Gamma}\Big(\Psi(\underline{x})\frac{\partial G}{\partial n}(\underline{x},\underline{x}_0)- G(\underline{x})\frac{\partial\Psi}{\partial n} (\underline{x},\underline{x}_0)\Big)\,d\underline{x} $              (18)

avec $ \underline{x}_0 \in \partial\Gamma\, $ La fonction de Green $ G (x,y)\, $ est égale à $ -\frac{j}{4}H_0^{(1)}\Big(k\left| \underline{x}-\underline{y}\right| \Big) $ La dérivée normale $ \frac{\partial G}{\partial n}(\underline{x},\underline{y}) $ est égale à $ \frac{j}{4}H_0^{(1)}\Big(k\left| \underline{x}-\underline{y}\right| \Big)\frac{\underline{x}-\underline{y}}{\left|\underline{x}-\underline{y}\right| }\underline{n} $ avec $ H_0^{(1)} $ et $ H_1^{(1)} $ les fonctions de Hankel de première espèce d’ordre 0 et 1.

Les conditions de radiation : conditions aux limites non-réfléchissantes ou absorbantes

Une frontière ouverte est une frontière absorbante. Elle peut être considérée comme une frontière réfléchissante dont le coefficient de réflexion KR est nul. A partir de l’équation (3), la condition aux limites peut s’écrire :

$ \frac{\partial \Psi}{\partial n}= jk \cos \theta \Psi $              (19)

Il est recommandé de prendre l’angle $ \theta\, $de l’onde incidente égal à zéro, c’est-à-dire de considérer que l’onde est sortante suivant la normale à la frontière. L’équation (19) devient alors :

$ \frac{\partial \Psi}{\partial n}= jk \Psi $                (20)


Une méthode utilisée communément pour les problèmes de propagation d’onde est d’écrire l’approximation parabolique pour l’onde qui est radiée à travers la frontière ouverte. Avec R le rayon de courbure local, avec un domaine de profondeur constante, l’approximation parabolique en coordonnées curvilignes s’écrit (Xu et al., 1996) ou (Panchang et al., 2000) :

$ \frac{\partial \Psi}{\partial n}(R,\theta)= \Big[jk-\frac{1}{2R}+\frac{j}{8kR^2}\Big]\Psi(R,\theta)+\frac{j}{2KR^2}\frac{\partial^2\Psi }{\partial\theta^2}(R,\theta) $               (21)

Cette formulation a l’avantage d’être locale et donc facile à implémenter. D’un autre côté, l’approximation parabolique s’appuie sur l’hypothèse que l’onde est diffractée selon la normale à la frontière. Mieux cette hypothèse est respectée, meilleure sera la précision de la formulation. Pour un rayon de courbure infini, l’équation (21) devient en coordonnées cartésiennes :

$ \frac{\partial \Psi}{\partial x}(x,y)= jk\Big[1+\frac{1}{2k^2}\frac{\partial^2}{\partial y^2}\Big]\Psi(x,y) $              (22)

On démontre que, pour la formulation (20), la réflexion d’une onde incidente sur la frontière est égale à :

$ K_R=\left|\frac{1-\cos \theta}{1+\cos \theta}\right| $              (23)

Pour la formulation (22) , elle est égale à :

$ K_R={\left|\frac{1-\cos \theta}{1+\cos \theta}\right|}^2 $              (24)

Couches dissipatives ou absorbantes

L’équation de pente douce avec dissipation W s’écrit (Dingemans, 1997) :

$ \nabla\big(cc_g\nabla\Psi\big)+k^2cc_g\Psi=-j\omega W \Psi $              (25)

$ {\tilde{k}}^2=k^2+j\omega\frac{W}{cc_g} $              (26)

Fig. 12 : Couche dissipative

A la frontière entre le domaine de calcul et la couche dissipative, la loi de Snell-Descartes s’applique :

$ k \sin \theta = \tilde{k} \sin \tilde{\theta} $              (27)

On démontre que le module du coefficient de réflexion à la frontière est égal à :

$ K_R=\frac{\sin \big(\theta - \tilde{\theta}\big)}{\sin \big(\theta + \tilde{\theta}\big)} $              (28)

On remarque que, pour avoir un coefficient de réflexion le plus faible possible, il faut diminuer le coefficient de dissipation W et donc élargir la couche dissipative. Une alternative consiste à augmenter progressivement ce coefficient. Wei et Kirby (1995) proposent une loi exponentielle :

$ W=W_0\frac{\exp {\big(\frac{x-x_s}{x+x_s}\big)}^n-1}{\exp {(1)^n}} $              (29)

avec une couche dite éponge située entre les plans $ x = x_s \, $ et $ x = x_1\, $

Outils personnels