В настоящее время я пишу код Matlab для решения проблемы теплопроводности на me sh, созданной с помощью программного обеспечения GM SH. В этом коде я использую метод с главным элементом для вычисления функций формы в каждом узле каждого четырехугольного элемента и после этого элементарных матриц проводимости. После того, как я вычислил элементарные матрицы проводимости, когда я проверил их собственные значения, я обнаружил, что одна из них равна нулю, что означает, что эти матрицы являются единичными, какими они не должны быть (мой профессор сказал мне это). Поэтому, возможно, моя функция для вычисления матрицы проводимости элементов неверна. Может кто-нибудь помочь мне заметить, что не так с моей функцией?
function [k_e]= compute_element_conductivity_matrix(x1,x2,x3,x4,y1,y2,y3,y4,kcond)
% Gauss weights
W=1; %W1=W2=W3=W4=W=1
% Number of integration points
n_int=4;
nodes=4;
% Integration points
ksi=[-1/sqrt(3) 1/sqrt(3) 1/sqrt(3) -1/sqrt(3)];
eta=[-1/sqrt(3) -1/sqrt(3) 1/sqrt(3) 1/sqrt(3)];
% Global coordinates
x=[x1;x2;x3;x4];
y=[y1;y2;y3;y4];
for (i =1:n_int)
% Shape functions
N1(i)=(1/4)*(1-ksi(i))*(1-eta(i));
N2(i)=(1/4)*(1+ksi(i))*(1-eta(i));
N3(i)=(1/4)*(1+ksi(i))*(1+eta(i));
N4(i)=(1/4)*(1-ksi(i))*(1+eta(i));
% Shape functions derivatives wrt eta
N1_eta(i)= (1/4)*(ksi(i)-1);
N2_eta(i)= (1/4)*(-ksi(i)-1);
N3_eta(i)= (1/4)*(ksi(i)+1);
N4_eta(i)= (1/4)*(1-ksi(i));
% Shape functions derivatives wrt ksi
N1_ksi(i)= (1/4)*(eta(i)-1);
N2_ksi(i)= (1/4)*(1-eta(i));
N3_ksi(i)= (1/4)*(eta(i)+1);
N4_ksi(i)= (1/4)*(-eta(i)-1);
% Computation of the Global derivatives of the shape functions and the Jacobian determinants
x_ksi(i)=N1_ksi(i)*x1+N2_ksi(i)*x2+N3_ksi(i)*x3+N4_ksi(i)*x4;
y_ksi(i)=N1_ksi(i)*y1+N2_ksi(i)*y2+N3_ksi(i)*y3+N4_ksi(i)*y4;
x_eta(i)=N1_eta(i)*x1+N2_eta(i)*x2+N3_eta(i)*x3+N4_eta(i)*x4;
y_eta(i)=N1_eta(i)*y1+N2_eta(i)*y2+N3_eta(i)*y3+N4_eta(i)*y4;
j(i)=x_ksi(i)*y_eta(i)-x_eta(i)*y_ksi(i);
j_inv(i)=inv(j(i));
N1_x(i)=((N1_ksi(i)*y_eta(i)-N1_eta(i)*y_ksi(i)))*j_inv(i);
N1_y(i)=(-(N1_ksi(i)*x_eta(i)-N1_eta(i)*x_ksi(i)))*j_inv(i);
N2_x(i)=((N2_ksi(i)*y_eta(i)-N2_eta(i)*y_ksi(i)))*j_inv(i);
N2_y(i)=(-(N2_ksi(i)*x_eta(i)-N2_eta(i)*x_ksi(i)))*j_inv(i);
N3_x(i)=((N3_ksi(i)*y_eta(i)-N3_eta(i)*y_ksi(i)))*j_inv(i);
N3_y(i)=(-(N3_ksi(i)*x_eta(i)-N3_eta(i)*x_ksi(i)))*j_inv(i);
N4_x(i)=((N4_ksi(i)*y_eta(i)-N4_eta(i)*y_ksi(i)))*j_inv(i);
N4_y(i)=(-(N4_ksi(i)*x_eta(i)-N4_eta(i)*x_ksi(i)))*j_inv(i);
% grad_N1=[N1_x ; N1_y];
% grad_N2=[N2_x ; N2_y];
% grad_N3=[N3_x ; N3_y];
% grad_N4=[N4_x ; N4_y];
%Element conductivity matrix
k_11(i)=(N1_x(i)*N1_x(i)+N1_y(i)*N1_y(i))*kcond;
k_12(i)=(N1_x(i)*N2_x(i)+N1_y(i)*N2_y(i))*kcond;
k_13(i)=(N1_x(i)*N3_x(i)+N1_y(i)*N3_y(i))*kcond;
k_14(i)=(N1_x(i)*N4_x(i)+N1_y(i)*N4_y(i))*kcond;
k_22(i)=(N2_x(i)*N2_x(i)+N2_y(i)*N2_y(i))*kcond;
k_23(i)=(N2_x(i)*N3_x(i)+N2_y(i)*N3_y(i))*kcond;
k_24(i)=(N2_x(i)*N4_x(i)+N2_y(i)*N4_y(i))*kcond;
k_33(i)=(N3_x(i)*N3_x(i)+N3_y(i)*N3_y(i))*kcond;
k_34(i)=(N3_x(i)*N4_x(i)+N3_y(i)*N4_y(i))*kcond;
k_44(i)=(N4_x(i)*N4_x(i)+N4_y(i)*N4_y(i))*kcond;
end
%Contibution of each integration points to the Element conductivity matrix
k_e=zeros(4);
k_e(1,1)=k_11(1)*W+k_11(2)*W+k_11(3)*W+k_11(4)*W;
k_e(1,2)=k_12(1)*W+k_12(2)*W+k_12(3)*W+k_12(4)*W;
k_e(2,1)=k_e(1,2);
k_e(1,3)=k_13(1)*W+k_13(2)*W+k_13(3)*W+k_13(4)*W;
k_e(3,1)=k_e(1,3);
k_e(1,4)=k_14(1)*W+k_14(2)*W+k_14(3)*W+k_14(4)*W;
k_e(4,1)=k_e(1,4);
k_e(2,2)=k_22(1)*W+k_22(2)*W+k_22(3)*W+k_22(4)*W;
k_e(2,3)=k_23(1)*W+k_23(2)*W+k_23(3)*W+k_23(4)*W;
k_e(3,2)=k_e(2,3);
k_e(2,4)=k_24(1)*W+k_24(2)*W+k_24(3)*W+k_24(4)*W;
k_e(4,2)=k_e(2,4);
k_e(3,3)=k_33(1)*W+k_33(2)*W+k_33(3)*W+k_33(4)*W;
k_e(3,4)=k_34(1)*W+k_34(2)*W+k_34(3)*W+k_34(4)*W;
k_e(4,3)=k_e(3,4);
k_e(4,4)=k_44(1)*W+k_44(2)*W+k_44(3)*W+k_44(4)*W;
k_e;
end