ECOLE POLYTECHNIQUE DE OUAGADOUGOU

UFR / FILIÈRE : IGSIT/ GM1

Année universitaire : 2025 – 2026

PROJET DE
FIABILITÉ D’UNE STRUCTURE 2D PAR LA MÉTHODE DES
ÉLÉMENTS FINIS,
SIMULATION DE MONTE CARLO ET MÉTHODE FORM

Présenté par : N’DO Cheick Oumar

Encadré par : Docteur ZIO

1
Résumé
Ce mémoire présente une étude de fiabilité appliquée à une structure plane soumise à une charge de
compression. La structure est modélisée par la méthode des éléments finis 2D. La défaillance est définie à
partir de la contrainte de Von Mises maximale comparée à une contrainte limite du matériau. Deux
approches probabilistes sont développées et comparées : la simulation de Monte Carlo, qui fournit une
estimation directe de la probabilité de défaillance, et la méthode FORM (First Order Reliability Method), qui
donne une approximation rapide via l’indice de fiabilité β. Le document détaille la formulation mécanique, la
construction de la fonction d’état limite, les algorithmes numériques, une implémentation Python
commentée, ainsi qu’une discussion sur les avantages, limites et conditions d’emploi des deux méthodes.

2
Introduction générale
En ingénierie, le calcul déterministe classique consiste à fixer une valeur unique pour chaque paramètre
d’entrée du modèle (propriétés du matériau, dimensions, chargement, contraintes admissibles, etc.) puis à
résoudre le problème mécanique. Cette démarche reste utile, mais elle ne tient pas compte de l’incertitude
inhérente aux systèmes réels. Or, dans la pratique, les matériaux présentent une dispersion de leurs
propriétés, les charges ne sont pas parfaitement connues, les conditions d’exploitation varient, et les mesures
expérimentales sont elles-mêmes entachées d’erreurs.

La fiabilité structurale répond à cette difficulté en représentant certains paramètres comme des variables
aléatoires et en quantifiant la probabilité qu’un critère de sécurité soit violé. Dans ce cadre, la notion centrale
est la fonction d’état limite g(X), qui sépare l’ensemble des états sûrs de l’ensemble des états défaillants. Le
calcul de la probabilité de défaillance P_f = P[g(X) ≤ 0] devient alors un objectif majeur.

Le présent projet traite ce problème sur un cas simple mais pédagogique : une barre verticale en acier,
encastrée à sa base et soumise à une charge axiale. Le modèle mécanique est résolu par éléments finis 2D. La
défaillance est évaluée à partir de la contrainte de Von Mises maximale. Deux méthodes probabilistes sont
ensuite mises en œuvre : Monte Carlo et FORM.

Chapitre 1 – Problématique, objectifs et énoncé reformulé

1.1 Contexte du projet

Le sujet initial demande : (i) un résumé de la méthode FORM et de son application, (ii) un résumé de la
méthode de Monte Carlo et une implémentation appliquée à un solide, puis (iii) la résolution du problème
mécanique par éléments finis 2D en y appliquant la fiabilité.

1.2 Énoncé reformulé et exploitable

On considère une barre rectangulaire verticale de hauteur L = 1,0 m et de largeur b = 0,05 m, en acier,
modélisée en deux dimensions dans l’hypothèse des contraintes planes. La base de la barre est encastrée, ce
qui impose les conditions aux limites u_x = u_y = 0 sur y = 0. Une force de compression F est appliquée sur
le bord supérieur. Le coefficient de Poisson est pris égal à ν = 0,30.

Le module d’Young E du matériau est modélisé comme une variable aléatoire uniforme E ~ U(E_min,
E_max), avec E_min = 190 GPa et E_max = 210 GPa. On choisit une contrainte limite du matériau σ_lim =
250 MPa. Le critère de défaillance est défini par :

g(X) = σ_lim − σ_VM,max(X)

où σ_VM,max(X) est la contrainte de Von Mises maximale dans la structure. La structure est dite défaillante
lorsque g(X) ≤ 0.

1.3 Objectifs
 formuler le problème mécanique en éléments finis 2D ;
 calculer les déplacements et la contrainte de Von Mises maximale ;
 construire la fonction d’état limite ;
 estimer la probabilité de défaillance par Monte Carlo ;

3
 estimer la même probabilité par la méthode FORM ;
 comparer les deux approches du point de vue coût / précision / interprétation.

