S'abonner à un flux RSS
 

Jean-Michel Tanguy/SujetENTPE2022/CALLOCH - DECHAVASSINE - DIANOUX : Différence entre versions

De Wikibardig
(Ordre 0)
(Ordre n)
 
(29 révisions intermédiaires par un utilisateur sont masquées)
Ligne 1 : Ligne 1 :
 
Pour retourner à la page énoncé du sujet : [[Utilisateur:Jean-Michel Tanguy/SujetENTPE2022]]
 
Pour retourner à la page énoncé du sujet : [[Utilisateur:Jean-Michel Tanguy/SujetENTPE2022]]
 
== Cadre de l'étude ==
 
== Cadre de l'étude ==
 +
===Contexte===
 +
Le changement climatique entraîne des bouleversements majeurs dans l'ensemble des conditions de vie de notre planète. Parmi eux, nous nous intéressons ici à l'élévation du niveau des océans, qui contribue au recul du trait de côte sur les littoraux. Ce phénomène est bien visible en France, comme l'illustre cette photographie:<br/>
 +
[[File:erosion_littoral_SudOuest.jpg|400px|érosion du trait de côte, Soulac-sur-Mer, Aquitaine, 3 février 2014|érosion du trait de côte, Soulac-sur-Mer, Aquitaine, 3 février 2014 ]]<br/>
 +
''Source:"Sud Ouest" / Julien Lestage''<br/>
 +
Cette érosion est notamment provoquée par la houle, mouvement ondulatoire de la surface de la mer provoqué par le vent loin des côtes.
 +
