9  Application numérique

Résolution par éléments finis (FEniCSx) du problème de concentration de contrainte autour d’un trou circulaire dans une plaque soumise à une traction uniforme.

9.1 Importation des librairies

Organisation des imports par catégorie : maillage (gmsh), éléments finis (dolfinx), calcul scientifique (numpy), visualisation (pyvista, matplotlib).

9.2 Géométrie et maillage

9.2.1 Création de la géométrie avec Gmsh

On définit une plaque carrée de côté \(L = 0,4\) m percée d’un trou central de rayon \(a = 0,02\) m. Un champ de taille de maille permet un raffinement local autour du trou (\(h_\text{trou} = 0,004\) m) pour capturer le gradient de contrainte, tandis que le reste de la plaque utilise une taille globale \(h_\text{global} = 0,04\) m.

Les bords sont identifiés géométriquement pour l’application des conditions aux limites : - left (tag 4) : bord \(x = 0\) - bottom (tag 1) : bord \(y = 0\) - right (tag 2) : bord \(x = L\) - top (tag 3) : bord \(y = L\)

Code 
L, R = 0.4, 0.02
h_global = 0.04
h_hole = 0.004
gdim = 2
fdim = 1

gmsh.initialize()
occ = gmsh.model.occ
mesh_comm = MPI.COMM_WORLD
model_rank = 0

if mesh_comm.rank == model_rank:
    plate = occ.addRectangle(0.0, 0.0, 0.0, L, L)
    notch = occ.addDisk(0.0, 0.0, 0.0, R, R)
    outDimTags, _ = occ.cut([(gdim, plate)], [(gdim, notch)])
    perf_plate = outDimTags[0][1]
    occ.synchronize()

    eps = 1e-3
    left = gmsh.model.getEntitiesInBoundingBox(
        -eps, -eps, -eps, eps, L + eps, eps, dim=fdim)[0][1]
    bottom = gmsh.model.getEntitiesInBoundingBox(
        -eps, -eps, -eps, L + eps, eps, eps, dim=fdim)[0][1]
    right = gmsh.model.getEntitiesInBoundingBox(
        L - eps, -eps, -eps, L + eps, L + eps, eps, dim=fdim)[0][1]
    top = gmsh.model.getEntitiesInBoundingBox(
        -eps, L - eps, -eps, L + eps, L + eps, eps, dim=fdim)[0][1]

    gmsh.model.addPhysicalGroup(gdim, [perf_plate], 1)
    gmsh.model.addPhysicalGroup(fdim, [bottom], 1)
    gmsh.model.addPhysicalGroup(fdim, [right], 2)
    gmsh.model.addPhysicalGroup(fdim, [top], 3)
    gmsh.model.addPhysicalGroup(fdim, [left], 4)

    boundaries = gmsh.model.getBoundary([(gdim, perf_plate)], oriented=False)
    hole_curves = []
    for dim, tag in boundaries:
        com = gmsh.model.occ.getCenterOfMass(dim, tag)
        is_outer = (abs(com[0]) < eps or abs(com[0] - L) < eps
                    or abs(com[1]) < eps or abs(com[1] - L) < eps)
        if not is_outer:
            hole_curves.append(tag)

    field = gmsh.model.mesh.field
    field.add("Distance", 1)
    field.setNumbers(1, "CurvesList", hole_curves)
    field.add("Threshold", 2)
    field.setNumber(2, "InField", 1)
    field.setNumber(2, "SizeMin", h_hole)
    field.setNumber(2, "SizeMax", h_global)
    field.setNumber(2, "DistMin", R)
    field.setNumber(2, "DistMax", 4 * R)
    field.setAsBackgroundMesh(2)

    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h_hole)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h_global)
    gmsh.model.mesh.generate(gdim)

    domain, cells, facets = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()
Info    : Meshing 1D...                                                                                      
Info    : [  0%] Meshing curve 1 (Ellipse)
Info    : [ 30%] Meshing curve 2 (Line)
Info    : [ 50%] Meshing curve 3 (Line)
Info    : [ 70%] Meshing curve 4 (Line)
Info    : [ 90%] Meshing curve 5 (Line)
Info    : Done meshing 1D (Wall 0.00159512s, CPU 0.00145s)
Info    : Meshing 2D...
Info    : Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.00348843s, CPU 0.0035s)
Info    : 302 nodes 607 elements

