Le changement climatique est devenu une préoccupation des citoyens du monde et un enjeu planétaire, imposant aux gouvernements d’intégrer la quantification de ses effets dans les politiques publiques. Plus particulièrement, l’exhaussement du niveau des océans constitue un enjeu majeur pour les populations situées le long des côtes. D’après l’OCDE « Il est désormais quasi certain que le niveau des mers augmentera d’au moins un mètre, et cette hausse pourrait même intervenir dans les 80 prochaines années selon certaines modélisations, avec à la clé des conséquences graves qui causeront des dommages aux infrastructures, la disparition de terres et des déplacements de populations ».
Le thème d’étude est la quantification de l’impact des houles sur le littoral. La probabilité de survenue d’événements extrêmes va considérablement augmenter et leurs impacts sur le littoral provoquer des désordres sur les secteurs les plus vulnérables entrainant le recul du trait de côte et la disparition de propriétés et d’ouvrages.
Afin de quantifier ces effets, les bureaux d’études ont besoin d’outils de simulation capables de travailler à grande échelle, rapides d’exécution et précis, qui permettent de dégrossir les problèmes et d’identifier des zones spécifiques où des outils de modélisation numériques beaucoup plus sophistiqués, et à faible résolution spatiale viendront préciser les impacts.
Le modèle le plus usité qui permet de modéliser les houles est le modèle de Berkhoff qui intègre la représentation des phénomènes de réfraction des fonds, le shoaling, la diffraction et la réflexion des structures côtières comme les digues ou les jetées [1].
Vidéo d'introduction par Wikhydro
Modélisation de la houle par homotopie : une brève introduction
Modèle de Berkhoff
Ce modèle a pour expression (en bi-dimensionnel):
$ \nabla.(CC_g\nabla \phi)+k^2CC_g\phi=0 $
où $ \phi $ est le potentiel, k est le nombre d’onde fonction de la profondeur H et de la fréquence $ \omega $ (T est la période), par la relation implicite $ \omega^2=gk \tanh(kH) $ , C est la célérité de l’onde, Cg est la célérité de groupe des vagues.
Pour simplifier le problème, nous nous placerons dans le domaine des ondes longues, ce qui signifie que $ C=C_g=\sqrt{gH} $.
L'évolution dans le temps de la hauteur de houle est donnée par : $ h(x,t)=\Re \left (\phi e^{-i\omega t} \right ) $.
Homotopie
Considérons l’équation différentielle suivante à résoudre :
$ L(u(x)=f(x)+N(u(x), x \in \Omega $
avec les conditions limites: $ B(u,u_n=0), x \in \Gamma $
où L est un opérateur linéaire, N un opérateur non-linéaire et f les termes complémentaires de l’équation.
L'homotopie de Liao [9] est définie de la manière suivante :
$ (1-p)\left[ L(U(x,t);p)-L(u_0(x,t)) \right ] +H(p)\left[ L(U(x,t);p)-N(U(x,t),p)-f(x) \right] $
où $ p \in \left[0,1\right] $ est un paramètre qui varie de 0 à 1, u0 est une estimation initiale de la solution.
Lorsque p=0 nous retrouvons la solution initiale U=u0 et lorsque p=1, nous retrouvons la solution exacte.
La transformation de p de 0 à 1 qui conduit U(x) de la solution estimée à la solution exacte provient de la transformation de U(x) en série de Taylor
$ U(x)=u_0(x)+\sum_{m=1}^\infty u_m(x)p^m $ avec $ u_m(x)=\dfrac{1}{m!}\dfrac{\partial^mU(x,t;p)}{\partial p^m}\Bigr]_{p=0} $
La mise en œuvre de l’homotopie nécessite de faire des hypothèses sur la fonction H(p) qui permet notamment un ajustement des non-linéarités.
Pour la méthode HAM « Homotopie Perturbation Method », que nous mettrons en œuvre: $ H(p)=1 $.
Objectif des travaux
L’objectif du présent projet est de résoudre cette EDP par une méthode semi-analytique : l’homotopie, assortie des conditions limites et des conditions initiales ad-hoc fonction du domaine étudié et du régime de houle de projet.
L’intérêt de la méthode est qu’elle permet de partir d’une solution connue relativement simple et de converger vers une solution complète. L’idée maitresse est d’introduire un paramètre p qui varie entre 0 et 1 et d’assurer une déformation continue entre une première estimation de la solution relativement simple (p=0) et la valeur finale de la solution (p=1) du système d’équations à résoudre. L’un des grands intérêts de la méthode est qu’elle permet d’atteindre la solution finale avec peu de termes. Par ailleurs, elle utilise la résolution formelle et peut donc être facilement programmée à l’aide d’outils de calcul formel tels que WSMAXIMA que nous utiliserons ici.
Travaux demandés
Différents cas concrets de géométries et de conditions aux limites devront être étudiés :
Mise en œuvre de la méthode d'homotopie sur une déformation de courbe
Chaque trinôme choisira 2 courbes (solution initiale et solution finale et mettra en œuvre la méthode d'homotopie pour passer de l'une à l'autre comme l'exemple présenté par BRETON-NGUYEN-SALLES, de manière à bien comprendre la dynamique de la méthode d'homotopie.
Cas n°1
Canal monodimensionnel plat de longueur L avec entrée par l'aval d'une onde de fréquence unitaire $ \phi=1 $ (condition de Dirichlet) et sortie libre amont $ \phi_x=ik\phi $ (condition de Robin).
Pour davantage de facilité de compréhension et de mise en œuvre, il est suggéré d'imposer ces 2 conditions limites en x=0. On peut montrer que ceci est équivalent au cas où les 2 conditions limites sont imposées aux 2 extrémités du domaine (x=0 et x=L).
Attention : Ceci n'est valable que pour ce cas bien particulier.
Conseils de mise en œuvre
La relation d'homotopie s'écrit en choisissant la dérivée seconde comme fonction auxiliaire linéaire et en partant d'une solution initiale nulle:
$ (1-p)\phi_{xx}+p(\phi_{xx}+k^2\phi)=0 $
En injectant la décomposition en série entière $ \phi(x,p)=\phi_0(x)+p\phi_1(x)+p^2\phi_2(x)+p^3\phi_3(x)+... $ et sa seconde dérivée:
$ \phi_{0,xx}(x)+p\phi_{1,xx}(x)+p^2\phi_{2,xx}(x)+p^3\phi_{3,xx}(x)+... $
Nous obtenons:
$ (1-p)(\phi_{0,xx}(x)+p\phi_{1,xx}(x)+p^2\phi_{2,xx}(x)+p^3\phi_{3,xx}(x)+...)+p[\phi_{0,xx}(x)+p\phi_{1,xx}(x)+p^2\phi_{2,xx}(x)+p^3\phi_{3,xx}(x)+...+k^2(\phi_0(x)+p\phi_0(x)+p^2\phi_2(x)+p^3\phi_3(x)+...)]=0 $
Il faut ensuite simplifier et écrire cette relation suivant les puissances de p croissantes. Cette relation étant valable quel que soit p, tous les coefficients devant les puissances de p sont donc nuls.
$ \phi_{0,xx}=0 $ soit $ \phi_0=Ax+B $.
Introduisant les conditions limites suivantes: $ \phi_0^0=1 $ et $ \phi_{0,x}^L=ik\phi_{0}^L $, il vient :
$ \phi_0=1+\dfrac{ik}{1-ikL}x $.
$ \phi_{1,xx}-\phi_{0,xx}+\phi_{0,xx}+k^2\phi_0=0 $ soit $ \phi_{1}=-k^2\int\phi_0dxdx +Ax+B $.
Introduisant les conditions limites suivantes: $ \phi_1^0=0 $ et $ \phi_{1,x}^L=ik\phi_{1}^L $,
il vient : $ \phi_1=- \dfrac{k^2L(k^2L^2+3ikL-3} {3(1-ikL)^2}x-k^2\left (\dfrac{ik} {6(1-ikL)}x^3+\dfrac{1} {2}x^2 \right ) $
A vous d'écrire la suite et de mettre en œuvre WXMAXIMA (voir un exemple de code source ci-dessous) avec les valeurs numériques suivantes:
$ k=\dfrac{1} {100} $ (nombre d'onde en m-1), $ H=40 $ (profondeur en m), $ c=\sqrt{gH} $ (célérité de l'onde en m/s),$ \lambda=\dfrac{2\pi}{k} $ (longueur d'onde en m) ,$ L=2\lambda $ (longueur du domaine en m).
- montrer que la solution par homotopie converge vers la solution analytique (choisir kL<<1)
- faire une animation de la solution dans le temps en incrémentant le nombre d'ordres d'homotopie
Cas n°2
Domaine monodimensionnel plat de longueur L avec entrée par l'aval d'une onde de fréquence unitaire et une condition de flux aval
$ \phi_x=ik(2-\phi) $ et réflexion totale amont $ \phi_x=0 $ .
- 3 trinômes utiliseront la méthode standard qui fonctionne très bien.
- 3 trinômes utiliseront la méthode de l’image virtuelle (superposition onde incidente en provenance de l’aval et onde réfléchie en provenance de l’amont) sur la base du cas n°1.
- montrer que la solution par homotopie converge vers la solution analytique
- faire une animation de la solution dans le temps en incrémentant le nombre d'ordres d'homotopie
Calcul explicite de $ \phi_1(x) $
$ \phi_1(x)=-k^2\int\phi_0(x)dxdx+Ax+B $ d'où: $ \phi_{1,x}(x)=-k^2\int\phi_0(x)dx+A $
Condition limite en x=0 :
- $ \phi_{1,x}(x=0)=-ik \phi_{1,x}(x=0) $ soit $ \left|-k^2\int\phi_0(x)dx+A\right|_{x=0}=ik\left|-k^2\int\phi_0(x)dxdx+Ax+B\right|_{x=0}==>A=-ikB $ (vrai quel que soit l'ordre)
Condition limite en x=L :
- $ \phi_{1,x}(x=L)=0 $ soit $ \left|-k^2\int\phi_0(x)dx+A\right|_{x=L}=0 ==>A=k^2\int\phi_0(x)dx==>A=2k^2L $
on en déduit: $ \phi_1(x)=-k^2x^2+2k^2Lx+2ikL $
... à vous de poursuivre les calculs et mettre en oeuvre le programme sous WSMAXIMA.
Cas n°3
Domaine monodimensionnel de longueur L avec pente du fond constante (s=cste) avec entrée par l'aval d'une onde de fréquence unitaire $ \phi=1 $ et sortie libre amont $ \phi_x=ik\phi $.
Nota : on représentera le terme de pente par un développement de Mac-Laurin en mettant en évidence la précision obtenue en fonction du nombre de terme pris en compte.
2 cas sont à étudier :
- 3 trinômes : on considèrera que $ k=k_0=cste $
- 3 autres trinômes : on considèrera que $ k=k_0\sqrt{\dfrac{H_0}{H(x)}} $
Dans chaque cas, les travaux à réaliser sont les suivants:
- détermination de la solution analytique pour chacun des cas ci-dessus, de l'équation de base transformée par le changement de variable $ z=H0-sx $. Cette équation est une une équation de type Bessel.
- résolution de cette équation par la méthode d'homotopie en utilisant les conditions limites du cas n°1
- superposition de la solution analytique et de la solution par homotopie pour les valeurs numériques ci-dessus et pour une pente du fond $ s=\dfrac{1}{200} $
- animation temporelle de la solution
- étude de sensibilité de la solution en faisant varier la valeur de k
Cas n°4
Vague sphérique générée par une source périodique sinusoïdale.
Nous traitons ici de l'évolution de la surface libre dans un domaine infini en grande profondeur. La source ponctuelle est appliquée autour d'un cercle de rayon r_0 centré sur un domaine circulaire de rayon R qui laisse sortir librement cette onde en r=R.
l'équation de Berkhoff se simplifie alors en équation de Helmholtz et s'exprime en coordonnées polaires avec les conditions suivantes:
$ \begin{cases} \Delta \phi + k^2\phi=0, \\ \phi^{r=r_0}=1, \\\phi_r^{r=R}=ik\phi^{r=R}. \end{cases} $
En coordonnées polaires, la relation ci-dessus s'écrit de manière simplifiée étant donné que le problème est caractérisé par une symétrie de révolution, donc est indépendant de $ \theta $.
$ \begin{cases} \phi_{rr}+\dfrac{1}{r}\phi_r + k^2\phi=0, \\ \phi^{r=r_0}=1, \\\phi_r^{r=R}=ik\phi^{r=R}. \end{cases} $
On choisira les valeurs numériques suivantes: $ r_0=1m $, $ R=100m $ et $ k=0.1m^{-1} $.
Les résultats sont affichés suivant un rayon de r0 à R:
- le premier graphique superpose la solution analytique avec la solution par homotopie
- le second consiste en une application temporelle avec une fréquence de $ \omega=\sqrt{gk} $ car nous sommes en profondeur infiniée.
- on produira également une animation 2D si possible
Documents à produire
Rendez-vous 0 :
Le cas 1 sera traité en entier et il vous sera demandé de préparer le RV 1 complètement.
Rendez-vous 1:
Pour le RV1: il convient de produire un premier document devant détailler les points suivants:
Avoir traité en entier le cas n°1 sur WIKHYDRO
- La solution analytique
- La solution par homotopie
- Une superposition de la solution analytique et de la solution par homotopie sur la base des valeurs numériques fournies
- Une étude de sensibilité de la solution en fonction du nombre d'onde
- Une réflexion portant sur l'analyse des résultats, les limites et l'intérêt de la méthode d'homotopie
De manière à démarrer de manière efficace, il vous sera demandé de fournir également une première étude théorique du cas n°2, notamment:
- La solution analytique
- La solution par homotopie
Pour la fin du cours
Pour chaque cas ci-dessus, il conviendra de fournir :
- La solution analytique
- La solution par homotopie
- Une superposition de la solution analytique et de la solution par homotopie sur la base des valeurs numériques fournies
- Une étude de sensibilité de la solution en fonction du nombre d'onde
- Une réflexion portant sur l'analyse des résultats, les limites et l'intérêt de la méthode d'homotopie
Plus globalement,
- Le rapport sera publié sur le wiki du CEREMA : WIKHYDRO sous la forme d’une page wiki.
- Une ou plusieurs vidéos à caractère pédagogique seront élaborées explicitant les fondements de la méthode homotopie, les résultats obtenus complétés par des propositions d’application de cette méthode dans le contexte du changement climatique.
Bibliographie
- Panchang V.G., Cushman-Roisin B., Pearce B.R., “Combined Refraction-Diffraction of Short-Waves in Large Coastal Regions”, Coastal Eng., 12 (1988), 133-156
- Biazar J. Azimi F., “He’s Homotopy Perturbation Method for Solving Helmholtz Equation”, Int. J. Cotemp. Math. Sciences, Vol. 3, 2008, n° 15, 739-744
- Rajarama R., Hariharan G., “Solving Helmholtz equation by the homotopy perturbation method”, Int. Jour. of Math. and Computer Applications Research (IJMCAR) Vol. 2 Issue 3 Sept. 2012 70-75.
- Biazar J., Eslami M., “A new technique for non-linear two-dimensionnal wave equations”, Sciencia Iranica B (2013) 20(2), 359-363
- Chun C., Jafari H., Kim Y., “Numerical method for the wave and nonlinear diffusion equations with the homotopy perturbation method”, Comp. and Math. With Applications 57 (2009) 1126-1231
- Prakash A., Goyal M., Gupta S., “Numerical simulation of space-fractional Helmholtz equation arising in seismic wave propagation, imaging and inversion, Pramana – J. Phys. (2019) 93:28
- Jafari H. Momani S., 3Solving Fractional diffusion and wave equations by modified homotopy perturbation method”, Physics Letters A 370 (2007) 388-396
- Yin X-B, Kumar S, Kumar D., “A modified homotopy analysis method for solution of fractional wave equations”, Advances in Mechanical Engineering 2015, Vol. 7(12) 1-8
- Shijun Liao, “Homotopy Analysis Method in Nonlinear Differential Equations”, Springer, 549 p
Quelques conseils
- Chaque trinôme disposera d'une page WIKHYDRO accessible avec les codes qui vont être fournis lors du RV 0:
-
- MEDIAWIKI/WIKHYDRO dispose d'un tuto (voir en bas à gauche de la page lorsque vous êtes loggés uniquement)
- MEDIAWIKI/WIKHYDRO dispose également d'une aide LATEX pour écrire les formules mathématiques
N'hésitez pas à écrire de belles pages, bien documentées et agréables à lire.
Exemple et code source
Voici un exemple simple de la résolution par homotopie de l'équation de convection linéaire. Il s'agit d'un transport sans déformation d'une sinusoïde.
$ h_t+h_x=0 $
La relation correspondante de l'homotopie s'écrit:
$ (1-p)(h_t-Z_{0,t})+p( h_t+h_x)=0 $
Il suffit d'injecter la décomposition en série de $ h(x,t)=\sum_{n=0}^{N}(p^nh_n(x,t)) $ et de ses dérivées (N est le nombre de termes de la série)
- La solution analytique est donnée par la théorie des caractéristiques.
- Solution initiale : $ Z(x)=Z0+sin(kx) $
Le programme WXMAXIMA est le suivant: (il suffit de copier-coller ces sources dans le logiciel)
kill(all)$ /* initialisation */
/* Parametres du cas a etudier */
Z0:1$ /* amplitude de l'onde */
L:1 $ /* longueur du domaine monodimensionnel */
k:2*%pi/L$ /* nombre d'onde */
t0:1.1$ /* temps fin de calcul */
Niter:20$ /* nombre d'ordres homotopie */
/* initialisation des dimensions des tables*/
array(h,Niter)$ /* hauteur de houle */
/* Solution initiale*/
h[0]:Z0+sin(k*x)$
/* boucle sur le nombre d'ordres de l'homotopie */
for n:1 thru Niter do
(print("n = ",n),
hh:-integrate(diff(h[n-1],x),t),
h[n]:ratsimp(hh));
/* Hauteur de houle cumulée sur les ordres d'homotopie*/
HT:sum(ev(h[n],t=t0),n,0,Niter)$
/* solution analytique donnée par les caractéristiques */
SA:Z0+sin(k*(x-t0))$
/* Solution analytique discrétisée */
np:11$interval:L/(np-1)$ /* nb de points de discrétisation */
SAD:makelist([(i-1)*interval,Z0+sin(k*(interval*(i-1)-t0))],i,np)$ /* Solution analytique */
/* tracé des courbes*/
wxplot2d ([h[0],[discrete,SAD],HT], [x, 0, L],[y, 0, 25/10],[style,lines,points,lines],[color,blue,red,green],[point_type,times],[legend, "Sol initiale","Sol analytique à t=0.1s", "Sol Homotopie à t=0.1s"], [xlabel, "Abscisse (m)"], [ylabel, "Hauteur (m)"], [title, "Convection lineaire ; k=2*pi"] );
|
S'abonner à un flux RSS
