10  5/ Solides discrets, modèle de treillis

Flambage (instabilité) d’une barre en compression

Objectifs

  • Gérer un système à connectivité fixe pour simuler un solide “macroscopique”
  • Modéliser un treillis linéaire élastique à température nulle
  • Ajuster des paramètres de l’énergie potentielle à une réponse macroscopique de la structure
  • Prendre l’initiative dans l’écriture du code

10.1 Système à connectivité fixe

La connectivité est la liste des particules qui ont une interaction. Dans le TP 1, aucune particule n’avait d’interaction avec d’autre particule. Dans le TP 2 et la première partie du TP 3, chaque particule avait une interaction avec toutes les autres particules du système. Dit autrement, chaque particule était connectée à chaque autre. Dans la seconde partie du TP 3, nous avons introduit un rayon de troncation pour les interactions, ainsi chaque particule n’interagit qu’avec ses proches voisines, elle n’est plus connectée à toutes les particules du système.

En revanche, dans les calculs avec rayon de troncation, la connexion d’une particule avec une autre était susceptible de changer : si deux particules prennent des trajectoires qui les éloignent trop, elles finissent par ne plus être connectées. Inversement, des particules éloignées peuvent être amenées à interagir. Dans un système à connectivité fixée, les connexions entre particules ne changent jamais : elles ne dépendent pas de la position. Autrement dit, les paires de particules restent identiques pendant la simulation, quel que soit l’éloignement des particules.

C’est notamment le cas de structures en treillis, comme celui ci-dessous :

Pont en treillis métallique, auteur: Leonard G., via Wikimedia

La Figure 10.1 montre une vue schématisée du pont à treillis ci-dessous. Ici, les “particules” sont les nœuds (numérotés de 0 à 7) reliés par des barres que l’on peut considérer élastiques linéaires de raideur \(k_0\).

Figure 10.1: Treillis schématisé avec noeuds numérotés

La connectivité des nœuds est une tableau qui indique quel nœud i est connecté à quel autre nœud j :

Connectivité pour le treills de la Figure 10.1
i j
0 1
0 5
1 5
1 2
2 5
2 6
2 7
2 3
3 7
3 4
4 7
5 6
6 7

Si i et j sont des tableaux Numpy, alors on peut calculer la matrice rij comme cela :

rij = positions[:, j] - positions[:, i]

Les étapes de calcul des forces restent inchangées : le calcul de dij, nij, fij est identique à ce que l’on a fait pour les TP précédents. En revanche la somme fi = fij.sum(axis=2) doit être effectuée par la fonction suivante :

def assemble_forces(fij, natoms, i, j):
    """Calcule la somme des forces sur j pour chaque atome i"""
    fi = np.zeros([fij.shape[0], natoms])
    for d in range(fij.shape[0]):
        fi[d] += 0.5 * np.bincount(j, weights=fij[d], minlength=natoms)
        fi[d] -= 0.5 * np.bincount(i, weights=fij[d], minlength=natoms)
    return fi

10.2 Énergie potentielle

Pour une barre entre i et j de longueur initiale \(L_{ij}\) et de raideur \(k_{ij}\), l’énergie potentielle harmonique s’écrit :

\[U_2(d_{ij}) = \frac{1}{2}k_{ij}(d_{ij} - L_{ij})^2\]

Notez que toutes les barres n’ont pas la même longueur initiale \(L_{ij}\), ni forcément la même raideur \(k_{ij}\).

Consigne

Sur le modèle des TP précédents, implémentez les fonctions suivantes :

def elastic_bar_force(dij, kij, Lij):
    return ...

def elastic_bar_energy(dij, kij, Lij):
    return ...

def compute_force_energy_bar(positions, i, j, kij, Lij):
    """Calcule les forces à partir d'une connectivité fixe donnée par i et j"""
    return ..., ...

Réutilisez le plus possible les fonctions développées lors des précédents TP.

10.3 Extension d’une barre en treillis

Figure 10.2: Barre en treillis en traction

L’objectif de cette partie du TP est de simuler l’extension d’une barre encastrée à gauche où l’on vient appliquer une force de traction à droite. Pour cela, on approxime le comportement de la barre par le treillis montré à la Figure 10.2. Comme on contrôle le comportement des barres du treillis mais pas le comportement “macroscopique”, ou global, de la structure, on cherche à accomplir deux tâches :

  • Calculer le comportement effectif de la barre encastrée (c’est-à-dire le module de Young et le coefficient de Poisson) pour une distribution de \(k_{ij}\) donnée,
  • Ajuster la distribution de raideurs \(k_{ij}\) pour obtenir les propriétés effectives désirées.

10.3.1 Positions initiales

On cherche à créer le tableau Numpy des positions pour une grille rectangulaire de 5 lignes fois 20 colonnes (nx = 20, ny = 5).

Consigne