9.2.2 Visualisation du maillage

Le maillage est affiché avec PyVista. Le raffinement concentrique autour du trou est visible.

9.3 Espace fonctionnel

On utilise un espace d’éléments finis vectoriel de type Lagrange d’ordre 2 (P2) pour approcher le champ de déplacement \(\mathbf{u} = (u_x, u_y)\).

V = fem.functionspace(domain, ("P", 2, (gdim,)))

9.4 Propriétés matériau

Le matériau est un acier élastique linéaire isotrope. En contraintes planes, la matrice de compliance \(\mathbf{S}\) relie le tenseur des déformations vectorialisées \((\varepsilon_{11}, \varepsilon_{22}, 2\varepsilon_{12})\) au tenseur des contraintes vectorialisées \((\sigma_{11}, \sigma_{22}, \sigma_{12})\).

E = 210.0
nu = 0.3
G = E / (2.0 * (1.0 + nu))

S = ufl.as_matrix([
    [1.0/E, -nu/E, 0.0],
    [-nu/E, 1.0/E, 0.0],
    [0.0,   0.0,   1.0/G]
])
C = ufl.inv(S)

9.4.1 Tenseurs de déformation et de contrainte

Le tenseur de déformation linéarisé \(\boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T)\) est défini via UFL. La contrainte s’en déduit par la loi de Hooke : \(\boldsymbol{\sigma} = \mathbf{C} : \boldsymbol{\varepsilon}\).

9.5 Conditions aux limites

On applique des conditions de Dirichlet pour bloquer les déplacements rigides : - Blocage horizontal (\(u_x = 0\)) sur le bord gauche (tag 4) - Blocage vertical (\(u_y = 0\)) sur le bord inférieur (tag 1)

9.6 Chargement

Une traction uniforme \(T = 10\) MPa est appliquée sur le bord supérieur (tag 3), sous forme d’une pression normale \(\mathbf{T} = T \,\mathbf{n}\).

ds = ufl.Measure("ds", domain=domain, subdomain_data=facets)
T = fem.Constant(domain, 10.0)
n = ufl.FacetNormal(domain)

9.7 Formulation variationnelle

La forme bilinéaire \(a(\mathbf{u}, \mathbf{v})\) représente l’énergie de déformation virtuelle, et la forme linéaire \(L(\mathbf{v})\) le travail des forces extérieures :

\[a(\mathbf{u}, \mathbf{v}) = \int_\Omega \boldsymbol{\sigma}(\mathbf{u}) : \boldsymbol{\varepsilon}(\mathbf{v}) \, d\Omega\] \[L(\mathbf{v}) = \int_{\Gamma_\text{top}} T \,\mathbf{n} \cdot \mathbf{v} \, d\Gamma\]

du = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
u = fem.Function(V, name="Displacement")

a_form = ufl.inner(stress(du), strain(u_)) * ufl.dx
L_form = ufl.dot(T * n, u_) * ds(3)

9.8 Résolution

Le problème linéaire est résolu avec le solveur PETSc par défaut (LU).

problem = fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
problem.solve()
Coefficient(FunctionSpace(Mesh(blocked element (Basix element (P, triangle, 1, equispaced, unset, False, float64, []), (2,)), 0), blocked element (Basix element (P, triangle, 2, gll_warped, unset, False, float64, []), (2,))), 2)

9.9 Post-traitement : projection des contraintes

Le tenseur des contraintes est projeté sur un espace DG d’ordre 0 (constant par cellule), puis sur un espace CG d’ordre 1 (continu) pour accéder aux valeurs nodales. On reconstruit la contrainte orthoradiale \(\sigma_{\theta\theta}\) par changement de base polaire :

\[\sigma_{\theta\theta} = \sigma_{xx}\sin^2\theta + \sigma_{yy}\cos^2\theta - 2\sigma_{xy}\sin\theta\cos\theta\]

9.10 Résultats au bord du trou

