To build assembly matrix from elemental matrix
조회 수: 2 (최근 30일)
이전 댓글 표시
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
Walter Roberson
2022년 8월 6일
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
2022년 8월 6일
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.
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Numbers and Precision에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!