clc; clear; close all;
%% Paramètres physiques
E = 2.3e9; % Module d'Young [Pa]
nu = 0.2; % Coefficient de Poisson
gamma_c = 27e3;% Poids volumique du béton [N/m³]
h_totale = 20; % Hauteur totale [m]
L = 1; % Épaisseur du mur [m]
el_size = 0.5; % Taille élément [m]
h = 4; % Largeur du mur [m]
%% Maillage T6
% Création du maillage triangulaire rectangulaire via Delaunay
x = 0:el_size:L;
y = 0:el_size:h_totale;
[X, Y] = meshgrid(x, y);
DT = delaunayTriangulation(X(:), Y(:));
nodes = [Link];
elements_T3 = [Link];
nn = size(nodes, 1);
nel = size(elements_T3, 1);
middle_nodes = [Link];
T6_elements = zeros(nel, 6);
for e = 1:nel
n = elements_T3(e, :);
T6_elements(e, 1:3) = n;
edges = [n(1), n(2); n(2), n(3); n(3), n(1)];
for i = 1:3
edge = sort(edges(i, :));
key = sprintf('%d-%d', edge(1), edge(2));
if ~isKey(middle_nodes, key)
coord = (nodes(edge(1), :) + nodes(edge(2), :)) / 2;
nodes = [nodes; coord];
new_node = size(nodes, 1);
middle_nodes(key) = new_node;
end
T6_elements(e, 3+i) = middle_nodes(key);
end
end
nn = size(nodes, 1);
elements = T6_elements;
%% Visualisation du maillage T6
figure;
subplot(1,2,1);
triplot(elements(:, 1:3), nodes(:, 1), nodes(:, 2), 'b');
axis([-2,5,0,20]);
subplot(1,2,2);
hold on;
triplot(elements(:, 1:3), nodes(:, 1), nodes(:, 2), 'b');
plot(nodes(:, 1), nodes(:, 2), 'r.');
title('Maillage T6 (avec milieux d''arêtes)');
axis equal;
%% Numérotation des nœuds
for i = 1:size(nodes,1)
text(nodes(i,1), nodes(i,2), sprintf('%d', i), 'FontSize', 6,
'Color', 'b');
end
% Numérotation des éléments
for e = 1:size(elements,1)
coords = nodes(elements(e,1:3), :); % Les 3 premiers nœuds sont les
sommets
centroid = mean(coords);
text(centroid(1), centroid(2), sprintf('%d', e), 'FontSize', 6,
'Color', 'r');
end
% Matrice constitutive D en contrainte plane
D = (E/(1-nu^2))*[1 nu 0; nu 1 0; 0 0 (1-nu)/2];
%% Assemblage global
K = sparse(2*nn, 2*nn);
F = zeros(2*nn,1);
% Points de Gauss triangle (ordre 2)
gp = [1/6 1/6; 2/3 1/6; 1/6 2/3];
w = [1/6 1/6 1/6];
for e = 1:nel
coords = nodes(elements(e,:),:);
Ke = zeros(12);
fe = zeros(12,1);
for i = 1:3
[N, dNdxi] = shapeT6(gp(i,1), gp(i,2));
J = coords' * dNdxi;
detJ = det(J);
dNdx = dNdxi / J;
B = computeB(dNdx);
Ke = Ke + (B' * D * B) * detJ * w(i) * h;
% CORRECTION: Calcul correct du vecteur forces
N_matrix = [N(1) 0 N(2) 0 N(3) 0 N(4) 0 N(5) 0 N(6) 0;
0 N(1) 0 N(2) 0 N(3) 0 N(4) 0 N(5) 0 N(6)];
body_force = [0; -gamma_c]; % Force volumique [0; -γ]
fe = fe + N_matrix' * body_force * detJ * w(i) * h;
end
% Connectivité DOFs
edof = [2*elements(e,:)-1; 2*elements(e,:)];
edof = edof(:)';
K(edof,edof) = K(edof,edof) + Ke;
F(edof) = F(edof) + fe;
end
%% Application des charges latérales
tol = 1e-6;
sides = [0, L];
for s = 1:2
side = sides(s);
edge_nodes = find(abs(nodes(:,1) - side) < tol);
[~, idx] = sort(nodes(edge_nodes,2));
sorted_nodes = edge_nodes(idx);
for k = 1:length(sorted_nodes)-1
n1 = sorted_nodes(k);
n2 = sorted_nodes(k+1);
y1 = nodes(n1,2); y2 = nodes(n2,2);
Le = abs(y2 - y1);
if side == 0
p1 = (200 - 10*y1) * 1e3;
p2 = (200 - 10*y2) * 1e3;
else
p1 = max(0, 100e3 * (-1 + y1/10));
p2 = max(0, 100e3 * (-1 + y2/10));
end
% Distribution sur les noeuds d'extrémité seulement
f1 = h * Le * (2*p1 + p2)/6;
f2 = h * Le * (p1 + 2*p2)/6;
F(2*n1-1) = F(2*n1-1) + f1;
F(2*n2-1) = F(2*n2-1) + f2;
end
end
%% Conditions aux limites
fixed_nodes = find(abs(nodes(:,2)) < 1e-6);
fixed_dofs = sort([2*fixed_nodes-1; 2*fixed_nodes]);
free_dofs = setdiff(1:2*nn, fixed_dofs);
%% Résolution
U = zeros(2*nn,1);
U(free_dofs) = K(free_dofs,free_dofs) \ F(free_dofs)*7e-2;
%% Visualisation
plotResults(nodes, elements, U, 1);
%% Post traitement
% Calcul des contraintes et déformations aux points de Gauss
strain = zeros(nel, 3);
stress = zeros(nel, 3);
centroids = zeros(nel, 2);
for e = 1:nel
nodes_e = elements(e,:); % les 6 nœuds de l'élément
coords = nodes(nodes_e,:); % coordonnées 6x2
Ue = reshape(U([2*nodes_e-1; 2*nodes_e]), 2, [])'; % 6x2 → reshape
en 12x1
% Points de Gauss pour intégration
gp = [1/6 1/6; 2/3 1/6; 1/6 2/3];
% Initialisation
strain_e = zeros(1,3);
stress_e = zeros(1,3);
for i = 1:3
[~, dNdxi] = shapeT6(gp(i,1), gp(i,2));
J = coords' * dNdxi;
dNdx = dNdxi / J;
B = computeB(dNdx);
strain_gp = B * Ue(:);
stress_gp = D * strain_gp;
strain_e = strain_e + strain_gp';
stress_e = stress_e + stress_gp';
end
strain(e,:) = strain_e / 3; % Moyenne sur les points de Gauss
stress(e,:) = stress_e / 3;
centroids(e,:) = mean(coords(1:3,:)); % Centre du triangle (sommets
seulement)
end
%% Interpolation nodale améliorée
stress_nodal = zeros(nn,3);
strain_nodal = zeros(nn,3);
count = zeros(nn,1);
for e = 1:nel
% Poids selon la proximité aux nœuds
for k = 1:6
n = elements(e,k);
stress_nodal(n,:) = stress_nodal(n,:) + stress(e,:);
strain_nodal(n,:) = strain_nodal(n,:) + strain(e,:);
count(n) = count(n) + 1;
end
end
for i = 1:nn
if count(i) > 0
stress_nodal(i,:) = stress_nodal(i,:) / count(i);
strain_nodal(i,:) = strain_nodal(i,:) / count(i);
end
end
%% Visualisation avec triangulation T3
tri = elements(:,1:3);
% Figure 1: Contraintes
figure('Name','Contraintes','Position',[100 100 1200 800])
subplot(1,3,1)
trisurf(tri, nodes(:,1), nodes(:,2), zeros(nn,1), stress_nodal(:,1)/1e6,
...
'FaceColor', 'interp', 'EdgeColor', 'none');
view(2), axis equal, title('\sigma_{xx} [MPa]'), colorbar
subplot(1,3,2)
trisurf(tri, nodes(:,1), nodes(:,2), zeros(nn,1), stress_nodal(:,2)/1e6,
...
'FaceColor', 'interp', 'EdgeColor', 'none');
view(2), axis equal, title('\sigma_{yy} [MPa]'), colorbar
subplot(1,3,3)
trisurf(tri, nodes(:,1), nodes(:,2), zeros(nn,1), stress_nodal(:,3)/1e6,
...
'FaceColor', 'interp', 'EdgeColor', 'none');
view(2), axis equal, title('\tau_{xy} [MPa]'), colorbar
% Figure 2: Déformations
figure('Name','Déformations','Position',[100 100 1200 800])
subplot(1,3,1)
trisurf(tri, nodes(:,1), nodes(:,2), zeros(nn,1), strain_nodal(:,1)*100,
...
'FaceColor', 'interp', 'EdgeColor', 'none');
view(2), axis equal, title('\epsilon_{xx} [%]'), colorbar
subplot(1,3,2)
trisurf(tri, nodes(:,1), nodes(:,2), zeros(nn,1), strain_nodal(:,2)*100,
...
'FaceColor', 'interp', 'EdgeColor', 'none');
view(2), axis equal, title('\epsilon_{yy} [%]'), colorbar
subplot(1,3,3)
trisurf(tri, nodes(:,1), nodes(:,2), zeros(nn,1), strain_nodal(:,3)*100,
...
'FaceColor', 'interp', 'EdgeColor', 'none');
view(2), axis equal, title('\gamma_{xy} [%]'), colorbar
%% Export des résultats vers Excel
filename = 'Resultats_MEF_T6.xlsx';
% Table des noeuds
T_noeuds = array2table([(1:nn)', nodes], ...
'VariableNames', {'Noeud_ID', 'X', 'Y'});
% Table de connectivité (160 lignes, 6 nœuds par élément)
T_connect = array2table([(1:nel)', elements], ...
'VariableNames', ['Element_ID', arrayfun(@(x) sprintf('Noeud_%d',x),
1:6, 'UniformOutput', false)]);
% Table des déplacements (405 lignes, 1 par nœud)
T_depl = array2table([(1:nn)', U([Link]nd), U([Link]nd)], ...
'VariableNames', {'Noeud_ID', 'Ux', 'Uy'});
% Tables contraintes/déformations (160 lignes)
T_contraintes = array2table([(1:nel)', stress/1e6], ...
'VariableNames', {'Element_ID', 'Sigma_xx_MPa', 'Sigma_yy_MPa',
'Tau_xy_MPa'});
T_deformations = array2table([(1:nel)', strain*100], ...
'VariableNames', {'Element_ID', 'Epsilon_xx_pourcent',
'Epsilon_yy_pourcent', 'Gamma_xy_pourcent'});
% Export vers Excel
writetable(T_noeuds, filename, 'Sheet', 'Noeuds');
writetable(T_connect, filename, 'Sheet', 'Connectivite');
writetable(T_depl, filename, 'Sheet', 'Deplacements');
writetable(T_contraintes, filename, 'Sheet', 'Contraintes');
writetable(T_deformations, filename, 'Sheet', 'Deformations');
% Affichage infos console
Ux = U([Link]nd);
Uy = U([Link]nd);
disp(['Résultats exportés vers ' filename]);
disp(['Déplacement maximal suivant Ox ' num2str(100*max(Ux)) ' cm']);
disp(['Déplacement maximal suivant Oy ' num2str(100*max(Uy)) ' cm']);
function [N, dNdxi] = shapeT6(L1, L2)
L3 = 1 - L1 - L2;
N = [L1*(2*L1-1); L2*(2*L2-1); L3*(2*L3-1); 4*L1*L2; 4*L2*L3;
4*L3*L1];
dNdxi = [4*L1-1, 0; 0, 4*L2-1; -4*L3+1, -4*L3+1; 4*L2, 4*L1; -4*L2,
4*(1-2*L2-L1); 4*(1-2*L1-L2), -4*L1];
end
function B = computeB(dNdx)
B = zeros(3,12);
for i = 1:6
B(1,2*i-1) = dNdx(i,1);
B(2,2*i) = dNdx(i,2);
B(3,2*i-1) = dNdx(i,2);
B(3,2*i) = dNdx(i,1);
end
end
function plotResults(nodes, elements, U, scale_factor)
Ux = U([Link]nd) * scale_factor;
Uy = U([Link]nd) * scale_factor;
displaced = nodes + [Ux Uy];
figure;
hold on;
% Tracé initial (en bleu)
for e = 1:size(elements,1)
plot(nodes(elements(e,[1 4 2 5 3 6 1]),1), nodes(elements(e,[1 4
2 5 3 6 1]),2), 'b-');
end
% Tracé déformé (en rouge)
for e = 1:size(elements,1)
plot(displaced(elements(e,[1 4 2 5 3 6 1]),1),
displaced(elements(e,[1 4 2 5 3 6 1]),2), 'r-');
end
axis equal;
title(sprintf('Déformation (échelle x%.1f)',scale_factor));
legend('Initial', 'Déformé');
xlabel('x [m]'); ylabel('y [m]');
end