De nombreuses études sont menées sur ce phénomène dans des laboratoires de recherche et utilisent pour cela des canaux ou bassins à houle (comme au [https://m2c.cnrs.fr/poles-de-competences/mesocosmes-et-modelisation-analogique/ laboratoire M2C-CNRS/Univ. de Rouen/Univ. de Caen] https://m2c.cnrs.fr/poles-de-competences/mesocosmes-et-modelisation-analogique/).
 +
On peut en voir un exemple ici:<br/>
 +
[[File:canalHoule_OuestFrance.jpeg|400px|Canal à houle de l'ESITC]]<br/>
 +
''Source: Ouest-France/ESITC/EDINBURGH DESIGNS LTD - 15/01/2018''<br/>
  
On cherche à modéliser une houle dans quatre cas différents. Pour cela, on utilise le modèle de Berkhoff, dont l'équation est :<math>\nabla.(CC_g\nabla \phi)+k^2CC_g\phi=0</math><br />
+
Le schéma ci-dessous permet de comprendre les principales étapes de formation de ce phénomène:<br/>
 +
[[File:schemahoule_magicsurfschool.png|400px|Schéma des phases de la houle]]<br/>
 +
''Source: site internet magicsurfschool.com/blog/formation-des-vagues/''<br/>
 +
La démarche présentée par la suite repose sur la modélisation de la houle par deux méthodes de résolution différentes:<br/>
 +
*une méthode analytique, notamment grâce au logiciel MAXIMA pour les équations les moins évidentes à résoudre analytiquement
 +
*une méthode d'homotopie, utilisant la décomposition en série entière de la solution
 +
 
 +
===Modélisation===
 +
On cherche à modéliser une houle dans quatre cas différents. Pour cela, on utilise le modèle de Berkhoff, dont l'équation est :<br/>
 +
:::<span style="color: #26B260"><math>\nabla.(CC_g\nabla \phi)+k^2CC_g\phi=0</math></span><br />
 
où <math> \phi </math> est le potentiel, <math>k</math> est le nombre d’onde, <math>C</math> est la célérité de l’onde et <math>C_g</math> est la célérité de groupe des vagues.<br />
 
où <math> \phi </math> est le potentiel, <math>k</math> est le nombre d’onde, <math>C</math> est la célérité de l’onde et <math>C_g</math> est la célérité de groupe des vagues.<br />
 
On se place dans le domaine des ondes longues, ce qui entraîne que <math>C = C_g = \sqrt{gH}</math>, avec g l'accélération de la pesanteur et H la profondeur du canal. L'équation de Berkhoff devient ainsi : <br />
 
On se place dans le domaine des ondes longues, ce qui entraîne que <math>C = C_g = \sqrt{gH}</math>, avec g l'accélération de la pesanteur et H la profondeur du canal. L'équation de Berkhoff devient ainsi : <br />
:::<math>\nabla^2 \phi + k^2\phi = 0</math>. <br />
+
:::<span style="color: #26B260">'''<math>\nabla^2 \phi + k^2\phi = 0</math>'''</span>. <br />
 
L'évolution de la hauteur de la houle en fonction du temps est donnée par : <math>h(x,t)=\Re \left (\phi e^{-i\omega t} \right )</math><br />
 
L'évolution de la hauteur de la houle en fonction du temps est donnée par : <math>h(x,t)=\Re \left (\phi e^{-i\omega t} \right )</math><br />
  
Ligne 53 : Ligne 71 :
  
 
On peut visualiser l'évolution de la solution homotopique selon la valeur de p grâce au logiciel Maxima:<br />
 
On peut visualiser l'évolution de la solution homotopique selon la valeur de p grâce au logiciel Maxima:<br />
[[File:ezgif.com-gif-maker.gif|400px]] <br />
+
[[File:GIF Cas 1 final (1).gif|400px]]
Sur le graphique ci-dessus, la solution analytique est représentée par les croix rouges. La solution homotopique est représentée par la courbe verte. On illustre ici la variation du paramètre p, de 0 (solution en bleue) à 1 (où elle est la plus proche de la solution analytique). La méthode par homotopie permet bien de retrouver la solution analytique.
+
<br />
 +
Sur le graphique ci-dessus, la solution analytique est en rouge, la solution homotopique est représentée par la courbe verte. On illustre ici la variation du paramètre p, de 0 (solution en bleue) à 1 (où elle est la plus proche de la solution analytique). La méthode par homotopie permet bien de retrouver la solution analytique.
  
 
== Cas n°2 : Canal plat fermé ==
 
== Cas n°2 : Canal plat fermé ==
Ligne 65 : Ligne 84 :
  
 
==== Solution par homotopie ====
 
==== Solution par homotopie ====
 +
L’équation homotopique à résoudre s’écrit : <math>(1−p)\phi_{xx}+p(\phi_{xx}+k^2\phi)=0 </math>.
 +
La résolution est très similaire au cas 1, seules les conditions aux limites changent.
 
===== Ordre 0 =====
 
===== Ordre 0 =====
 
+
Elle devient donc <math>\phi_{0,xx}=0 \ soit \ \phi_0=Ax+B. </math>
 +
A l’aide des conditions aux limites, on obtient <math> \phi_0 = 2. </math>
 +
===== Ordre n =====
 +
A chaque ordre, on a <math> \phi_n(x)=−k^2\int \phi_{n-1}(x)dx+A_nx+B_n </math>
 +
On a les conditions aux limites suivantes : <math> \phi_{n,x}^0=−ik\phi_{n,x}^0 \ et \ \phi_{n,x}^L=0 </math>
 +
Il advient alors <math> A_n=(k^2 \iint \phi_{n-1}(x)dxdx)_{x=L} et B_n={i \over k}A_n. </math>
 +
On programme ensuite la résolution par récurrence sur Maxima.
  
 
== Cas n°3 : Canal en pente ==
 
== Cas n°3 : Canal en pente ==
  
 
===Résolution analytique===
 
===Résolution analytique===
On s'intéresse à présent à un canal dont le fond est en pente. Le phénomène étudié est toujours unidimensionnel, avec les mêmes conditions aux limites que pour le cas n°1. <br />
+
On s'intéresse à présent à un canal dont le fond est en pente. Le phénomène étudié est toujours unidimensionnel, avec les mêmes conditions aux limites que pour le cas n°1, soit: <br />
 +
*En aval, onde de fréquence unitaire <math>\phi=1</math> (condition de Dirichlet) et
 +
*En sortie libre amont <math>\phi_x=ik\phi</math> (condition de Robin).<br />
 
[[File:pente_canal.jpeg|400px]] <br />
 
[[File:pente_canal.jpeg|400px]] <br />
 
Revenons à la modélisation de la houle par l'équation de Berkhoff:
 
Revenons à la modélisation de la houle par l'équation de Berkhoff:
Ligne 77 : Ligne 106 :
 
Le terme CCg n'est pas constant et dépend à présent de la pente du fond du canal :<math>CCg=g.H(x)</math>, avec <math>H(x)=H_0-sx</math><br />
 
Le terme CCg n'est pas constant et dépend à présent de la pente du fond du canal :<math>CCg=g.H(x)</math>, avec <math>H(x)=H_0-sx</math><br />
 
On obtient alors l'équation suivante: <br />
 
On obtient alors l'équation suivante: <br />
<math>H(x)\frac{\partial^2\Phi}{\partial x^2}-s\frac{\partial \Phi}{\partial x}+k^2H(x)\Phi</math> <math>(1)</math>
+
<math>H(x)\frac{\partial^2\Phi}{\partial x^2}-s\frac{\partial \Phi}{\partial x}+k^2H(x)\Phi=0</math> <math>(1)</math>
 
<br />
 
<br />
 
Deux cas à examiner: k constant / k non constant.<br />
 
Deux cas à examiner: k constant / k non constant.<br />
*Si k est constant (<math>k=k_0</math>), <math>(1)</math> devient <math>H(x)\frac{\partial^2\Phi}{\partial x^2}-s\frac{\partial \Phi}{\partial x}+k_0^2H(x)\Phi</math><br />
+
*'''Si k est constant''' (<math>k=k_0</math>), <math>(1)</math> devient <math>H(x)\frac{\partial^2\Phi}{\partial x^2}-s\frac{\partial \Phi}{\partial x}+k_0^2H(x)\Phi=0</math><br />
Par changement de variable (<math>z=H(x)=H_0-sx</math>), l'équation s'exprime finalement <math>z\frac{\partial^2\Phi}{\partial z^2}-s\frac{\partial \Phi}{\partial z}+\frac{k_0}{s}^2z\Phi</math><br />
+
Par changement de variable (<math>z=H(x)=H_0-sx</math>), l'équation s'exprime finalement <math>z\frac{\partial^2\Phi}{\partial z^2}-s\frac{\partial \Phi}{\partial z}+\frac{k_0}{s}^2z\Phi=0</math><br />
En posant <math>\alpha ^2=(\frac{k_0}{s})^2</math>, on obtient une équation de type Bessel: <math>z\frac{\partial^2\Phi}{\partial z^2}-s\frac{\partial \Phi}{\partial z}+\alpha ^2 z\Phi</math><br />
+
En posant <math>\alpha ^2=(\frac{k_0}{s})^2</math>, on obtient une équation de type Bessel: <math>z\frac{\partial^2\Phi}{\partial z^2}-s\frac{\partial \Phi}{\partial z}+\alpha ^2 z\Phi=0</math><br />
*Avec k non constant <math>(k=k_0\sqrt{\frac{H_0}{H(x)}})</math>, (1) devient <math>H(x)\frac{\partial^2\Phi}{\partial x^2}-s\frac{\partial \Phi}{\partial x}+k_0^2H_0\Phi</math> <br />
+
*'''Avec k non constant''' <math>(k=k_0\sqrt{\frac{H_0}{H(x)}})</math>, (1) devient <math>H(x)\frac{\partial^2\Phi}{\partial x^2}-s\frac{\partial \Phi}{\partial x}+k_0^2H_0\Phi=0</math> <br />
 
Alors, en effectuant le changement de variable <math>z=H(x)</math> (1) s'exprime<br />
 
Alors, en effectuant le changement de variable <math>z=H(x)</math> (1) s'exprime<br />
 
<math>z\frac{\partial^2 \Phi}{\partial z^2}+\frac{\partial \Phi}{\partial z}+(\frac{k_0}{s})^2 H_o\Phi=0</math><br />
 
<math>z\frac{\partial^2 \Phi}{\partial z^2}+\frac{\partial \Phi}{\partial z}+(\frac{k_0}{s})^2 H_o\Phi=0</math><br />
En posant <math>\alpha ^2=(\frac{k_0}{s})^2 H_0</math>, on obtient une équation de type Bessel: <math>z\frac{\partial^2\Phi}{\partial z^2}+\frac{\partial \Phi}{\partial z}+\alpha ^2 z\Phi</math><br />
+
En posant <math>\alpha ^2=(\frac{k_0}{s})^2 H_0</math>, on obtient une équation de type Bessel: <math>z\frac{\partial^2\Phi}{\partial z^2}+\frac{\partial \Phi}{\partial z}+\alpha ^2 z\Phi=0</math><br />
 
La forme générale de la solution de ce type d'équation est: <math>\Phi(z)=C_0.J_0(2\alpha\sqrt z)+C_1.Y_0(2\alpha \sqrt z)</math>, où <br />
 
La forme générale de la solution de ce type d'équation est: <math>\Phi(z)=C_0.J_0(2\alpha\sqrt z)+C_1.Y_0(2\alpha \sqrt z)</math>, où <br />
*<math>J_0</math> est une fonction de Bessel de première espèce (définie en 0), et
+
*<math>J_0</math> est une fonction de Bessel de première espèce (définie en 0),
*<math>Y_0</math> une fonction de Bessel de seconde espèce (non définie en 0).<br />
+
*<math>Y_0</math> une fonction de Bessel de seconde espèce (non définie en 0)
 +
* <math>C_0</math> et <math>C_1</math> sont deux constantes à déterminer grâce aux conditions aux limites.<br />
 
Afin de tirer parti des conditions aux limites du problème, exprimons la dérivée de <math>\Phi</math>: <br />
 
Afin de tirer parti des conditions aux limites du problème, exprimons la dérivée de <math>\Phi</math>: <br />
 
<math>\Phi_z(z)=\frac{\alpha}{\sqrt z}C_0.J_1(2\alpha\sqrt z)-\frac{\alpha}{\sqrt z}C_1.Y_1(2\alpha \sqrt z)</math><br />
 
<math>\Phi_z(z)=\frac{\alpha}{\sqrt z}C_0.J_1(2\alpha\sqrt z)-\frac{\alpha}{\sqrt z}C_1.Y_1(2\alpha \sqrt z)</math><br />
Ligne 118 : Ligne 148 :
 
\left\{  
 
\left\{  
 
\begin{array}{l }
 
\begin{array}{l }
C_0=\frac{Y_1^L-Y_0^L}{D}\\
+
C_0=i*\frac{Y_1^L-Y_0^L}{D}\\
 
\\
 
\\
 
C_1=\frac{-(J_1^L-i.J_0^L)}{D}
 
C_1=\frac{-(J_1^L-i.J_0^L)}{D}
Ligne 124 : Ligne 154 :
 
\right.
 
\right.
 
</math>
 
</math>
 +
<br/>
 +
<br/>
 +
<br/>
 +
Finalement, en remplaçant ces expressions de <math>C_0</math> et <math>C_1</math> dans l'expression de <math>\Phi</math>, on obtient:<br/>
 +
<br/>
 +
 +
<math>\Phi(x)=\frac{Y_1^L-Y_0^L}{J_0(Y_1^L-Y_0^L)-Y_0^0(J_1^L-i.J_0^L)}.J_0(2\alpha\sqrt{H_0-s*x} )+\frac{-(J_1^L-i.J_0^L)}{J_0(Y_1^L-Y_0^L)-Y_0^0(J_1^L-i.J_0^L)}.Y_0(2\alpha \sqrt{H_0-s*x})</math>
  
 
===Résolution homotopique===
 
===Résolution homotopique===
Ligne 162 : Ligne 199 :
 
== Cas n°4 : Vague sphérique dans un domaine infini ==
 
== Cas n°4 : Vague sphérique dans un domaine infini ==
 
On étudie dans ce cas une surface dans un domaine de profondeur considérée comme infinie. Cette surface est agité par une onde originaire d'un cercle de rayon <math>r_0</math>. Le domaine est circulaire et l'onde sort librement en <math>r=R</math>. Cela se traduit par les conditions aux limites suivantes : <br />
 
On étudie dans ce cas une surface dans un domaine de profondeur considérée comme infinie. Cette surface est agité par une onde originaire d'un cercle de rayon <math>r_0</math>. Le domaine est circulaire et l'onde sort librement en <math>r=R</math>. Cela se traduit par les conditions aux limites suivantes : <br />
<math>\phi^{r_0} = 1</math> et <math>\phi_r^R = ik\phi</math>
+
[[File:schemaCas4.png|400px]]
 +
 
 +
<math>\phi^{r_0} = 1</math> <br/>
 +
<math>\phi_r^R = ik\phi</math><br/>
 +
En reprenant l'équation de Berhoff, que l'on exprime en coordonnées sphériques on obtient l'équation suivante:
 +
<math>\Phi_{rr}+\frac{1}{r}\Phi_r+k^2\Phi=0</math> avec les conditions aux limites:<br/>
 +
*<math>\Phi^{r_0}=1</math>
 +
*<math>\Phi_r^{R}=i*k*\Phi_R</math>
 +
Il s'agit d'une équation de type Bessel
 +
===Résolution analytique===
 +
La forme générale de la solution de ce type d'équation est:<br/>
 +
<math>\Phi=A.J_0(r)+ B.Y_0(r)</math><br/>
 +
Les conditions aux limites permettent de déterminer A et B:
 +
* <math>\Phi^{r_0}=1 <=> B=\frac{1-A.J_0(r_0)}{Y_0(r_0)}</math>
 +
*<math>\Phi_r^{R}=i*k*\Phi_R <=> A=\frac{Y_1(R)+i*k*Y_0(R)}{J_0(r_0)*Y_1(R)+i*k*J_0(r_0)Y_0(R)-J_1(R)Y_0(r_0)-i*k*J_0(R)Y_0(r_0)}</math>
 +
<br/>
 +
Finalement, <math>\Phi=\frac{Y_1(R)+i*k*Y_0(R)}{J_0(r_0)*Y_1(R)+i*k*J_0(r_0)Y_0(R)-J_1(R)Y_0(r_0)-i*k*J_0(R)Y_0(r_0)}.J_0(r)+ \frac{1-\frac{Y_1(R)+i*k*Y_0(R)}{J_0(r_0)*Y_1(R)+i*k*J_0(r_0)Y_0(R)-J_1(R)Y_0(r_0)-i*k*J_0(R)Y_0(r_0)}.J_0(r_0)}{Y_0(r_0)}.Y_0(r)</math>
 +
 
 +
 
 +
===Résolution par homotopie===
 +
 
 +
*Résolution homotopique
 +
Dans le cas de coordonnées sphérique, l'équation utilisée pour la résolution par homotopie (pour la méthode HAM) s'écrit :
 +
<math>(1−p) \Phi_{rr}+p∗(\Phi_{rr}+\frac{1}{r}\Phi_r+k^2Φ) =0</math>
 +
De la même manière que pour le cas 3, on décompose <math>\Phi</math> et ses dérivées en série entière. En les remplaçant dans l'équation, on aboutit à l'expression des différents ordres d'homotopie :
 +
====Ordre 0====
 +
On a : <math> \phi_{0,rr} = 0 </math> donc <math>\phi_0 = Ar+B</math>
 +
Les conditions aux limites donnent : \phi_0 = \frac {{ik} {1+ik(r0 – R)} r + 1 – frac{{ikr0} {1+ik(r0 – R)}
 +
 
 +
====Ordre n====
 +
On cherche une relation de récurrence entre <math>\Phi_n</math> et <math> \Phi_{n-1}</math>.
 +
L’équation de l’homotopie donne :
 +
<math>\phi_{n,rr} + \frac {1}{r} \phi_{n-1,r} + k^2\phi_{n-1} = 0</math>
 +
On a ainsi : <math>\phi_n = - \iint {1 \over r} \phi_{n-1,r} dr^2 – k^2 \iint \phi_{n-1} dr^2 + A_n r + B_n</math>
 +
On détermine les constantes.
 +
<math>\phi_n^{r0} = ik \phi _n^R</math> donne <math> B_n = (\iint \frac {1} {r} \phi_{n-1,r} dr^2)_{r = r0}
 +
+ k^2 (\iint \phi_{n-1}dr^2)_{r = r0} – A_nr0</math><br/>
 +
 +
<math>\phi_{n,r}^R = ik\phi_n^R </math>donne <br/>
 +
<math>A_n = \frac{1} {1 + ik(r0 – R)} (\int \frac {1} {r} \phi_{n-1,r}dr + k^2\int \phi_{n-1}dr)_{r = R} + ik (\iint  {1 \over r} \phi_{n-1,r}dr^2 + k^2\iint \phi_{n-1}dr^2)_{r = r0} – ik (\iint \frac {1} {r} \phi_{n-1,r}dr^2 + k^2\iint \phi_{n-1}dr^2)_{r = R}</math>
 +
 
 +
Finalement, <math>\phi_n = - \iint {1 \over r} \phi_{n-1,r} dr^2 – k^2 \iint \phi_{n-1} dr^2 + A_n r + B_n</math>
 +
 
 +
 
 +
 
 +
 
 +
Nous n’avons pas pu aboutir à une animation satisfaisante, pour deux raisons, la principale étant probablement un problème d’expression qui nous a échappé. Ci-dessous, le résultat que nous obtenons. On constate que la solution homotopique ne correspond pas à la solution analytique, même à l’entrée du domaine.
 +
 
 +
== Analyse de sensibilité ==
 +
*Cas n°2
 +
{| class="wikitable"
 +
|-
 +
! k=0.1!! k=0.5!! k=1
 +
|-
 +
| [[File:Cas2sens_k01.png|400px]]|| [[File:Cas2sens_kL05.png|400px]]|| [[File:Cas2sens_kL1.png|400px]]
 +
|}
 +
Pour avoir un résultat cohérent il faut obligatoirement avoir kL<1. De plus, plus kL est petit, plus on constate que les solutions analytiques et homotopiques convergent
 +
*Cas n°3
 +
{| class="wikitable"
 +
|-
 +
! k=0.1!! k=0.5!! k=1
 +
|-
 +
| [[File:Cas3sens_k01_Ordre5.png|400px]] || [[File:Cas3sens_k1_Ordre5.png|400px]] ||[[File:Cas3sens_k05_Ordre5.png|400px]]
 +
|}
 +
On arrive donc aux même conclusions que dans le cas 2, à savoir qu’il faut avoir kL<1, et le plus petit possible.

Version actuelle en date du 14 juin 2022 à 15:24

Pour retourner à la page énoncé du sujet : Utilisateur:Jean-Michel Tanguy/SujetENTPE2022

Sommaire

[modifier] Cadre de l'étude

[modifier] Contexte

Le changement climatique entraîne des bouleversements majeurs dans l'ensemble des conditions de vie de notre planète. Parmi eux, nous nous intéressons ici à l'élévation du niveau des océans, qui contribue au recul du trait de côte sur les littoraux. Ce phénomène est bien visible en France, comme l'illustre cette photographie:
érosion du trait de côte, Soulac-sur-Mer, Aquitaine, 3 février 2014
Source:"Sud Ouest" / Julien Lestage
Cette érosion est notamment provoquée par la houle, mouvement ondulatoire de la surface de la mer provoqué par le vent loin des côtes. De nombreuses études sont menées sur ce phénomène dans des laboratoires de recherche et utilisent pour cela des canaux ou bassins à houle (comme au laboratoire M2C-CNRS/Univ. de Rouen/Univ. de Caen https://m2c.cnrs.fr/poles-de-competences/mesocosmes-et-modelisation-analogique/). On peut en voir un exemple ici:

Canal à houle de l'ESITC

Source: Ouest-France/ESITC/EDINBURGH DESIGNS LTD - 15/01/2018

Le schéma ci-dessous permet de comprendre les principales étapes de formation de ce phénomène:
Schéma des phases de la houle
Source: site internet magicsurfschool.com/blog/formation-des-vagues/
La démarche présentée par la suite repose sur la modélisation de la houle par deux méthodes de résolution différentes:

  • une méthode analytique, notamment grâce au logiciel MAXIMA pour les équations les moins évidentes à résoudre analytiquement
  • une méthode d'homotopie, utilisant la décomposition en série entière de la solution

[modifier] Modélisation

On cherche à modéliser une houle dans quatre cas différents. Pour cela, on utilise le modèle de Berkhoff, dont l'équation est :

$ \nabla.(CC_g\nabla \phi)+k^2CC_g\phi=0 $

$ \phi $ est le potentiel, $ k $ est le nombre d’onde, $ C $ est la célérité de l’onde et $ C_g $ est la célérité de groupe des vagues.
On se place dans le domaine des ondes longues, ce qui entraîne que $ C = C_g = \sqrt{gH} $, avec g l'accélération de la pesanteur et H la profondeur du canal. L'équation de Berkhoff devient ainsi :

$ \nabla^2 \phi + k^2\phi = 0 $.

L'évolution de la hauteur de la houle en fonction du temps est donnée par : $ h(x,t)=\Re \left (\phi e^{-i\omega t} \right ) $

Nous allons résoudre cette équation en utilisant la méthode d'homotopie, plus précisément la méthode HAM "Homotopie Perturbation Method". Le principe est de pouvoir passer, par une transformation continue, d'un milieu complexe dont il est difficile de connaitre le comportement à un milieu plus simple.
Elle consiste à résoudre l'équation $ (1-p)L(\phi-\psi)+pH(p)({d^2 \phi \over {dx^2}} + k^2 \phi)=0 $, avec L un opérateur linéaire (ici on prendra la dérivée seconde par rapport à x) et $ \psi $ la solution initiale. p est un paramètre qui varie entre 0 et 1, il représente la progression de la solution initiale ($ p = 0 $) à la solution finale ($ p = 1 $). La méthode HAM implique que $ H(p) = 1 $

[modifier] Cas n°1 : Canal plat ouvert

On considère un canal à fond plat, ouvert aux deux extrémités, de longueur L et de profondeur H. On adopte la notation $ \phi_x $ pour la dérivée de $ \phi $ par rapport à x et $ \phi^x $ pour la valeur de $ \phi $ en x.
En aval, il entre une onde de fréquence unitaire $ \phi=1 $ (condition de Dirichlet) et sortie libre amont $ \phi_x=ik\phi $ (condition de Robin).
Diagramme canal cas 1.png

[modifier] Résolution analytique

Soit à résoudre : $ \phi_{xx} + k^2\phi = 0 $
La solution est donc de la forme : $ \phi = Ae^{-ikx} + B e^{ikx} $
Les conditions aux limites de notre problème sont : $ \phi^0 = 1 $ et $ \phi_x^L = \phi^L $
Donc $ B = 1 $ et $ A = 0 $
Ainsi, $ \phi = e^{ikx} $. L'enjeu est maintenant de retrouver cette solution en utilisant la méthode d'homotopie.


[modifier] Résolution homotopique

Avec une solution initiale nulle, la relation d'homotopie devient :

$ (1-p)\phi_xx + k^2\phi = 0 $
On décompose $ \phi $ en somme de fonctions : $ \phi = \sum_{i=0}^{+\infty} p^i\phi_i $ Cette relation est vrai quel que soit p, donc on peut établir que les coefficients de chaque puissance de p sont nuls aussi. On décompose par ordre, c'est-à-dire par puissance de p.

  • Ordre 0

$ \phi_{0,xx} = 0 $ d'où $ \phi_0 = Ax+B $
On sait que :
$ \phi_0(x = 0) = 1 $ donc $ B = 1 $
et $ \phi_{0,x}(x = L) = ik\phi_0(x =L) $ donc $ A = ik(AL + 1) $ donc $ A = {ik \over{1-ikL}} $
Ainsi, $ \phi_0(x) = {ik\over{1-ikL}} x + 1 $

  • Ordre 1

$ \phi_{1,xx} + k^2\phi_0 = 0 \iff \phi_1 = -k^2\iint \phi_0 dx^2 + Ax + B $
Nos conditions aux limites sont toujours : $ \phi_0^1=0 \ et \ \phi^L_{1,x}=ik\phi^L_1 $. Il vient : $ \phi_1= −k^2L{(k^2L^2+3ikL−3) \over {3(1−ikL)^2}}x−k^2({ik\over{6(1−ikL)}}x^3+ \dfrac {1} {2}x^2) $

  • Ordre n

On a : $ \forall n, \phi_{n,xx} + k^2\phi_{n-1} = 0 $
Ainsi, $ \phi_n = -k^2\iint \phi_{n-1} dx^2 +A_nx +B $. (B = 0 )
Comme $ \phi_{n,x}^L = ik\phi_{n-1}^L, A_n = {1\over {1-ikL}}(k^2\int_0^L\phi_{n-1}dx^2 - ik^3\iint_0^L\phi_{n-1}dx^2) $
Finalement, on obtient :
$ h(x,t) = \sum p^n\Re(\phi_ne^{-i\omega t}) $


On peut visualiser l'évolution de la solution homotopique selon la valeur de p grâce au logiciel Maxima:
GIF Cas 1 final (1).gif


Sur le graphique ci-dessus, la solution analytique est en rouge, la solution homotopique est représentée par la courbe verte. On illustre ici la variation du paramètre p, de 0 (solution en bleue) à 1 (où elle est la plus proche de la solution analytique). La méthode par homotopie permet bien de retrouver la solution analytique.

[modifier] Cas n°2 : Canal plat fermé

On étudie un canal similaire à celui du cas précédent, mais fermé en x=L. Cela se traduit par les conditions aux limites suivantes :
$ \phi_x=ik(2-\phi) $ en x=0 et $ \phi_x = 0 $ en x=L

[modifier] Solution analytique

On cherche $ \phi $ sous la forme $ \phi = Ae^{-ikx}+Be^{ikx} $, car il y a réflexion totale.
Les conditions aux limites donnent respectivement B=1 et $ A = e^{2ikL} $
On a ainsi $ \phi = e^{ik(2L-x)}+e^{ikx} $

[modifier] Solution par homotopie

L’équation homotopique à résoudre s’écrit : $ (1−p)\phi_{xx}+p(\phi_{xx}+k^2\phi)=0 $. La résolution est très similaire au cas 1, seules les conditions aux limites changent.

[modifier] Ordre 0

Elle devient donc $ \phi_{0,xx}=0 \ soit \ \phi_0=Ax+B. $ A l’aide des conditions aux limites, on obtient $ \phi_0 = 2. $

[modifier] Ordre n

A chaque ordre, on a $ \phi_n(x)=−k^2\int \phi_{n-1}(x)dx+A_nx+B_n $ On a les conditions aux limites suivantes : $ \phi_{n,x}^0=−ik\phi_{n,x}^0 \ et \ \phi_{n,x}^L=0 $ Il advient alors $ A_n=(k^2 \iint \phi_{n-1}(x)dxdx)_{x=L} et B_n={i \over k}A_n. $ On programme ensuite la résolution par récurrence sur Maxima.

[modifier] Cas n°3 : Canal en pente

[modifier] Résolution analytique

On s'intéresse à présent à un canal dont le fond est en pente. Le phénomène étudié est toujours unidimensionnel, avec les mêmes conditions aux limites que pour le cas n°1, soit:

  • En aval, onde de fréquence unitaire $ \phi=1 $ (condition de Dirichlet) et
  • En sortie libre amont $ \phi_x=ik\phi $ (condition de Robin).

Pente canal.jpeg
Revenons à la modélisation de la houle par l'équation de Berkhoff: $ \nabla (CCg\Phi) + k^2.CCg\Phi=0 $
Le terme CCg n'est pas constant et dépend à présent de la pente du fond du canal :$ CCg=g.H(x) $, avec $ H(x)=H_0-sx $
On obtient alors l'équation suivante:
$ H(x)\frac{\partial^2\Phi}{\partial x^2}-s\frac{\partial \Phi}{\partial x}+k^2H(x)\Phi=0 $ $ (1) $
Deux cas à examiner: k constant / k non constant.

  • Si k est constant ($ k=k_0 $), $ (1) $ devient $ H(x)\frac{\partial^2\Phi}{\partial x^2}-s\frac{\partial \Phi}{\partial x}+k_0^2H(x)\Phi=0 $

Par changement de variable ($ z=H(x)=H_0-sx $), l'équation s'exprime finalement $ z\frac{\partial^2\Phi}{\partial z^2}-s\frac{\partial \Phi}{\partial z}+\frac{k_0}{s}^2z\Phi=0 $
En posant $ \alpha ^2=(\frac{k_0}{s})^2 $, on obtient une équation de type Bessel: $ z\frac{\partial^2\Phi}{\partial z^2}-s\frac{\partial \Phi}{\partial z}+\alpha ^2 z\Phi=0 $

  • Avec k non constant $ (k=k_0\sqrt{\frac{H_0}{H(x)}}) $, (1) devient $ H(x)\frac{\partial^2\Phi}{\partial x^2}-s\frac{\partial \Phi}{\partial x}+k_0^2H_0\Phi=0 $

Alors, en effectuant le changement de variable $ z=H(x) $ (1) s'exprime
$ z\frac{\partial^2 \Phi}{\partial z^2}+\frac{\partial \Phi}{\partial z}+(\frac{k_0}{s})^2 H_o\Phi=0 $
En posant $ \alpha ^2=(\frac{k_0}{s})^2 H_0 $, on obtient une équation de type Bessel: $ z\frac{\partial^2\Phi}{\partial z^2}+\frac{\partial \Phi}{\partial z}+\alpha ^2 z\Phi=0 $
La forme générale de la solution de ce type d'équation est: $ \Phi(z)=C_0.J_0(2\alpha\sqrt z)+C_1.Y_0(2\alpha \sqrt z) $, où

  • $ J_0 $ est une fonction de Bessel de première espèce (définie en 0),
  • $ Y_0 $ une fonction de Bessel de seconde espèce (non définie en 0)
  • $ C_0 $ et $ C_1 $ sont deux constantes à déterminer grâce aux conditions aux limites.

Afin de tirer parti des conditions aux limites du problème, exprimons la dérivée de $ \Phi $:
$ \Phi_z(z)=\frac{\alpha}{\sqrt z}C_0.J_1(2\alpha\sqrt z)-\frac{\alpha}{\sqrt z}C_1.Y_1(2\alpha \sqrt z) $
Les conditions aux limites de ce cas sont identique à celles du cas 1:

$ \left\{  \begin{array}{l} \Phi(x=0)=1 \\ \Phi_x^L=i.k_L\Phi^L \end{array} \right. $

Comme $ z=H(x)=H_0-sx $, on a $ -s\Phi=i.k_L\Phi $ et $ \frac{\partial \Phi}{\partial x}=-s\Phi_z=i.k_L\Phi $
Alors $ C_0 $ et $ C_1 $ sont solution du système:
$ \left\{ \begin{array}{l l} C_0.J_0(2\alpha\sqrt {H_0})+C_1Y_0(2\alpha\sqrt {H_0})=1 & (condition\ en\ x=0)\\ C_0.J_1(2\alpha\sqrt {H_L})+C_1Y_1(2\alpha\sqrt {H_L})=i.C_0.J_0(2\alpha\sqrt {H_L})+C_1Y_0(2\alpha\sqrt {H_L}) & (condition\ en\ x=L) \end{array} \right. $
L'écriture de ce système se simplifie ainsi:
$ \left\{ \begin{array}{l l} C_0.J_0^0+C_1Y_0^0=1 & (condition\ en\ x=0)\\ C_0.J_1^L+C_1Y_1^L=i.C_0.J_0^L+C_1Y_0^L & (condition\ en\ x=L) \end{array} \right. $
En posant $ D=J_0(Y_1^L-Y_0^L)-Y_0^0(J_1^L-i.J_0^L) $, on obtient les valeurs de $ C_0 $ et $ C_1 $ suivantes:
$ \left\{ \begin{array}{l } C_0=i*\frac{Y_1^L-Y_0^L}{D}\\ \\ C_1=\frac{-(J_1^L-i.J_0^L)}{D} \end{array} \right. $


Finalement, en remplaçant ces expressions de $ C_0 $ et $ C_1 $ dans l'expression de $ \Phi $, on obtient:

$ \Phi(x)=\frac{Y_1^L-Y_0^L}{J_0(Y_1^L-Y_0^L)-Y_0^0(J_1^L-i.J_0^L)}.J_0(2\alpha\sqrt{H_0-s*x} )+\frac{-(J_1^L-i.J_0^L)}{J_0(Y_1^L-Y_0^L)-Y_0^0(J_1^L-i.J_0^L)}.Y_0(2\alpha \sqrt{H_0-s*x}) $

[modifier] Résolution homotopique

La formulation générale de l'homotopie s'écrit $ (1-p)L(\Phi)+p[Equation]=0 $
Dans le cas considéré, en choisissant la dérivée seconde comme opérateur linéaire, cela se traduit par:
$ (1-p)\Phi_{xx}+p[(1-\epsilon x)\Phi_{xx}-\epsilon\Phi_x+k^2\Phi]=0 $$ \epsilon=\frac{s}{H_0} $
En écrivant la décomposition de $ \Phi $:
$ \Phi=\Phi_0+p.\Phi_1+p^2\Phi_2+...+p^n\Phi_n+... $
...et ses dérivées successives:
$ \Phi_{,z}=\Phi_{0,z}+p.\Phi_{1,z}+p^2\Phi_{2,z}+...+p^n\Phi_{n,z}+... $
...que l'on remplace dans l'équation, on obtient l'expression suivante:

$ \begin{array}{l l} &\Phi_{0,x}.p^0\\ + &(\Phi_{1,xx} -\Phi_{0,xx}+ (1-\epsilon x)\Phi_{0,xx} -\epsilon\Phi_{0,x}+k_0^2 ).p^1\\ +&...\\ +&(\Phi_{n+1,xx} -\Phi_{n,xx}+ (1-\epsilon x)\Phi_{n,xx} -\epsilon\Phi_{n,x}+k_0^2 ).p^{n+1} \\ =0 & \end{array} $

Tous les coefficients devant les puissances de p devant être nulles, on en déduit les différents ordres d'homotopie:

[modifier] Ordre 0

$ \Phi_{0,xx}=0 $ donc $ \Phi=A_0.x+B_0 $
Puisque l'on a les mêmes conditions aux limites que dans le cas 1, on peut écrire:

  • En x=0 $ \Phi_0(x=0)=1 <=> B=1 $
  • En x=L $ \Phi_0,x(x=L)=ik(x)\Phi_0(x=L) <=> A=\frac{ik(x)}{1–ik(x)L} $



Finalement, $ \Phi=\frac{ik(x)}{1–ik(x)L}*x+1 $

[modifier] Ordre 1

$ \Phi_{1,xx}-\epsilon(\Phi_{0,xx}+\Phi_{0,x}+k_0^2\Phi_0)=0 $ Donc $ \Phi_1= \int \int ( \epsilon x\Phi_{0,xx}+\epsilon \Phi_{0,x}-k_0^2\Phi_0) \,dx\,dx +A_1x+B_1 $

[modifier] Ordre n

On généralise l'expression précédente et l'on obtient $ \Phi_{n+1,xx}-\epsilon(\Phi_{n,xx}+\Phi_{n,x}+k_0^2\Phi_n)=0 $ $ \Phi_{n+1}= \int \int ( \epsilon x\Phi_{n,xx}+\epsilon \Phi_{n,x}-k_0^2\Phi_n) \,dx\,dx +A_n.x+B_n $

[modifier] Cas n°4 : Vague sphérique dans un domaine infini

On étudie dans ce cas une surface dans un domaine de profondeur considérée comme infinie. Cette surface est agité par une onde originaire d'un cercle de rayon $ r_0 $. Le domaine est circulaire et l'onde sort librement en $ r=R $. Cela se traduit par les conditions aux limites suivantes :
SchemaCas4.png

$ \phi^{r_0} = 1 $
$ \phi_r^R = ik\phi $
En reprenant l'équation de Berhoff, que l'on exprime en coordonnées sphériques on obtient l'équation suivante: $ \Phi_{rr}+\frac{1}{r}\Phi_r+k^2\Phi=0 $ avec les conditions aux limites:

  • $ \Phi^{r_0}=1 $
  • $ \Phi_r^{R}=i*k*\Phi_R $

Il s'agit d'une équation de type Bessel

[modifier] Résolution analytique

La forme générale de la solution de ce type d'équation est:
$ \Phi=A.J_0(r)+ B.Y_0(r) $
Les conditions aux limites permettent de déterminer A et B:

  • $ \Phi^{r_0}=1 <=> B=\frac{1-A.J_0(r_0)}{Y_0(r_0)} $
  • $ \Phi_r^{R}=i*k*\Phi_R <=> A=\frac{Y_1(R)+i*k*Y_0(R)}{J_0(r_0)*Y_1(R)+i*k*J_0(r_0)Y_0(R)-J_1(R)Y_0(r_0)-i*k*J_0(R)Y_0(r_0)} $


Finalement, $ \Phi=\frac{Y_1(R)+i*k*Y_0(R)}{J_0(r_0)*Y_1(R)+i*k*J_0(r_0)Y_0(R)-J_1(R)Y_0(r_0)-i*k*J_0(R)Y_0(r_0)}.J_0(r)+ \frac{1-\frac{Y_1(R)+i*k*Y_0(R)}{J_0(r_0)*Y_1(R)+i*k*J_0(r_0)Y_0(R)-J_1(R)Y_0(r_0)-i*k*J_0(R)Y_0(r_0)}.J_0(r_0)}{Y_0(r_0)}.Y_0(r) $


[modifier] Résolution par homotopie

  • Résolution homotopique

Dans le cas de coordonnées sphérique, l'équation utilisée pour la résolution par homotopie (pour la méthode HAM) s'écrit : $ (1−p) \Phi_{rr}+p∗(\Phi_{rr}+\frac{1}{r}\Phi_r+k^2Φ) =0 $ De la même manière que pour le cas 3, on décompose $ \Phi $ et ses dérivées en série entière. En les remplaçant dans l'équation, on aboutit à l'expression des différents ordres d'homotopie :

[modifier] Ordre 0

On a : $ \phi_{0,rr} = 0 $ donc $ \phi_0 = Ar+B $ Les conditions aux limites donnent : \phi_0 = \frac {{ik} {1+ik(r0 – R)} r + 1 – frac{{ikr0} {1+ik(r0 – R)}

[modifier] Ordre n

On cherche une relation de récurrence entre $ \Phi_n $ et $ \Phi_{n-1} $. L’équation de l’homotopie donne : $ \phi_{n,rr} + \frac {1}{r} \phi_{n-1,r} + k^2\phi_{n-1} = 0 $ On a ainsi : $ \phi_n = - \iint {1 \over r} \phi_{n-1,r} dr^2 – k^2 \iint \phi_{n-1} dr^2 + A_n r + B_n $ On détermine les constantes. $ \phi_n^{r0} = ik \phi _n^R $ donne $ B_n = (\iint \frac {1} {r} \phi_{n-1,r} dr^2)_{r = r0} + k^2 (\iint \phi_{n-1}dr^2)_{r = r0} – A_nr0 $

$ \phi_{n,r}^R = ik\phi_n^R $donne
$ A_n = \frac{1} {1 + ik(r0 – R)} (\int \frac {1} {r} \phi_{n-1,r}dr + k^2\int \phi_{n-1}dr)_{r = R} + ik (\iint {1 \over r} \phi_{n-1,r}dr^2 + k^2\iint \phi_{n-1}dr^2)_{r = r0} – ik (\iint \frac {1} {r} \phi_{n-1,r}dr^2 + k^2\iint \phi_{n-1}dr^2)_{r = R} $

Finalement, $ \phi_n = - \iint {1 \over r} \phi_{n-1,r} dr^2 – k^2 \iint \phi_{n-1} dr^2 + A_n r + B_n $



Nous n’avons pas pu aboutir à une animation satisfaisante, pour deux raisons, la principale étant probablement un problème d’expression qui nous a échappé. Ci-dessous, le résultat que nous obtenons. On constate que la solution homotopique ne correspond pas à la solution analytique, même à l’entrée du domaine.

[modifier] Analyse de sensibilité

  • Cas n°2
k=0.1 k=0.5 k=1
Cas2sens k01.png Cas2sens kL05.png Cas2sens kL1.png

Pour avoir un résultat cohérent il faut obligatoirement avoir kL<1. De plus, plus kL est petit, plus on constate que les solutions analytiques et homotopiques convergent

  • Cas n°3
k=0.1 k=0.5 k=1
Cas3sens k01 Ordre5.png Cas3sens k1 Ordre5.png Cas3sens k05 Ordre5.png

On arrive donc aux même conclusions que dans le cas 2, à savoir qu’il faut avoir kL<1, et le plus petit possible.

Outils personnels