TP_Geologie.sagews

TP sur un problème de géologie.

Le but de ce TP est d'aborder un problème de "géologie". Nous allons modéliser la propagation d'un front d'onde dans le sol. La difficulté sera de tenir compte des inhomogénéités du terrain pour décrire la propagation de l'onde. Ce TP se décomposera en $2$ sous parties

1. La définition du milieu,
2. et le calculs effectifs de la propagation de l'onde.

La définition du milieu

Le terrain que nous allons chercher à modéliser sera représentée par une fonction à $2$ coordonnées $x$ et $y$ (la profondeur), qui à chaque position associe une densité. Dans notre modèle, cette densité indiquera la vitesse $v(x,y)$ de propagation des ondes en ce point. On supposera que la vitesse varie entre 0. (aucune propagation) et 1. (vitesse maximale). Par exemple, pour un terrain parfaitement homogène, on définira une fonction comme ci-après :

def milieu_homogene(x,y):
    #milieu parfaitement homogène
    return 1.0

La fonction ci-dessus renvoie pour toutes positions la valeur constante $1$, donc dans tout le milieu l'onde se propage à vitesse constante. Cependant, ce modèle est assez peu réaliste car le sol est généralement structuré par couches. On parle de milieu stratifié.

Écrire une fonction qui définit un milieu contenant $3$ couches dans lesquelles nous avons trois vitesses constantes différentes. On supposera que $y$ varie entre 0 et -1, on découpera l'espace en trois couches égales avec les vitesses : 0.35, 0.75 et 1. La vitesse varie donc uniquement selon la profondeur $y$.

Il est possible d'afficher la valeur de la vitesse dans le milieu à l'aide la fonction density_plot. Si on exécute la fonction sur le milieu homogène on obtient un restangle uni.

density_plot(milieu_homogene,(-5,5),(-1,0))

Utilisez density_plot pour afficher la densité de votre milieu à trois couches.

La fonction density_plot peut prendre un argument cmap qui permet de faire varier les couleurs utilisées. cmap peut prendre les valeurs suivantes.

import matplotlib.cm as cm
for k in cm.datad:
    print k

density_plot(milieu_homogene,(-5,5),(-1,0), cmap = "Blues")

Testez différentes valeurs pour cmap pour afficher votre densité

Une manière de définir un obstacle est de considérer une zone dans laquelle la vitesse de l'onde est nulle. Dans ce cas, l'onde ne se propage plus et est arrêtée par l'obstacle.
Ecrivez deux fonctions obstacle_rectangle et obstacle_rond qui prennent en argument x et y et retournent True si le point est à l'intérieur de l'obstacle et et False sinon. Les obstacles doivent être un rectangle et un rond, la taille et la position sont laissées libres.

Ecrivez quatre fonctions qui correspondent aux quatre milieux suivants :

• milieu homogène avec obsctacle rectangle

• milieu homogène avec obsctale rond

• milieu 3 couches avec obstacle rectangle

• milieu 3 couches avec obstacle rond

La vitesse doit être de 0 à l'intérieur de l'obstacle et suivre la règle habituelle du milieu sinon. Pensez à utiliser les fonctions écrites précédemment !!

Affichez chaque milieu avec density_plot

La propagation d'une ondes

Nous allons maintenant modéliser la propagation d'une onde dans nos différents domaines. Pour simplifier, on supposera que la position du front d'onde est décrite par : $$(x(t),y(t)) = r(t) (cos(\theta),sin(\theta))$$ où l'angle $\theta$ ne dépend pas de $t$ et $r(t)$ est une distance postivie $r(t)>0$. Comme les ondes se propagent dans le sol, donc vers le bas, on supposera $\theta \in [\pi, 2\pi]$ (ainsi, la coordonnée $y$ est toujours négative).

Dans une direction $\theta \in [\pi,2\pi]$ fixée, on suppose que la vitesse $v$ vérifie $$r'(t) = v\left(r(t)cos(\theta), r(t)sin(\theta) \right).$$

Ecrire une fonction vitesse qui prend en paramètre un angle theta, un milieu (sous la forme de fonction python) et un rayon r et qui retourne la vitesse au point donné.

Soit $f(r)$ la vitesse au rayon $r$ pour un angle et un milieu donné. L'équation différentielle $r'(t) = f(r)$ ne peut pas être résolue exactement, car la fonction $f$ est a priori quelconque. Nous allons donc utiliser une méthode numérique de résolution. L'idée est d'utiliser l'algorithme suivant : $$\dfrac{r(t_{n+1}) - r(t_{n})}{\delta_t} = f(r(t_n))$$ où $t_n = n \delta_t$, et $\delta_t>0$ est un paramètre petit. L'idée de la méthode est de remplacer la dérivée $r'(t)$ par un taux d'accroissement. Remarquer que quand $\delta_t \rightarrow 0$, ce taux d'accroissement tend vers la dérivée $r'(t_n)$. Dans l'expression ci-dessus, on peut isoler $r(t_{n+1})$ comme suit : $$r(t_{n+1}) = \delta_t f(r(t_n)) + r(t_{n}).$$ Étant donnée $r(t_0=0)=0$, on peut calculer la suite des $r(t_{n})$ avec l'expression ci-dessus jusque $t_N = T$, où $T$ est un paramètre choisi.

À l'aide de la méthode décrite (méthode d'Euler explicite), écrire une fonction prenant en paramètre le milieu, l'angle theta, le temps $T$ et le pas $\delta_t$ et qui retourne la liste des valeurs $r(t_n)$ jusqu'à $t = T$. (Pour les tests, on choisira $T=1$ et $\delta_t = 0.1$.)

On considère un émetteur placé sur le sol en coordonnée $(0,0)$. Il émet une onde qui va se propager dans toutes les directions $\theta \in [\pi,2\pi]$. Avec ce que nous avons développé jusqu'ici, le but est de représenter l'onde se propageant dans un milieu choisi.

Écrire une fonction qui prend en argument un milieu et un temps $T$ de simulation, et qui renvoie une animation du front d'onde se propageant dans le milieu considéré (tester avec les différents milieux définis).

Indication : Pour calculer la solution dans 'toutes' les directions $\theta$, considérer simplement $50$ valeurs de $\theta$ reparties équitablement dans l'intervalle $[\pi,2\pi]$.

︠544ccc2d-10c6-41c0-8d3f-b2bedb7e24c8︠
︠53170748-686c-40e3-a179-73051262c9d4︠
︠40867a13-1fab-4cc4-a235-c58d4a7c1dd2︠