1.4 Hypothèse importante sur le caractère aléatoire

Il faut signaler un point méthodologique important : si l’on ne rend aléatoire que le module d’Young E tout
en imposant une force F déterministe, la contrainte peut devenir peu sensible à E en élasticité linéaire. Le
déplacement dépend alors de E, mais pas nécessairement la contrainte finale. Pour cette raison, le rapport
distingue deux cas :

 Cas A : version conforme à la consigne — E aléatoire et F déterministe ;

 Cas B : version enrichie et plus réaliste — E aléatoire et F aléatoire, par exemple F ~ N(600 kN, (40
kN)^2).

Chapitre 2 – Modélisation mécanique par éléments finis 2D

2.1 Géométrie, maillage et degrés de liberté

La géométrie retenue est un rectangle [0, b] × [0, L]. Le domaine est discrétisé en triangles linéaires à trois
nœuds (éléments CST, Constant Strain Triangle). Chaque nœud porte deux degrés de liberté : le déplacement
horizontal u_x et le déplacement vertical u_y.

Figure 1 – Schéma du solide étudié : barre verticale encastrée et soumise à une force de compression F.

2.2 Hypothèse de contraintes planes

Le problème est traité en contraintes planes. La matrice constitutive du matériau isotrope est :

D = E/(1−ν²) · [[1, ν, 0], [ν, 1, 0], [0, 0, (1−ν)/2]]

4
Cette matrice relie le vecteur des déformations ε = [ε_x, ε_y, γ_xy]^T au vecteur des contraintes σ = [σ_x,
σ_y, τ_xy]^T.

2.3 Élément triangulaire linéaire

Pour un triangle de sommets (x1,y1), (x2,y2), (x3,y3), les déplacements sont interpolés linéairement à partir
des valeurs nodales. La matrice de déformation B est constante sur l’élément, ce qui simplifie l’assemblage.
La matrice de rigidité élémentaire est donnée par :

K_e = t · A · B^T D B

où A est l’aire du triangle et t l’épaisseur du modèle 2D.

2.4 Assemblage global et conditions aux limites

Les matrices élémentaires sont assemblées dans une matrice globale K. Le chargement est représenté par un
vecteur global f, construit en répartissant la force de compression F sur les nœuds du bord supérieur.
L’encastrement de la base est imposé en bloquant les deux composantes de déplacement de tous les nœuds
situés sur y = 0.

Le système à résoudre est alors :

Ku=f

2.5 Contrainte de Von Mises et critère de défaillance

Après calcul du vecteur de déplacement global u, on détermine les déformations puis les contraintes dans
chaque élément. En contraintes planes, la contrainte équivalente de Von Mises s’écrit :

σ_VM = √(σ_x² − σ_xσ_y + σ_y² + 3τ_xy²)

On extrait ensuite la valeur maximale σ_VM,max sur l’ensemble des éléments. Le critère de sécurité est
formulé par la fonction d’état limite :

g(X) = σ_lim − σ_VM,max(X)

Chapitre 3 – Fiabilité structurale et fonction d’état limite

3.1 Définition
La fiabilité d’une structure est la probabilité que celle-ci reste dans le domaine sûr. Si l’on note D = {X :
g(X) ≤ 0} le domaine de défaillance, la probabilité de défaillance est :

P_f = P[g(X) ≤ 0]

La difficulté principale vient du fait que g(X) dépend ici d’un calcul éléments finis. Il n’existe donc
généralement pas d’expression analytique simple de P_f, d’où l’intérêt des méthodes numériques.

3.2 Variables aléatoires du projet

Variable Loi Paramètres Rôle physique

E Uniforme [190 ; 210] GPa Rigidité du matériau

5
 F (option enrichie) Normale μ = 600 kN, σ = 40 kN Charge appliquée

σ_lim Déterministe 250 MPa Seuil de rupture /

sécurité

3.3 Interprétation physique de g(X)

Lorsque g(X) > 0, la contrainte équivalente reste inférieure à la limite admissible : la structure est considérée
sûre. Lorsque g(X) = 0, la structure se trouve sur la frontière de défaillance. Enfin, lorsque g(X) < 0, le
niveau de contrainte dépasse le seuil admissible : on considère que la structure est défaillante.

Chapitre 4 – Méthode de Monte Carlo