La fonction make_crystal du TP 3 peut être modifiée pour générer la structure de la Figure 10.2. Créez une fonction make_bar à partir du code de make_crystal. Vérifiez que les positions initiales sont correctes avec plt.scatter(*positions).

Calculez l’espacement entre les nœuds pour que la barre fasse 200 cm de long.

10.3.2 Conditions aux limites

On a deux conditions aux limites :

  1. Le déplacement est nul pour les nœuds bleus (à gauche de la barre)
  2. Une force horizontale est appliquée aux nœuds rouges (à droite de la barre)

Pour gérer ces conditions aux limites, il faut pouvoir identifier des groupes de nœuds. Une façon possible de faire est de tester la position en \(x\) des nœuds. Par exemple, pour obtenir les nœuds proches de \(x = 1\), on pourra écrire :

node_group = np.abs(positions[0, :] - 1) < 1e-6

Le tableau node_group est alors un tableau de booléens qui vaut True pour les nœuds proches de \(x = 1\).

Consigne

Identifiez, avec des tableaux nodes_left et nodes_right les nœuds à droite et à gauche de la barre.

Une fois les nœuds identifiés, les conditions aux limites sont relativement simples :

  1. Un déplacement bloqué dans une direction est équivalent à une vitesse nulle et une force nulle dans cette direction, quel que soit le temps \(t\)
  2. Une force appliquée est simplement ajoutée, à chaque pas de temps, à la force que subit la particule (comme la force de gravité)

10.3.3 Connectivité initiale

La connectivité des nœuds est simple à énoncer : chaque nœud est connecté à ses 4 plus proches voisins et ses 4 voisins en diagonales. En revanche, elle est plus difficile à transcrire en code. C’est pourquoi on utilise de nouveau les listes de voisins, avec la fonction suivante :

from matscipy.neighbours import neighbour_list

def get_neighbour_info(positions, cutoff, data='ij',
                       domain=None, periodicity=None):
    """
    Calcule les listes de voisins (optionellement les distances).

    Exemple d'utilisation:
    >>> i, j = get_neighbour_info(positions, rc)
    """
    if domain is None:
        domain = np.max(positions, axis=1) - np.min(positions, axis=1)
        domain[domain == 0] = np.max(domain)

    domain = np.asanyarray(domain)

    # On complète la matrice du domaine
    if domain.ndim == 1:
        domain = np.diag(domain)

    if domain.shape == (2, 2):
        extended_domain = np.eye(3)

        extended_domain[:2, :2] = domain
        domain = extended_domain

    # Non-périodique par défaut
    if periodicity is None:
        periodicity = np.array([False] * 3)

    periodicity = np.asanyarray(periodicity)

    # On complète la périodicité
    if periodicity.shape[0] == 2:
        periodicity = np.concatenate((periodicity, [False]))

    # On rajoute une coordonnée en 2d
    if positions.shape[0] == 2:
        full_positions = np.vstack((positions, np.zeros(positions.shape[1])))
    else:
        full_positions = positions

    # Calcul des voisins avec matscipy
    neigh_data = list(neighbour_list(data, positions=full_positions.T,
                                     cell=domain, pbc=periodicity,
                                     cutoff=float(cutoff)))

    # Transpose rij
    D_index = data.find('D')
    if D_index != -1:
        neigh_data[D_index] = neigh_data[D_index].T

    return neigh_data
Consigne

Utilisez la fonction ci-dessus pour créer la connectivité initiale et calculer la longueur initiale des barres, \(L_{ij}\). Il faudra choisir cutoff judicieusement.

10.3.4 Extension, onde et équilibre

Dans un premier temps, on choisira la raideur des barres \(k_{ij}\) = 10 MN/m pour toutes les barres, et \(m\) = 0.2 kg.

Consigne

Estimez le pas de temps critique pour ce système (expliquez votre calcul).

Avec une vitesse initiale nulle, faites la simulation d’extension d’une barre avec force appliquée sur les nœuds à droite \(\vec f_A = f_A \vec e_x\) et \(f_A\) = 0.1 MN/m1. Comme pour les TP précédents, calculez l’énergie cinétique, potentielle et mécanique. L’énergie mécanique est-elle constante ? Pourquoi ?

Une observation des trajectoires sur Ovito met en évidence une onde de choc qui se propage de la droite vers la gauche. Afin d’éviter les vibrations “thermiques” ou dues aux ondes de chocs dans le système, on utilise le thermostat de Langevin avec une température nulle. Calculez \(\tau\) pour obtenir un amortissement critique (cf. l’oscillateur harmonique amorti vu au semestre précédent dans le cours de vibrations).

Avec le thermostat de Langevin, faites une simulation jusqu’à ce que l’énergie cinétique moyenne par nœud soit petite comparée à l’énergie potentielle (choisissez la tolérance). Notez que si la vitesse initiale des nœuds est nulle, l’énergie cinétique à la première itération est nulle. Il faut donc écrire une condition d’arrêt du type :

for step in range(nsteps):
    if Ec < ... and step > 100:
        break

