Transcript corrigé - IC Lycée Paul Constans

PT – IC
27 novembre 2014
C ORRIGÉ DU D EVOIR N°1
Durée : 2 heures. Tout document interdit, y compris la calculatrice.
D’après un sujet 0 des Mines.
Modélisation numérique des milieux granulaires.
I. Présentation des milieux granulaires
Introduction. L’étude des écoulements granulaires est un domaine scientifique en pleine expansion. Assez
solides pour soutenir le poids d’un immeuble, ils peuvent couler comme de l’eau dans un sablier ou être
transportés par le vent pour sculpter les dunes et les déserts.
La nature offre ainsi les exemples les plus spectaculaires de phénomènes et de structures où interviennent les
milieux granulaires : dunes de sable, plages s’étirant le long des côtes, éboulis, avalanches de neige, figures
d’érosion. . .
II. Calcul de la forme d’un tas par automates cellulaires
L’objectif de cette partie est la simulation numérique d’un tas de grain par un modèle appelé automate cellulaire. Cet automate se présente sous la forme d’un quadrillage 2D dont chaque case peut être occupée ou vide.
Un grain y est symbolisé par une case occupée. La configuration des cases, qu’on appelle état de l’automate,
évolue au cours du temps selon certaines règles très simples. Dans ce qui suit, nous allons simuler la formation d’un demi-tas de sable, situé à droite de l’axe de symétrie vertical. Le tas complet s’obtient en faisant la
symétrie par rapport à l’axe vertical.
Le problème. Une tour de hauteur H (voir figure 1) est une colonne de H cases pleines consécutives, et
symbolise un empilement de H grains à la même abscisse. On note h le nombre de cases de la tour dont les
voisines de droite sont vides. Si h > 0, cela signifie que la tour de droite a une hauteur de H − h. Si h > 1, on
détermine arbitrairement un nombre n de grains du sommet de la tour qui vont tomber sur celle de droite
avec l’instruction python :
n = floor((h + 2.0) / 2 * random()) + 1
où floor est la partie entière importée du module math et random() est une fonction (ne prenant pas d’argument) importée du module random qui renvoie de manière aléatoire un flottant de l’intervalle semi-ouvert
[0; 1[. Dans l’exemple figure 1, lors du passage de l’état 1 à l’état 2, on a n = h = 2. Lors du passage de l’état 2
à l’état 3, h = 3 et n = 1. L’état 3 est stable : plus aucun grain ne tombe car la différence de taille h entre une
colonne et sa voisine de droite n’excède jamais 1.
État 1
État 2
État 3
Grains qui vont tomber
Grains pouvant tomber mais restant immobiles
Grains stables
F IGURE 1 – Exemple d’évolution de l’automate cellulaire
page 1/ 7
PT – IC
27 novembre 2014
C ORRIGÉ DU D EVOIR N°1
1. Quel est le type de n ? Déterminer un encadrement de n pour h > 1.
n est un entier. De plus, pour tout x ∈ R, on a x − 1 < bxc É x d’où
h +2
h +2
random() < n É
random() + 1
2
2
h +2
h
+1 ⇒ 1 É n É +1 .
2
2
En particulier, pour h > 1, n É h (heureusement !)
Or random() ∈ [0; 1[ donc 1 É n <
2. Écrire une fonction nommée calcul_n(h) qui prend h comme argument et qui renvoie le nombre n de
grains qui vont tomber sur la colonne suivante.
def calcul_n(h):
"""Renvoie le nombre de grains qui vont tomber sur la colonne suivante."""
if h > 1:
return floor((h + 2.0) / 2 * random()) + 1
else:
return 0
Pour représenter l’automate en Python, il est inutile de considérer une représentation à 2 dimensions : stocker
la hauteur H de chaque colonne suffit. Le demi-tas est complètement défini par un tableau stockant la hauteur
de chaque colonne. Le nombre de colonnes est fixé au début de la simulation. L’évolution du demi-tas que
l’on propose est la suivante :
• au départ chaque colonne est vide (et a donc pour hauteur 0) ;
• périodiquement, on dépose un grain sur la première colonne (de gauche) ;
• une étape d’actualisation du demi-tas consiste à faire évoluer le nombre de grains comme dans la figure 1 :
on parcourt le demi-tas de gauche à droite en déplaçant un certain nombre de grains (éventuellement
zéro) d’une colonne vers sa voisine de droite suivant les règles définies au dessus. À l’extremité droite du
demi-tas il n’y a rien : les grains qui tombent sont perdus.
Rappel : Lorsque l’on passe un tableau en paramètre d’une fonction à Python, on passe uniquement l’emplacement mémoire où se trouve le tableau. Modifier une entrée dans la fonction modifie donc « globalement »
le tableau. Par exemple :
>>>
...
...
>>>
>>>
>>>
[3,
def f(T):
T[0]=3
tableau = [8, 9, 10]
f(tableau) # rien ne s'affiche car f ne renvoie rien
print(tableau)
9, 10]
3. Définir une fonction initialisation(nbCols) prenant en entrée un entier nbCols > 0 et renvoyant un
tableau représentant le demi-tas au début de la simulation, c’est à dire un tableau de taille nbCols rempli
de zéros.
def initialisation(nbCols):
"""Renvoie un tableau représentant le demi-tas au début de la simulation."""
return [0] * nbCols
4. Écrire la fonction d’actualisation actualise(T) qui prend en entrée un tableau représentant un demi-tas
et le fait évoluer en utilisant les règles définies précédemment. Cette fonction doit renvoyer le nombre de
grains perdus (tombés de la dernière colonne) lors de cette étape d’actualisation.
page 2/ 7
PT – IC
C ORRIGÉ DU D EVOIR N°1
27 novembre 2014
def actualise(T):
"""actualise le demi-tas T et renvoie le nombre de grains perdus."""
for col in range(len(T)-1): # Le comportement de la dernière est différent.
h = T[col] - T[col+1]
n = calcul_n(h)
T[col] -= n
T[col+1] += n
h = T[-1] # La colonne suivante est vide.
n = calcul_n(h)
T[-1] -= n
return n
5. Écrire le bloc d’instructions du programme principal
nbCols = int(input("Entrez le nombre de colonnes : "))
demi_tas = initialisation(nbCols)
perdu = 0
nb_tour = 0
while perdu < 1000:
if nb_tour % 10 == 0:
demi_tas[0] += 1
perdu += actualise(demi_tas)
nb_tour += 1
On rappelle que pour tracer un graphique en utilisant le module
matplotlib.pyplot, il suffit de donner deux tableaux X et Y de
tailles égales représentant les abscisses et les ordonnées, suivant la
syntaxe suivante :
import matplotlib.pyplot as plt
plt.plot(X,Y)
6. Écrire une fonction dessine(T) prenant en paramètre un demitas et qui trace l’allure du tas comme dans la figure 2. Celui-ci
est obtenu en faisant une symétrie du demi-tas par rapport à
F IGURE 2 – Un tas avec nbCols = 100
sa première colonne (Si nbCols est le nombre de colonnes du
demi-tas, on obtient donc un tas à 2nbCols − 1 colonnes). On
pourra utiliser list(range(a,b)) pour obtenir un tableau contenant les entiers de a à b − 1.
def dessine(T):
"""Dessine l'allure du tas correspondant au demi-tas T."""
nbCols = len(T)
x = list(range(-nbCols+1, nbCols))
y = T[-1:0:-1] + T # Ici, le slicing est très efficace !
plt.plot(x, y)
plt.show()
page 3/ 7
PT – IC
27 novembre 2014
C ORRIGÉ DU D EVOIR N°1
III. Méthode des éléments distincts (DEM)
Calcul de
Pour remédier à certaines limitations des modèles par autol’accélération
mates cellulaires, il est possible de modéliser les milieux granulaires par la méthode des éléments distincts DEM issue
Calcul des
des modèles moléculaires et présentée par Cundall en 1971.
forces
Les grains sont considérés comme des corps isolés et les inIntégrations
teractions comme des forces extérieures. Le mouvement de
des équation
chaque grain est donné par les équations de la dynamique.
du mouvement
La boucle principale simplifiée de l’algorithme DEM est pré- Incrément du
sentée ci-contre. Nous limitons l’étude à des éléments 2D temps t = t + dt
(disques) et tous les grains sont supposés semblables : ils ont
Sauvegarde
donc tous les mêmes caractéristiques.
de l’état
Un des principaux inconvénients de la méthode est la détermination des paramètres physiques qui permettront une
F IGURE 3 – Principe général du DEM
bonne corrélation entre la simulation et l’expérimentation.
Différents laboratoires de recherche ont créé une base de données commune qui leur permet de partager
des jeux de paramètres et de les noter de 1 à 5 selon les résultats obtenus lors des simulations. Une fonction
d’agrégation permet de calculer la note moyenne. Il est donc possible d’effectuer des requêtes de recherche
dans la relation « granular_base » définie par la table ci-après dont les lignes sont données en exemple :
Labo
MSC
LPMDI
LPGP
LMGC
PMMH
..
.
Geom
poly
poly
sphere
sphere
quartz
..
.
Mat
Densite
verre
acier
metal
verre
verre
..
.
2.1
6
4.7
1.56
1.46
..
.
granular_base
R
Kn
Gamma
4
1
2
..
.
7
3.10
5.106
106
2.107
7.108
..
.
20
1
10
5
0.5
..
.
Kt
6
2.10
4.105
3.104
107
2.108
..
.
Mu
Note
Votant
...
0.5
0.7
0.2
0.3
0.4
..
.
2
3
3
4
2
..
.
2
6
3
5
1
..
.
...
...
...
...
...
..
.
F IGURE 4 – Extrait de la table « granular_base »
7. Écrire la requête en langage SQL qui va permettre de trouver les valeurs de paramètres possibles (Densite, R, Kn, Gamma, Kt, et Mu) pour des sphères en verre avec une note strictement supérieure à 3 sur un
échantillon d’au moins 3 votants.
SELECT Densite, R, Kn, Gamma, Kt, Mu
FROM granular_base
WHERE Geom='sphere' AND Mat='verre' AND Note>3 AND Votant>=3 ;
Après une phase d’initialisation qui permet de récupérer les différentes valeurs issues de la base de données,
on exécute un code Python (hors programme et non reproduit) qui va permettre une manipulation agréable
des grains. On dispose ainsi de la faculté de manipuler des objets de type « grain » aussi simplement que les
types prédéfinis de Python. On crée un nouveau grain g en effectuant g=grain(x, y, m) où x et y désignent
les coordonnées de la position du centre du grain et m sa masse. Étant donné un grain g, on peut accéder à
sa masse, sa position, sa vitesse et la force qui le concerne par g.masse, g.pos, g.vit et g.force qui sont
chacun un couple de flottants (initialisé à (0,0) pour la vitesse et la force). Par exemple :
grain_i = grain(0, 1, 1)
grain_j = grain(1, 2.7, 1)
(Xi, Yi) = grain_i.pos
grain_j.force = (0, -9.8)
#
#
#
#
Un premier grain de coordonnées (0,1) de masse 1
Un autre grain de coordonnées (1,2.7) de masse 1
On a alors Xi = 0 et Yi = 1
Modifie la force qui s'exerce sur grain_j
Pour commencer la simulation, on crée un tas avec la fonction initDEM ci-dessous, où la fonction sqrt est la
fonction racine carrée importée du module math :
page 4/ 7
PT – IC
C ORRIGÉ DU D EVOIR N°1
27 novembre 2014
def initDEM(R, n, M):
"Renvoie une liste de grains de rayon R, de masse M formant un tas à n niveaux."
tas = [] # Initialisation de la liste
for i in range(n):
for j in range(n-i):
tas.append(grain(0 + i*R + j*2*R, R + i * sqrt(3)*R, M))
return tas
Par exemple, pour obtenir un tas de grains de rayon 1 et de masse 1 (unités S.I.) agencés en 5 niveaux puis
déplacer le premier grain d’une unité de longueur selon la direction #»
x par :
tas = initDEM(1, 5, 1)
(X0, Y0) = tas[0].pos
tas[0].pos = (X0+1, Y0)
8. Justifier et esquisser la forme du « tas » obtenu par initDEM(1, 5, 1).
La « distance entre
p deux couches » est la hauteur d’un triangle équilatéral de côté
2R, c’est à dire R 3 :
π/3
2R
p
R 3
2R
p
R 3
R
9. Les positions sont calculées de façon à ce que théoriquement les grains soient empilés sans aucune interpénétration, mais est-ce possible numériquement ?
Non : les erreurs d’arrondis font qu’il y aura forcément des interpénétration.
Pour modéliser les interactions entre les particules, nous allons utiliser la méthode des solides déformables
« sphères molles » : le calcul des forces et des moments sera réalisé en considérant que les disques sont indéformables mais peuvent s’interpénétrer légèrement avec δN et δT le déplacement normal et tangentiel à partir
du contact. Ainsi, pour modéliser les forces normales, on utilise un modèle de ressort (raideur k N ) associé à
une dissipation visqueuse (coefficient γn ) permettant de reproduire une collision inélastique. Pour les forces
tangentielles, un ressort (raideur k T ) couplé à un patin (limite du glissement µ) permet de modéliser la force
de friction (voir figure 5 ).
F IGURE 5 – Sphères molles avec interpénétration exagérée
page 5/ 7
PT – IC
27 novembre 2014
C ORRIGÉ DU D EVOIR N°1
Dans la suite du problème, nous nous intéressons à la résolution du Principe Fondamental de la Dynamique
en résultante selon #»
x et #»
y pour chaque grain i :
d 2 xi
d 2 yi
=
F
et
m
.
= Fi y − m i .g
(1)
i
x
i
dt2
dt2
On considère qu’on dispose de la fonction calcul_Fnt(gi, gj) qui prend en arguments deux grains g i et
g j et qui renvoie (F j i n , F j i t ) où F j i n et F j i t correspondent respectivement aux composantes algébriques nor#»
males (direction #»
n ) et tangentielles (direction t ) de la force exercée par le grain g j sur le grain g i . Si les deux
grains ne sont pas en contact alors la fonction renvoie (0, 0).
mi .
10. À partir du paramétrage fourni figure 5, écrire une fonction angle(gi, gj) qui prend en paramètre deux
grains g i et g j et qui renvoie la valeur de cos(α) et sin(α) avec α = ( #»
x , #»
n ).
def angle(gi, gj):
"""Renvoie (cos(al), sin(al)) où al est l'angle (x, n)."""
(Xi, Yi) = gi.pos
(Xj, Yj) = gj.pos
(dx, dy) = (Xj - Xi, Yj - Yi) # Coordonnées du vecteur CiCj
n = sqrt(dx ** 2 + dy ** 2)
# sa norme
return (dx / n, dy / n)
11. En déduire une fonction somme_int(tas) prenant en argument un tableau de grains représentant un tas
et qui, pour chaque grain g i , calcule la somme des interactions Fi x et Fi y (selon #»
x et #»
y ) exercées par tous
les autres grains sur g i , et affecte cette force résultante au champ force du grain g i . On utilisera la fonction
calcul_Fnt.
def somme_int(tas):
"""met à jour les champs gi.force de tous les grains du tas."""
for i in range(len(tas)):
gi = tas[i]
(Fx, Fy) = (0, 0)
for j in range(len(tas)):
if j != i:
gj = tas[j]
cos_al, sin_al = angle(gi, gj)
Fjin, Fjit = calcul_Fnt(gi, gj)
Fx += Fjin * cos_al - Fjit * sin_al
Fy += Fjin * sin_al + Fjit * cos_al
gi.force = (Fx, Fy)
Le schéma d’intégration utilisé dans la suite sera celui de l’algorithme de Verlet « saute- mouton » (leapfrog)
qui possède certains avantages par rapport à la méthode d’Euler notamment en termes de réversibilité du
temps, de précision et de stabilité. On calcule les positions des grains aux temps t = 0, ∆t , 2∆t , . . .où ∆t est le
pas de temps.
Les vitesses sont calculées et mémorisées pour des temps intermédiaires, t = ∆t /2, 3∆t /2, . . .On utilise les vitesses intermédiaires pour
∆t
∆t
déterminer les positions aux instants k∆t (voir figure ). On note :
v
v
k−1/2
• x k la position du grain au temps t = k∆t ;
x k−1
• v k+1/2 la vitesse au temps t = (k + 1/2)∆t ;
k+1/2
xk
∆t
• a k l’accélération au temps t = k∆t .
page 6/ 7
x k+1
t
PT – IC
C ORRIGÉ DU D EVOIR N°1
27 novembre 2014
12. Donner une formule de calcul de v k+1/2 en fonction de v k−1/2 , a k et ∆t ainsi que x k+1 en fonction de x k ,
v k+1/2 et ∆t .
On a v k+1/2 = v k−1/2 + a k × ∆t et x k+1 = x k + v k+1/2 × ∆t .
13. Écrire une fonction etape_Verlet(tas, Dt) qui prend en argument un tableau tas représentant un tas
à un certain instant, le pas de temps utilisé Dt et qui met à jour le tableau tas conformément au schéma
d’intégration pour obtenir la situation à l’instant suivant. On utilisera la fonction somme_int définie précédemment. On prendra aussi en compte la pesanteur telle que définie dans l’équation de la dynamique (1).
On rappelle que chaque grain a été créé avec une masse qui peut être lue par tas[i].masse. On pourra
considérer que l’on dispose d’une variable globale g désignant l’accélération de la pesanteur.
def etape_Verlet(tas, Dt):
"""met à jour le tableau tas conformément au schéma d'intégration pour
obtenir la situation à l'instant suivant, le pas de temps étant Dt."""
somme_int(tas)
for grain in tas:
(Fx, Fy) = grain.force
(Vx, Vy) = grain.vit
(x, y) = grain.pos
(ax, ay) = (Fx / grain.masse, -g + Fy / grain.masse)
Vx += ax * Dt ; Vy += ay * Dt
x += Vx * Dt ; y += Vy * Dt
grain.pos = (x, y)
grain.vit = (Vx, Vy)
14. Écrire une fonction simulation qui prend en argument les paramètres permettant la génération d’un
tas initial par initDEM, la durée Tmax pour laquelle on simule l’écoulement du tas et le nombre de pas de
simulation nb_pas et qui exécute la simulation et renvoie le tas obtenu à la fin.
def simulation(R, n, M, Tmax, nb_pas):
"""Renvoie le tas obtenu à la fin de la simulation sur un tas n couches de
grains de rayon R et de masse M, pour un temps de 0 à Tmax en nb_pas."""
tas = initDEM(R, n, M)
Dt = Tmax / nb_pas
t = 0
while t < Tmax:
etape_Verlet(tas, Dt)
t += Dt
return tas
IV. Optimisation de l’algorithme
La simulation des modèles par éléments distincts (DEM) est très gourmande en temps processeur et en allocation mémoire, il est donc nécessaire d’optimiser le code.
15. Quelle est la partie de l’algorithme proposé qui est la plus pénalisante ? Évaluer sa complexité.
C’est clairement la fonction somme_int dont la complexité est en O(n 2 ) où n est le nombre de grain.
16. Proposer (sans coder) des améliorations possibles afin de diminuer cette complexité et donner une estimation du gain ainsi obtenu.
On pourrait stocker le tas par couches, chaque couche pouvant interagir uniquement avec les deux
p
couches adjacentes. Vu la forme du tas de départ, un couche contient O( n) grains, et on obtient donc
p
une complexité en O(n n).
page 7/ 7