4.1 Principe
La méthode de Monte Carlo est conceptuellement la plus simple. On génère N réalisations aléatoires des
variables d’entrée, on résout le problème mécanique pour chaque réalisation, puis on compte la fréquence
des cas défaillants.

Si I_i = 1 lorsque g(X_i) ≤ 0 et 0 sinon, alors :

P_f ≈ (1/N) Σ I_i

Cette méthode est robuste et peu sensible à la forme de g(X), mais elle peut devenir coûteuse si chaque
évaluation de g(X) nécessite un calcul EF et si la probabilité de défaillance est faible.

4.2 Algorithme détaillé

1. Générer un échantillon aléatoire de E (et éventuellement de F).
2. Construire la matrice de rigidité K(E).
3. Appliquer les conditions aux limites et résoudre le système K u = f.
4. Calculer la contrainte de Von Mises dans chaque élément.
5. Extraire σ_VM,max puis évaluer g(X).
6. Compter le nombre de cas tels que g(X) ≤ 0.
7. Estimer P_f et, si nécessaire, l’intervalle de confiance.

6
 Figure 2 – Illustration typique de la convergence de l’estimateur Monte Carlo de la probabilité de défaillance.

4.3 Avantages et limites

 Avantage : très générale, facile à comprendre et à programmer.
 Avantage : donne une référence numérique fiable si N est suffisamment grand.
 Limite : temps de calcul élevé lorsque le modèle EF est coûteux.
 Limite : inefficace pour les probabilités très faibles sans techniques d’accélération (importance sampling,
subset simulation, etc.).

Chapitre 5 – Méthode FORM (First Order Reliability Method)

5.1 Idée générale

FORM cherche le point le plus probable de défaillance, appelé point de conception, sur la surface g(X)=0.
Pour cela, les variables physiques X sont transformées dans un espace de variables normales réduites U. La
surface de défaillance est ensuite approchée localement par son plan tangent au voisinage du point de
conception.

L’indice de fiabilité β est la distance entre l’origine de l’espace réduit et le point de conception u*.
L’approximation FORM de la probabilité de défaillance est alors :

P_f ≈ Φ(−β)

7
 Figure 3 – Interprétation géométrique de FORM dans l’espace normal réduit.

5.2 Transformation des variables

Dans une implémentation simple, une variable uniforme E ~ U(E_min, E_max) peut être remplacée par une
variable normale équivalente de même moyenne et même écart-type :

μ_E = (E_min + E_max)/2

σ_E = (E_max − E_min)/√12

u = (E − μ_E) / σ_E

Cette approximation permet de mettre en place une FORM de premier niveau sans alourdir excessivement le
projet.

5.3 Procédure FORM simplifiée pour ce projet

8. Choisir un point initial X^(0), généralement la moyenne des variables.
9. Évaluer g(X^(k)).
10. Calculer numériquement le gradient de g au point courant.
11. Mettre à jour le point de conception dans l’espace réduit.
12. Répéter jusqu’à convergence de β ou du point de conception.
13. En déduire P_f ≈ Φ(−β).

8
5.4 Interprétation de β
Plus β est grand, plus la structure est éloignée du domaine de défaillance. À titre indicatif : β ≈ 1 correspond
à une structure faiblement fiable, β ≈ 3 correspond à un niveau de fiabilité élevé, et β > 4 indique
généralement une structure très sûre vis-à-vis du critère étudié.

Chapitre 6 – Implémentation Python détaillée

Le script Python présenté ci-dessous est structuré en blocs : génération du maillage, assemblage éléments
finis, résolution mécanique, calcul de la contrainte de Von Mises, simulation de Monte Carlo et
approximation FORM. Chaque fonction est commentée afin de pouvoir être expliquée pendant la
soutenance.

# ============================================================
# PROJET : Fiabilité d'une structure 2D par EF + Monte Carlo + FORM
# ============================================================

import numpy as np
from math import erf

# ------------------------------------------------------------
# 1) Fonction utilitaire : CDF de la loi normale centrée réduite
# ------------------------------------------------------------
def normal_cdf(x):
return 0.5 * (1.0 + erf(x / [Link](2.0)))

# ------------------------------------------------------------
# 2) Génération d'un maillage triangulaire d'un rectangle
# ------------------------------------------------------------
def generate_mesh(L=1.0, b=0.05, nx=6, ny=20):
xs = [Link](0.0, b, nx + 1)
ys = [Link](0.0, L, ny + 1)
nodes = [Link]([[x, y] for y in ys for x in xs], dtype=float)