De cette façon, on force la simulation d’un certain nombre de pas de temps (ici 100, à modifier) avant de vérifier la condition sur l’énergie cinétique. Une autre possibilité pour vérifier l’équilibre est de calculer la norme des forces : à l’équilibre celle-ci est nulle. On pourra donc comparer la norme des forces à une tolérance.

Calculez l’élongation de la barre à partir du déplacement moyen des nœuds à droite de la barre.

Astuce

Afin de tester l’implémentation, on pourra faire une simulation plus simple avec nx = 20 et ny = 1, ce qui est une simple chaine de ressort, dont la raideur équivalente est connue. Il est donc possible d’en déduire le déplacement et l’élongation analytiquement et de comparer avec les valeurs obtenues par simulation.

10.3.5 Module d’élasticité

On cherche à ajuster \(k_{ij}\) pour obtenir le module de Young effectif de la barre \(E_\mathrm Y\) que l’on désire. Pour cela, il nous faut mesurer le module de Young pour la barre avec tous les \(k_{ij}\) = 10 MN/m.

Consigne

Faites varier la force \(f_A\) de 0 kN/m à 1000 kN/m, et affichez sur un graphe l’élongation \(\varepsilon_l\) en abscisse et la contrainte (force totale divisée par la largeur de la barre) en ordonnée. Estimez le module de Young effectif à partir de ce graphe.

Quelle valeur de \(k_{ij}\) doit-on utiliser pour obtenir \(E_\mathrm{Y} \approx\) 200 GPa ? Montrez sur le même graphe la courbe élongation-contrainte avec les raideurs ajustées.

10.3.6 Coefficient de Poisson

Le coefficient de Poisson \(\nu\) donne la relation entre l’élongation longitudinale \(\varepsilon_l\) de la barre (dans la direction de la force) et l’élongation transversale \(\varepsilon_t\) (dans la direction perpendiculaire à la force) :

\[ \varepsilon_t = -\nu \varepsilon_l\]

Consigne

Mesurez l’élongation transversale de la barre et estimez le coefficient de Poisson. De quelles barres faut-il ajuster la raideur pour changer le coefficient de Poisson ? Justifiez.

10.4 Évaluation

L’évaluation intermédiaire (comptant pour 30% de la note finale) se fait sur la base de deux productions scientifiques à partir du contenu de ce TP :

  • Un rapport de TP (15%)
  • Une présentation orale (15%)

Le travail est à faire en binôme.

10.4.1 Rapport

Le rapport (de 3 à 5 pages environ) devra :

  • Énoncer les objectifs scientifiques du TP
  • Expliquer la méthodologie employée pour atteindre ces objectifs (modèle et méthodes de simulation, valeurs des paramètres)
  • Présenter et discuter les résultats (figures)
  • Conclure

Servez-vous des questions et consignes dans ce sujet de TP pour écrire le rapport (avec vos propres mots). Les figures doivent être lisibles2, avec nom des axes, symboles des quantités, unités et légende si plusieurs courbes sont montrées.

Consigne

Le rapport est à uploader sur Moodle au format PDF obligatoirement, avant la présentation orale.

10.4.2 Présentation

La présentation (5 minutes + 2 minutes pour les questions) devra suivre la même structure que le rapport, cependant, la transmission d’une information à l’oral ne se fait pas de la même manière qu’à l’écrit. Quelques conseils :

  • Évitez les phrases complètes sur vos slides,
  • Évitez, si possible, d’utiliser des symboles dans les figures,
  • Toutes les figures sur un slide doivent être commentées à l’oral,
  • Les figures doivent être lisibles depuis le fond (attention à la taille de la police d’écriture), il faut souvent modifier les figures du rapport pour une présentation,
  • Entrainez-vous à respecter le temps de présentation.

Les présentations se font devant la classe, afin de :

  • Recevoir des commentaires sur votre propre présentation
  • Apprendre des autres présentations

Pour celles et ceux qui le souhaiteraient, il est possible de faire la présentation en anglais.

Consigne

Le fichier de la présentation est à uploader sur Moodle au format PDF obligatoirement, avant la présentation orale.

  1. Comme on fait une simulation 2D d’un solide, on divise la force par l’épaisseur du solide dans la dimension \(z\), supposée petite par rapport aux autres dimensions. La “force” appliquée à l’extrémité est donc une force par unité de longueur.↩︎

  2. Il existe deux familles de fichiers pour les images : les images bitmap (ou raster), qui sont un tableau de pixels, et les images vectorielles, qui contiennent la définition de formes géométriques. Ces dernières sont à préférer dans le rapport, puisqu’elles ne sont pas pixelisées. Si le rapport est rédigé sous Word ou LibreOffice, le format SVG doit être utilisé (plt.savefig('figure.svg')). Si le rapport est rédigé sous LaTeX (avec Overleaf par exemple), le format à utiliser est le PDF (plt.savefig('figure.pdf')).↩︎