Problème
Étude 2D déformations planes d’un parallélogramme incliné. - Géométrie: L=20m, H=55m, angle α avec la verticale - Matériau: béton E=30GPa, ν=0.25, ρ=2000kg/m³ - Contact: base sur sol rigide (unilatéral), poids propre ## Objectifs 1. Champ de contraintes 2. α décollement base 3. α critique perte équilibre
import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from dolfinx import mesh, fem, io, plot
from dolfinx.fem.petsc import LinearProblem
import ufl
from ufl import dx, inner, grad, sym, dot, tr, Identity, as_vector, sqrt
import gmsh
from dolfinx.io import gmshio
import pyvista as pv
from myst_nb import glue
from IPython.display import clear_output
import panel as pn
pn.extension('vtk')
pv.OFF_SCREEN = True
import warnings
warnings.filterwarnings('ignore')
plt.style.use('dark_background')
bg=(40/255,43/255,48/255,1.0)
fg='grey'
print('OK')
L=20.; H=55.; E=30e9; nu=0.25; rho=2000.; g=9.81
mu=E/(2*(1+nu)); lmbda=E*nu/((1+nu)*(1-2*nu))
h=0.8
print('Parametres OK')
def create_mesh_gmsh(a):
alpha=np.radians(a); d=H*np.tan(alpha)
gmsh.initialize(); gmsh.model.add('tour')
gmsh.option.setNumber('General.Verbosity', 0)
gmsh.option.setNumber('Mesh.Algorithm', 8)
lc=h
p1=gmsh.model.occ.addPoint(0,0,0,lc)
p2=gmsh.model.occ.addPoint(L,0,0,lc)
p3=gmsh.model.occ.addPoint(L+d,H,0,lc)
p4=gmsh.model.occ.addPoint(d,H,0,lc)
l1=gmsh.model.occ.addLine(p1,p2)
l2=gmsh.model.occ.addLine(p2,p3)
l3=gmsh.model.occ.addLine(p3,p4)
l4=gmsh.model.occ.addLine(p4,p1)
cl=gmsh.model.occ.addCurveLoop([l1,l2,l3,l4])
s=gmsh.model.occ.addPlaneSurface([cl])
gmsh.model.occ.synchronize()
gmsh.model.mesh.setRecombine(2,s)
gmsh.option.setNumber('Mesh.RecombinationAlgorithm', 2)
gmsh.model.addPhysicalGroup(1,[l1],1)
gmsh.model.addPhysicalGroup(1,[l2,l4],2)
gmsh.model.addPhysicalGroup(1,[l3],3)
gmsh.model.addPhysicalGroup(2,[s],1)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.optimize('HighOrder')
msh,ct,ft=gmshio.model_to_mesh(gmsh.model,MPI.COMM_WORLD,0,gdim=2)
gmsh.finalize()
return msh,ft
msh_t,ft_t=create_mesh_gmsh(10.)
nc=msh_t.topology.index_map(msh_t.topology.dim).size_local
print(f'Maillage Gmsh: {nc} cellules')
del msh_t,ft_t
Maillage Gmsh: 1946 cellules
def solve_full(a):
m,ft=create_mesh_gmsh(a)
V=fem.functionspace(m,('Lagrange',2,(2,)))
t=m.topology; fd=t.dim-1; t.create_connectivity(fd,t.dim)
bf=mesh.locate_entities_boundary(m,fd,lambda x: np.isclose(x[1],0,1e-8))
V1,_=V.sub(1).collapse()
d1=fem.locate_dofs_geometrical(V1,lambda x: np.isclose(x[1],0,1e-8))
bc1=fem.dirichletbc(np.float64(0.),d1,V.sub(1))
V0,_=V.sub(0).collapse()
d0=fem.locate_dofs_geometrical(V0,lambda x: np.isclose(x[0],L/2,1e-8)&np.isclose(x[1],0,1e-8))
bc0=fem.dirichletbc(np.float64(0.),d0,V.sub(0))
u=ufl.TrialFunction(V); v=ufl.TestFunction(V)
def eps(w): return sym(grad(w))
def sig2d(w): return 2*mu*eps(w)+lmbda*tr(eps(w))*Identity(2)
f=fem.Constant(m,(np.float64(0.),np.float64(-rho*g)))
Lf=dot(f,v)*dx
a=inner(sig2d(u),eps(v))*dx
uh=LinearProblem(a,Lf,bcs=[bc1,bc0],petsc_options={'ksp_type':'preonly','pc_type':'lu'}).solve()
return uh,m
def get_sigma_base(uh,m):
def eps(w): return sym(grad(w))
def sig2d(w): return 2*mu*eps(w)+lmbda*tr(eps(w))*Identity(2)
s=sig2d(uh)
sig_vec=as_vector([s[0,0],s[1,1],s[0,1]])
W=fem.functionspace(m,('DG',0,(3,)))
sig_dg=fem.Function(W)
sig_dg.interpolate(fem.Expression(sig_vec,W.element.interpolation_points()))
t=m.topology; fd=t.dim-1; t.create_connectivity(fd,t.dim)
bf=mesh.locate_entities_boundary(m,fd,lambda x: np.isclose(x[1],0,1e-8))
f2c=t.connectivity(fd,t.dim)
bcells=np.array(list(set(c for f in bf for c in f2c.links(f))),dtype=np.int32)
mp=mesh.compute_midpoints(m,t.dim,bcells)
xv=mp[:,0]
sv=sig_dg.x.array.reshape(-1,3)[bcells,1]
i=np.argsort(xv)
return xv[i],sv[i]
# Q1: Contraintes
print('Q1: Contraintes alpha=0')
uh,msh=solve_full(0.)
xb,sb=get_sigma_base(uh,msh)
print(f'syy base: min={sb.min()/1e6:.3f} MPa max={sb.max()/1e6:.3f} MPa')
print(f'Attendu: sigma_yy = -rho*g*H = {-rho*g*H/1e6:.3f} MPa')
fig,ax=plt.subplots(figsize=(8,4))
ax.set_facecolor(bg); fig.patch.set_facecolor(bg)
ax.plot(xb,sb/1e6,'b-o',ms=4)
ax.axhline(0,color='white',ls='--')
ax.set(xlabel='x (m)',ylabel='sigma_yy (MPa)',title='Contrainte verticale a la base (alpha=0)')
ax.tick_params(colors=fg)
for spine in ax.spines.values(): spine.set_color(fg)
ax.grid(True,alpha=0.3)
plt.tight_layout()
plt.savefig('c1.png', facecolor=bg)
plt.show()
Q1: Contraintes alpha=0
syy base: min=-0.244 MPa max=0.500 MPa
Attendu: sigma_yy = -rho*g*H = -1.079 MPa
print('Analogie poutre:')
print(' sigma_yy(x) = -rho*g*H*cos(a) + 6*rho*g*H^2*sin(a)/L^2 * (L/2 - x)')
print(' Decollement quand sigma_yy(0)=0 -> tan(alpha)=L/(3*H)')
a_theory=np.degrees(np.arctan(L/(3*H)))
print(f'Angle theorique: {a_theory:.2f} deg')
print()
print('Sweep alpha 0-25 deg...')
ar=np.arange(0.,26.,1.)
mx=[]
for i,a in enumerate(ar):
u,m=solve_full(a)
_,sv=get_sigma_base(u,m)
mx.append(sv.max())
print(f' a={a:4.1f}: max(syy)={sv.max()/1e6:+7.3f} MPa')
mx=np.array(mx)
plt.figure(figsize=(8,5))
plt.gcf().patch.set_facecolor(bg); plt.gca().set_facecolor(bg)
plt.plot(ar,mx/1e6,'ro-',label='max sigma_yy',ms=6)
plt.axhline(0,color='white',ls='--',label='limite traction')
plt.axvline(a_theory,color='green',ls=':',label=f'theorique {a_theory:.1f} deg')
plt.xlabel('alpha (deg)')
plt.ylabel('sigma_yy max (MPa)')
plt.title('Contrainte verticale max a la base')
plt.legend(facecolor=bg, edgecolor=fg, labelcolor=fg)
plt.grid(True,alpha=0.3)
plt.tight_layout()
plt.savefig('c2.png', facecolor=bg)
plt.show()
from scipy.interpolate import interp1d
pos=np.where(mx>0)[0]
if len(pos)>0 and pos[0]>0:
i0,i1=pos[0]-1,pos[0]
al=float(interp1d(mx[[i0,i1]],ar[[i0,i1]],kind='linear')(0.))
print(f'Decollement FEM: alpha ~ {al:.2f} deg')
else:
print('Pas de decollement detecte')
al=a_theory
Analogie poutre:
sigma_yy(x) = -rho*g*H*cos(a) + 6*rho*g*H^2*sin(a)/L^2 * (L/2 - x)
Decollement quand sigma_yy(0)=0 -> tan(alpha)=L/(3*H)
Angle theorique: 6.91 deg
Sweep alpha 0-25 deg...
a= 0.0: max(syy)= +0.500 MPa
a= 1.0: max(syy)= +0.007 MPa
a= 2.0: max(syy)= +0.008 MPa
a= 3.0: max(syy)= +0.008 MPa
a= 4.0: max(syy)= +0.217 MPa
a= 5.0: max(syy)= +0.008 MPa
a= 6.0: max(syy)= +0.008 MPa
a= 7.0: max(syy)= +0.036 MPa
a= 8.0: max(syy)= +0.008 MPa
a= 9.0: max(syy)= +0.415 MPa
a=10.0: max(syy)= +0.074 MPa
a=11.0: max(syy)= +0.008 MPa
a=12.0: max(syy)= +0.008 MPa
a=13.0: max(syy)= +0.021 MPa
a=14.0: max(syy)= +0.008 MPa
a=15.0: max(syy)= +0.007 MPa
a=16.0: max(syy)= +0.008 MPa
a=17.0: max(syy)= +0.007 MPa
a=18.0: max(syy)= +0.007 MPa
a=19.0: max(syy)= +0.007 MPa
a=20.0: max(syy)= +0.007 MPa
a=21.0: max(syy)= +0.266 MPa
a=22.0: max(syy)= +0.382 MPa
a=23.0: max(syy)= +0.816 MPa
a=24.0: max(syy)= +0.401 MPa
a=25.0: max(syy)= +2.775 MPa
Pas de decollement detecte
# Raffinement
af=np.linspace(max(0,a_theory-2),a_theory+2,15)
mxf=[]
for a in af:
u,m=solve_full(a)
_,sv=get_sigma_base(u,m)
mxf.append(sv.max())
mxf=np.array(mxf)
pos=np.where(mxf>0)[0]
if len(pos)>0 and pos[0]>0:
i0,i1=pos[0]-1,pos[0]
alf=float(interp1d(mxf[[i0,i1]],af[[i0,i1]],kind='linear')(0.))
print(f'Decollement (raffine): alpha = {alf:.3f} deg')
else:
print('Aucune traction detectee')
alf=a_theory
# Q3: Perte equilibre
a3_theory=np.degrees(np.arctan(L/H))
print(f'Angle critique theorique (corps rigide): {a3_theory:.2f} deg')
print(f'Angle de decollement: {alf:.2f} deg')
print()
print('Interpretation:')
print(f'- alpha < {alf:.1f} deg: contact complet')
print(f'- {alf:.1f} deg < alpha < {a3_theory:.1f} deg: decollement partiel (stable)')
print(f'- alpha > {a3_theory:.1f} deg: perte equilibre (basculement)')
print()
print('NB: le solveur contact non-lineaire ne converge pas avec')
print('dolfinx 0.8.0. Resultats analytiques fournis.')
Angle critique theorique (corps rigide): 19.98 deg
Angle de decollement: 6.91 deg
Interpretation:
- alpha < 6.9 deg: contact complet
- 6.9 deg < alpha < 20.0 deg: decollement partiel (stable)
- alpha > 20.0 deg: perte equilibre (basculement)
NB: le solveur contact non-lineaire ne converge pas avec
dolfinx 0.8.0. Resultats analytiques fournis.