def nid(i, j):

return j * (nx + 1) + i

elems = []
for j in range(ny):
for i in range(nx):
n1, n2 = nid(i, j), nid(i + 1, j)
n3, n4 = nid(i, j + 1), nid(i + 1, j + 1)
[Link]([n1, n2, n4])
[Link]([n1, n4, n3])
return nodes, elems

# ------------------------------------------------------------
# 3) Matrice B et aire d'un triangle CST
# ------------------------------------------------------------
def triangle_B_matrix(x1, y1, x2, y2, x3, y3):
A = 0.5 * [Link]([Link]([[1, x1, y1], [1, x2, y2], [1, x3, y3]]))
if A <= 0:
raise ValueError("Triangle dégénéré ou orientation incorrecte")

b1, b2, b3 = y2-y3, y3-y1, y1-y2

c1, c2, c3 = x3-x2, x1-x3, x2-x1

9
 B = (1/(2*A))*[Link]([
[b1, 0, b2, 0, b3, 0],
[0, c1, 0, c2, 0, c3],
[c1, b1, c2, b2, c3, b3]
], dtype=float)
return B, A

# ------------------------------------------------------------
# 4) Loi de comportement : contraintes planes
# ------------------------------------------------------------
def constitutive_matrix(E, nu):
coef = E / (1.0 - nu**2)
return coef * [Link]([
[1.0, nu, 0.0],
[nu, 1.0, 0.0],
[0.0, 0.0, (1.0 - nu) / 2.0]
], dtype=float)

# ------------------------------------------------------------
# 5) Assemblage de la matrice globale K
# ------------------------------------------------------------
def assemble_system(nodes, elems, E, nu, thickness=1.0):
ndof = 2 * len(nodes)
K = [Link]((ndof, ndof))
D = constitutive_matrix(E, nu)

for e in elems:
n1, n2, n3 = e
x1, y1 = nodes[n1]
x2, y2 = nodes[n2]
x3, y3 = nodes[n3]

B, A = triangle_B_matrix(x1, y1, x2, y2, x3, y3)

ke = thickness * A * (B.T @ D @ B)

dofs = [2*n1,2*n1+1, 2*n2,2*n2+1, 2*n3,2*n3+1]

for i in range(6):
for j in range(6):
K[dofs[i], dofs[j]] += ke[i,j]
return K

# ------------------------------------------------------------
# 6) Vecteur de charge et conditions aux limites
# ------------------------------------------------------------
def build_load_vector(nodes, F, L):
f = [Link](2*len(nodes))
top_nodes = [i for i,(x,y) in enumerate(nodes) if [Link](y, L)]
for i in top_nodes:
f[2*i+1] += -F / len(top_nodes)
return f

def apply_boundary_conditions(K, f, nodes):

fixed = []
for i,(x,y) in enumerate(nodes):
if [Link](y, 0.0):
[Link]([2*i, 2*i+1])

all_dofs = [Link]([Link][0])
free = np.setdiff1d(all_dofs, fixed)
Kff = K[np.ix_(free, free)]
ff = f[free]

10
 return Kff, ff, free

# ------------------------------------------------------------
# 7) Résolution EF et calcul de Von Mises
# ------------------------------------------------------------
def solve_fem(E, F, L=1.0, b=0.05, nu=0.30, thickness=1.0, nx=6, ny=20):
nodes, elems = generate_mesh(L=L, b=b, nx=nx, ny=ny)
K = assemble_system(nodes, elems, E, nu, thickness)
f = build_load_vector(nodes, F, L)
Kff, ff, free = apply_boundary_conditions(K, f, nodes)

u = [Link]([Link][0])
u[free] = [Link](Kff, ff)

D = constitutive_matrix(E, nu)
vm_max = 0.0
for e in elems:
n1, n2, n3 = e
x1,y1 = nodes[n1]
x2,y2 = nodes[n2]
x3,y3 = nodes[n3]
B, A = triangle_B_matrix(x1,y1,x2,y2,x3,y3)
dofs = [2*n1,2*n1+1, 2*n2,2*n2+1, 2*n3,2*n3+1]
ue = u[dofs]
strain = B @ ue
sx, sy, txy = D @ strain
vm = [Link](sx**2 - sx*sy + sy**2 + 3*txy**2)
vm_max = max(vm_max, vm)

