To build assembly matrix from elemental matrix
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
here is the code below. i tried to assemble elemental matrix using a for loop but the result ends with a matrix having a determinant 0. i need to use iteration to solve the resulting equation but the convergence of the matrix doesn't allow. i dont really know where i make the mistake whether if it is on the code or somewhere else. can i get any default assemblying code or can anyone help to assess my code and idetify where i made a mistake. the analysis is based on finite element method.
%% global matrix
clear all
clc
G=[112 28 29;
153 107 113;
% 150 110 108;
% 174 142 149;
138 89 91;
% 172 170 165;
63 161 168;
% 157 122 158;
% 175 115 154;
66 67 33]
%% stiffness matrix
syms z e
mu=0.28
density=0.05
for i=1:3
if i==1
s(1)=1-z-e;
elseif i==2
s(2)=z;
else
s(3)=e;
end
end
s(1)
s(2)
s(3)
%% jacobian matrix
jj=[diff(s(1),z) diff(s(2),z) diff(s(3),z);
diff(s(1),e) diff(s(2),e) diff(s(3),e)]
nodeco=[32.92 17.76;30.76 20.4;29.2 16.77]
jacobian=jj*nodeco
invJ=inv(jacobian)
D=det(jacobian)
%% stiffness matrix of velocity
for i=1:3
for j=1:3
ku(i,j)=int(int((invJ(2,1)*diff(s(i),z)+invJ(2,2)*diff(s(j),e))*(invJ(2,1)*diff(s(i),z)+invJ(2,2)*diff(s(j),e)),e,0,1-z),z,0,1)*D*mu;
end
end
ku
%% stiffness for pressure
for i=1:3
for j=1:3
kp(i,j)=int(int(s(i)*(invJ(1,1)*diff(s(j),z)+invJ(1,2)*diff(s(j),e)+invJ(2,1)*diff(s(j),z)+invJ(2,2)*diff(s(j),e)),e,0,1-z),z,0,1)*D/density;
end
end
kp
%%assembled velocity coefficient
sum=zeros(168);
for i=1:5
K=zeros(168);
for j=1:3
% [M(j,1),M(j,2), M(j,3)]==[K(G(i,j),G(i,1)) K(G(i,j),G(i,2)) K(G(i,j),G(i,3))]
K(G(i,j),G(i,1))=ku(j,1);
K(G(i,j),G(i,2))=ku(j,2);
K(G(i,j),G(i,3))=ku(j,3);
end
sum=sum+K;
end
kuu=sum
kuuu=kuu([28 29 33 63 66 67 89 91 107 112 113 138 153 161 168],[28 29 33 63 66 67 89 91 107 112 113 138 153 161 168])
%%assembled pressure coefficient
sum=zeros(168);
for i=1:5
K=zeros(168);
for j=1:3
% [M(j,1),M(j,2), M(j,3)]==[K(G(i,j),G(i,1)) K(G(i,j),G(i,2)) K(G(i,j),G(i,3))]
K(G(i,j),G(i,1))=kp(j,1);
K(G(i,j),G(i,2))=kp(j,2);
K(G(i,j),G(i,3))=kp(j,3);
end
sum=sum+K;
end
kpp=sum
kppp=kpp([28 29 33 63 66 67 89 91 107 112 113 138 153 161 168],[28 29 33 63 66 67 89 91 107 112 113 138 153 161 168])
%% source term
syms B Ui double
B=20;
conductivity=100;
for i=1:3
f(i,1)=int(int((s(i)*(B^2)*Ui),e,0,1-z),z,0,1)*(conductivity*D)/density;
end
f
%% assemble of source term
sum=zeros([168,1]);
for i=1:5
K=sym(zeros([168,1])) ;
for j=1:3
K(G(i,j),1)=f(j,1);
end
sum=sum+K;
end
kff=sum;
kfff=kff([28 29 33 63 66 67 89 91 107 112 113 138 153 161 168],1)
% double(kfff)
2 commentaires
Walter Roberson
le 6 Août 2022
syms B Ui double
Caution: that would create a symbolic variable named double . There is no double modifier for symbolic variables. There are modifiers such as real and positive
Walter Roberson
le 6 Août 2022
sum=zeros([168,1]);
It is recommended to avoid using sum as the name of a variable, as it conflicts with the MATLAB function sum() . It is very common for people to use sum as a variable name and then a few lines further down try to use sum() as a function.
Réponses (0)
Voir également
Catégories
En savoir plus sur Numbers and Precision dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!