On filtre les nœuds situés sur le bord du trou (\(r \approx a\)) et on calcule \(\sigma_{\theta\theta}\) en ces points. La comparaison avec la solution de Kirsch permet de valider le modèle.

RÉSULTATS AU BORD DU TROU (r = 0.02) :
MAX : 30.5652 MPa
MIN : -10.4058 MPa

9.11 Visualisation du champ de contrainte \(\sigma_{rr}\) et \(\sigma_{r\theta}\)

sigma_rr_all = (
    sig_vals[:, 0] * cos_t**2
    + sig_vals[:, 1] * sin_t**2
    + 2 * sig_vals[:, 2] * sin_t * cos_t
)
grid.point_data["sigma_rr"] = sigma_rr_all

sigma_rt_all = (
    (sig_vals[:, 1] - sig_vals[:, 0]) * sin_t * cos_t
    + sig_vals[:, 2] * (cos_t**2 - sin_t**2)
)
grid.point_data["sigma_rt"] = sigma_rt_all

p = pv.Plotter(shape=(1, 2))

grid.point_data["sigma_rr"] = sigma_rr_all
grid.point_data["sigma_rt"] = sigma_rt_all

p.subplot(0, 0)

p.add_mesh(
    grid,
    scalars="sigma_rr",
    show_edges=False,
    cmap="inferno",
    scalar_bar_args={"title": "sigma_rr (MPa)","title_font_size": 24, "color": "grey", "label_font_size": 18, "font_family": "arial" }
)

p.add_text(
    "Contrainte de radiale",
    font_size=12,
    color="grey",
    position="upper_edge",
)

p.add_axes(color="grey")
p.show_bounds(color="grey")
p.set_background(bg)

p.subplot(0, 1)

p.add_mesh(
    grid,
    scalars="sigma_rr",
    show_edges=False,
    cmap="inferno",
    scalar_bar_args={"title": "sigma_rr (MPa)","title_font_size": 24, "color": "grey", "label_font_size": 18, "font_family": "arial" }
)

p.add_text(
    "Contrainte de cisaillement",
    font_size=12,
    color="grey",
    position="upper_edge",
)

p.add_axes(color="grey")
p.show_bounds(color="grey")
p.set_background(bg)

p.link_views()
p.show()

9.12 Visualisation du champ de contrainte \(\sigma_{\theta\theta}\)

Le champ de contrainte orthoradiale est affiché sur l’ensemble du domaine. La concentration au sommet du trou (\(\sigma_{\theta\theta} \approx 30\) MPa) est clairement visible.

9.13 Comparaison avec la solution théorique

La solution de Kirsch donne les valeurs exactes au bord du trou :

  • Point A (\(\theta = \pi/2\)) : \(\sigma_{\theta\theta} = 3\sigma = 30\) MPa
  • Point B (\(\theta = 0\)) : \(\sigma_{\theta\theta} = -\sigma = -10\) MPa

Le facteur de concentration de contrainte théorique est \(K_t = 3\).

Code 
sigma_th = 10.0
max_theo = 3 * sigma_th
min_theo = -sigma_th

if np.any(mask_border):
    max_num = np.max(sigma_tt_border)
    min_num = np.min(sigma_tt_border)
    err_max = abs(max_num - max_theo) / abs(max_theo) * 100
    err_min = abs(min_num - min_theo) / abs(min_theo) * 100
    print(f"Point A (θ=π/2) : Kirsch = {max_theo} MPa, FEM = {max_num:.4f} MPa, erreur = {err_max:.2f}%")
    print(f"Point B (θ=0)   : Kirsch = {min_theo} MPa, FEM = {min_num:.4f} MPa, erreur = {err_min:.2f}%")
    print(f"Facteur K_t     : Kirsch = 3.00, FEM = {max_num/sigma_th:.4f}")

9.14 Conclusion

La simulation par éléments finis reproduit avec précision la solution analytique de Kirsch :

  • L’erreur relative sur la contrainte maximale est inférieure à \(2\%\).
  • Le facteur de concentration \(K_t = 3\) est correctement capturé.
  • La méthode FEM, validée sur ce cas de référence, peut être étendue à des géométries plus complexes (entailles, fissures) sans solution analytique.