return u, vm_max

# ------------------------------------------------------------
# 8) Fonction d'état limite
# ------------------------------------------------------------
def limit_state(E, F=6e5, sigma_lim=250e6):
_, sigma_vm_max = solve_fem(E=E, F=F)
return sigma_lim - sigma_vm_max

# ------------------------------------------------------------
# 9) Monte Carlo
# ------------------------------------------------------------
def monte_carlo_pf(N=2000, Emin=190e9, Emax=210e9,
F_mean=6e5, F_std=None, sigma_lim=250e6, seed=123):
rng = [Link].default_rng(seed)
n_fail = 0

for _ in range(N):
E = [Link](Emin, Emax)
F = F_mean if F_std is None else [Link](F_mean, F_std)
g = limit_state(E, F=F, sigma_lim=sigma_lim)
if g <= 0:
n_fail += 1

return n_fail / N

# ------------------------------------------------------------
# 10) FORM simplifiée 1D
# ------------------------------------------------------------
def form_1d(Emin=190e9, Emax=210e9, F=6e5, sigma_lim=250e6):
mu_E = 0.5*(Emin + Emax)
sig_E = (Emax - Emin)/[Link](12)

11
 def g(E):
return limit_state(E, F=F, sigma_lim=sigma_lim)

grid = [Link](Emin, Emax, 25)

vals = [Link]([g(E) for E in grid])

if [Link](vals > 0):

return {"beta": [Link], "pf": 0.0}
if [Link](vals <= 0):
return {"beta": -[Link], "pf": 1.0}

# recherche d'un intervalle de changement de signe

a, b = None, None
for i in range(len(grid)-1):
if vals[i] * vals[i+1] <= 0:
a, b = grid[i], grid[i+1]
break

# dichotomie
for _ in range(60):
m = 0.5*(a+b)
if g(a)*g(m) <= 0:
b = m
else:
a = m

E_star = 0.5*(a+b)
beta = abs((E_star - mu_E)/sig_E)
pf = normal_cdf(-beta)
return {"E*": E_star, "beta": beta, "pf": pf}

Chapitre 7 – Exploitation, résultats et interprétation

7.1 Plan d’exploitation conseillé

14. Exécuter le code avec F déterministe et E aléatoire pour vérifier le fonctionnement global du solveur EF.
15. Tracer la variation de la contrainte de Von Mises maximale en fonction de E.
16. Lancer Monte Carlo avec N = 1000, puis 5000, puis 10000 pour observer la convergence de P_f.
17. Introduire ensuite l’incertitude sur F afin d’obtenir un problème de fiabilité plus riche.
18. Comparer P_f obtenu par Monte Carlo à celui estimé par FORM.

7.3 Résultats
Dans le cas où seule la rigidité E est aléatoire, il est possible que la variabilité de la contrainte soit faible,
surtout si le chargement est imposé en force. Dans ce cas, la probabilité de défaillance peut être quasi nulle
ou peu sensible à la dispersion de E. En revanche, si la charge F est également aléatoire, l’effet sur
σ_VM,max devient plus marqué et l’étude de fiabilité prend tout son sens.

Monte Carlo doit être présenté comme la méthode de référence, car elle ne suppose pas de linéarisation de la
surface de défaillance. FORM, elle, fournit une approximation plus rapide et permet d’interpréter la fiabilité
via l’indice β.

12
Conclusion générale
Ce projet met en évidence le lien direct entre la mécanique des structures, la modélisation numérique et
l’analyse probabiliste. La méthode des éléments finis fournit la réponse mécanique de la structure, la
simulation de Monte Carlo permet d’estimer directement la probabilité de défaillance, et la méthode FORM
propose une approximation plus économique en temps de calcul. Le travail montre également qu’un projet
de fiabilité doit être soigneusement formulé : le choix des variables aléatoires et de la fonction d’état limite
conditionne fortement la pertinence physique des résultats.

Comme prolongements possibles, on peut envisager : l’introduction de plusieurs variables aléatoires, un

maillage plus fin, une formulation 3D, une méthode FORM multidimensionnelle complète de type Hasofer–
Lind–Rackwitz–Fiessler, ou encore des méthodes de simulation avancées pour les faibles probabilités de
défaillance.

Bibliographie
 Amiri, Mehdi_ Jackson, Christopher_ Modarres, Mohammad - Probabilistic

13