Матрица податливости

Есть КЭ-модель. Синим показаны закрепления в узлах, фиолетовым master-узлы. 

С помощью метода подструктур можно выписать матрицу податливости выбранных узлов (показано фиолетовым), которая описывает податливость всей системы в выбранном направлении. В теории, если умножить эту матрицу с вектором сил, то получим вектор деформаций в выбранных master-узлах.

Выписал матрицы податливости, координаты узлов и их эффективные площади для разных сеток (см. files.rar):

1. линейные элементы 10х10
2. линейные элементы 20х20
3. квадратичные элементы 10х10
4. квадратичные элементы 20х20

В Matlab перемножаю матрицы податливости с соответствующими векторами сил, не забывая учитывать весовые коэффициенты и эффективные площади узлов

clc; close all; clear
 
[Ulin10, philin10, Zlin10] = calc_disp('C_mat_lin10.txt','Node_info1_lin10.txt');
figure
contourf(philin10, Zlin10,Ulin10,40,'EdgeColor','none')
colorbar; colormap jet; ylabel('Y [m]'); xlabel('\phi [°]'); title('linear 10x10')
 
[Ulin20, philin20, Zlin20] = calc_disp('C_mat_lin20.txt','Node_info1_lin20.txt');
figure
contourf(philin20, Zlin20,Ulin20,40,'EdgeColor','none')
colorbar; colormap jet; ylabel('Y [m]'); xlabel('\phi [°]'); title('linear 20x20')
 
[Uq10, phiq10, Zq10] = calc_disp('C_mat_q10.txt','Node_info1_q10.txt','Node_info2_q10.txt');
figure
contourf(phiq10, Zq10,Uq10,40,'EdgeColor','none')
colorbar; colormap jet; ylabel('Y [m]'); xlabel('\phi [°]'); title('quad 10x10')
 
[Uq20, phiq20, Zq20] = calc_disp('C_mat_q20.txt','Node_info1_q20.txt','Node_info2_q20.txt');
figure
contourf(phiq20, Zq20,Uq20,40,'EdgeColor','none')
colorbar; colormap jet; ylabel('Y [m]'); xlabel('\phi [°]'); title('quad 10x10')
 
figure
hold on; grid on
plot(philin10(21,:),Ulin10(21,:),'-ob','LineWidth',1.5,'DisplayName','linear10x10')
plot(philin20(21,:),Ulin20(21,:),'--sr','LineWidth',1.5,'DisplayName','linear20x20')
plot(phiq10(21,:),Uq10(21,:),'-dg','LineWidth',1.5,'DisplayName','quad10x10')
plot(phiq20(21,:),Uq20(21,:),'--^k','LineWidth',1.5,'DisplayName','quad20x20')
legend; xlabel('\phi [°]'); ylabel('U');
 
figure
hold on; grid on
plot(Zlin10(:,21),Ulin10(:,21),'-ob','LineWidth',1.5,'DisplayName','linear10x10')
plot(Zlin20(:,21),Ulin20(:,21),'--sr','LineWidth',1.5,'DisplayName','linear20x20')
plot(Zq10(:,21),Uq10(:,21),'-dg','LineWidth',1.5,'DisplayName','quad10x10')
plot(Zq20(:,21),Uq20(:,21),'--^k','LineWidth',1.5,'DisplayName','quad20x20')
legend; xlabel('Y [m]'); ylabel('U');
 
function [U_int, phi_int, Z_int] = calc_disp(c_mat_filename,nodeinfo1_filename,nodeinfo2_filename)
if nargin == 3
    
    N = readmatrix(c_mat_filename,'DecimalSeparator','.');
    [NID, phi,z,A] = readvars(nodeinfo1_filename,'NumHeaderLines',1,...
        'ExpectedNumVariables',4,'DecimalSeparator','.');
    [NIDeck, ~,~,~,~] = readvars(nodeinfo2_filename,'NumHeaderLines',1,...
        'ExpectedNumVariables',5,'DecimalSeparator','.');
    
    F = ones(length(phi),1);
    for i = 1:length(NIDeck)
        F(NID==NIDeck(i)) = -1/12;
    end
    
    F(F==1) = 1/3;
    F = F.*1000.*A;
    
elseif nargin == 2
    N = readmatrix(c_mat_filename,'DecimalSeparator','.');
    [phi,z,A] = readvars(nodeinfo1_filename,'NumHeaderLines',1,...
        'ExpectedNumVariables',3,'DecimalSeparator','.');
    
    F = ones(length(phi),1)*1e3/4.*A;
end
 
u = N*F;
 
[phi_int, Z_int] = meshgrid(linspace(min(phi),max(phi),41),...
    linspace(min(z),max(z),41));
U_int = griddata(phi,z,u,phi_int,Z_int,'linear');
end

Результат:

В итоге для линейных элементов перемещения выглядят более-менее правдоподобно (учитывая что сетка очень грубая). Но для квадратичных элементов перемещения получаются "зубастыми". И тут вопрос, то ли проблема при экспорте матрицы из ANSYS, то ли в последующей обработке в матлаб?

Результат для полной модели из ANSYS: