csgt0
精度有问题？没程序和数据很难说明问题
z770428
可是刚度矩阵中所有元素都不大也不太小，你是说MATLAB精度有问题么
zhuomec
刚度矩阵主对角元素应该是大于零的，你是不是算错了？
z770428
谢谢，肯定错了，能运行得到刚度矩阵，但是K(138，138)，K(139，139)=0，K(140,140)=0.其他对角线都大于0，为什么呢？
csgt0
引用回帖:
z770428 at 2012-09-17 12:35:37
谢谢，肯定错了，能运行得到刚度矩阵，但是K(138，138)，K(139，139)=0，K(140,140)=0.其他对角线都大于0，为什么呢？
上程序看看
z770428
引用回帖:
csgt0 at 2012-09-17 13:21:44
上程序看看...
function [globalKT] = conductMatrix(omega,globalTDOF,iter,enrElem)
% This function calculates the global stiffness matrix for the desired
% discontinuities defined by the user supplied input.
global CONNEC CRACK DOMAIN MAT NODES PHI PSI XYZ
DOMAIN(1) = 5;DOMAIN(2) = 20;DOMAIN(3) = 0.2;DOMAIN(4) = 0.2; MAT(1) = 5.5E10; MAT(2) = 2.1E10;MAT(3) = 9.7E9;MAT(4) = 3.46;MAT(5) = 0.35; MAT(6) = 1; MAT(8) = 6.3E-6;MAT(9) = 2.0E-5;MAT(10) = 0.25;
iter=1;globalTDOF=147;omega=pi/4;enrElem=[ 46 47 51 52 53 56 57 58 59 61 62 63 64 67 68 69];
nXElem = DOMAIN(1); % Number of elements in the x-direction
nYElem = DOMAIN(2); % Number of elements in the y-direction
lXElem = DOMAIN(3); % Length of elements in the x-direction
E1 = MAT(1); % Young's modulus for the matrix
E2 = MAT(2); % Poisson's ratio for the matrix
G12 = MAT(3); % Young's modulus for the fiber
k11 = MAT(4); % Poisson's ratio for the fiber
k22 = MAT(5); % Plane stress or plane strain
nuxy12 = MAT(10);
plane = MAT(11);
nCT = size(PHI,2);
% Number of crack tips
% Initialize the FE stiffness matrix
if iter > 1, globalTDOF = globalTDOF+16*iter; end % Initialize extra space for growing
globalKT = sparse(globalTDOF,globalTDOF); % Define the global K
% Create thermal conductivity constant matrix
if plane == 1; % Plane stress
h = MAT(6); % Plane stress thickness
c1 = k11; % Constant for elastic constant matrix
c2 = k22; % Constant for elastic constant matrix
C = h*[c1 0;...
0 c2 ];
end
% Stiffness of matrix
m = size(CRACK,1); % Determine number of data points defining crack
if m > 0
if nCT == 1
xCT = CRACK(m,1); % X-coordinate of crack tip
yCT = CRACK(m,2); % Y-coordinate of crack tip
elseif nCT == 2
xCT = [CRACK(1,1) CRACK(m,1)]; % X-coordinates of crack tips
yCT = [CRACK(1,2) CRACK(m,2)]; % Y-coordinates of crack tips
end
end
% Initialize the variables which will create the traditional sparse matrix
nIndexT = 0; % Initialize index
uenrElem = nXElem*nYElem-length(enrElem); % Number of unenriched elements
allRowT = ones(uenrElem*32,1); % Row indices
allColT = ones(uenrElem*32,1); % Column indices
allValT = zeros(uenrElem*32,1); % Stiffness matrix values
for iElem = 1nYElem*nXElem)
N1 = CONNEC(iElem,2); % Node 1 for current element
N2 = CONNEC(iElem,3); % Node 2 for current element
N3 = CONNEC(iElem,4); % Node 3 for current element
N4 = CONNEC(iElem,5); % Node 4 for current element
NN = NODES([N1 N2 N3 N4]',; % Nodal data for current element
CTN = nnz(NN(:,4)); % Number of nodes with crack tip enrichment
HEN = nnz(NN(:,2)); % Number of nodes with Heaviside enrichment % Number of inclusion nodes
NEN = HEN+CTN; % Number of enriched nodes
localKT = 0; % Initialize stiffness for current element
local = [N1 N2 N3 N4]; % Traditional index locations
iLoc = 5; % Next index location
if (NEN == 0) % Unenriched nodes
% Traditional element
X1 = XYZ(N1,2); X2 = XYZ(N2,2); X3 = XYZ(N3,2); X4 = XYZ(N4,2); % Nodal x-coordinates
Y1 = XYZ(N1,3); Y2 = XYZ(N2,3); Y3 = XYZ(N3,3); Y4 = XYZ(N4,3); % Nodal y-coordinates
xyz = [X1 Y1;X2 Y2;X3 Y3;X4 Y4]; % Nodal coordinate matrix
[gp,gw,J] = subDomain(3,[],xyz,0,0,[]);
xi = gp(:,1);
eta = gp(:,2);
N = 1/4*[(1-xi).*(1-eta) (1+xi).*(1-eta) (1+xi).*(1+eta) (1-xi).*(1+eta)];
for i = 1:length(gp)
xi = gp(i,1); eta = gp(i,2); % Gauss points
W = gw(i); % Gauss weights
Ji = [J(i,1) J(i,2);J(i,3) J(i,4)]; % Jacobian of subdomain
detJ = det(Ji); % Determinant of the Jacobian
Nx = 2/lXElem*1/4*[-(1-eta);1-eta;1+eta;-(1+eta)]; % Derivative of shape functions with respect to x
Ny = 2/lXElem*1/4*[-(1-xi);-(1+xi);1+xi;1-xi]; % Derivative of shape functions with respect to y
Bu = [Nx(1) Nx(2) Nx(3) Nx(4);Ny(1) Ny(2) Ny(3) Ny(4)];
localKT = localKT + W*Bu'*C*Bu*detJ;
end
elseif NEN > 0; % Enriched element
X1 = XYZ(N1,2); X2 = XYZ(N2,2); X3 = XYZ(N3,2); X4 = XYZ(N4,2); % Nodal x-coordinates
Y1 = XYZ(N1,3); Y2 = XYZ(N2,3); Y3 = XYZ(N3,3); Y4 = XYZ(N4,3); % Nodal y-coordinates
if NEN == 4 % Fully enriched element
xyz = [X1 Y1;X2 Y2;X3 Y3;X4 Y4]; % Nodal coordinate matrix
if numel(PSI) == 0, PN = [0 0 0 0]; else
PN = [ PSI(N1) PSI(N2) PSI(N3) PSI(N4)]; % Nodal crack level set values
end
if HEN == 4 % Full Heaviside enrichment
[gp,gw,J] = subDomain(3,PN,xyz,0,0,[]);
elseif CTN == 4 % Full crack tip enrichmnet
[gp,gw,J] = subDomain(7,PN,xyz,1,0,[]);
else % Full heaviside/crack tip enrichment
[gp,gw,J] = subDomain(7,PN,xyz,0,0,[]);
end
% Partially enriched element
else
PN = [ PSI(N1) PSI(N2) PSI(N3) PSI(N4)]; % Nodal crack level set values
[gp,gw,J] = subDomain(3,PN,xyz,0,0,[]);
xi = gp(:,1);
eta = gp(:,2);
N = 1/4*[(1-xi).*(1-eta) (1+xi).*(1-eta) (1+xi).*(1+eta) (1-xi).*(1+eta)];
[gp,gw] = gauss(6,'QUAD');
J = [];
end
for i = 1:length(gp)
xi = gp(i,1); eta = gp(i,2); % Gauss points
W = gw(i); % Gauss weights
if isempty(J) == 0
Ji = [J(i,1) J(i,2);J(i,3) J(i,4)]; % Jacobian of subdomain
detJ = det(Ji); % Determinant of the Jacobian
else
detJ = lXElem/2*lXElem/2; % Determinant of the Jacobian
end
N = 1/4*[(1-xi)*(1-eta);(1+xi)*(1-eta);... % Shape functions
(1+xi)*(1+eta);(1-xi)*(1+eta)];
Nx = 2/lXElem*1/4*[-(1-eta);1-eta;1+eta;-(1+eta)]; % Derivative of shape functions with respect to x
Ny = 2/lXElem*1/4*[-(1-xi);-(1+xi);1+xi;1-xi]; % Derivative of shape functions with respect to y
X1 = XYZ(N1,2); X2 = XYZ(N2,2); X3 = XYZ(N3,2); X4 = XYZ(N4,2); % Nodal x-coordinates
Y1 = XYZ(N1,3); Y2 = XYZ(N2,3); Y3 = XYZ(N3,3); Y4 = XYZ(N4,3); % Nodal y-coordinates
Xgp = N(1)*X1+N(2)*X2+N(3)*X3+N(4)*X4; % The global X for the current gauss point
Ygp = N(1)*Y1+N(2)*Y2+N(3)*Y3+N(4)*Y4; % The global Y for the current gauss point
Benr = [];
Bt = [Nx(1) Nx(2) Nx(3) Nx(4) ;... % 温度场一样，名称改为Bt
Ny(1) Ny(2) Ny(3) Ny(4)];
index = 1;
for iN = 1:4
if NN(iN,2) ~= 0
psi1 = PSI(N1); % Psi level set value at node 1
psi2 = PSI(N2); % Psi level set value at node 2
psi3 = PSI(N3); % Psi level set value at node 3
psi4 = PSI(N4); % Psi level set value at node 4
psi = N(1)*psi1+N(2)*psi2+N(3)*psi3+N(4)*psi4; % Psi level set value at current gauss point
Hgp = sign(psi); % Heaviside value at current gauss point
Hi = NN(iN,3); % Nodal Heaviside value
H = Hgp-Hi; % Shifted Heaviside value
Ba = [Nx(iN)*H ;... % 温度场分析同
Ny(iN)*H];
% Benr(:,index) = Ba; index = index+1; Benr =zeros(1000, 2000); Benr(:,indexindex+1)) = Ba; Benr =[]; index = index+1;
Benr = Ba;
if (i == length(gp))
local(iLoc) = (NN(iN,2));
iLoc = iLoc+1;
end
elseif NN(iN,4) ~= 0
if nCT == 1
X = Xgp-xCT; % Horizontal distance from crack tip to gauss point
Y = Ygp-yCT; % Vertical distance from crack tip to gauss point
CCS = [cos(omega) sin(omega);-sin(omega) cos(omega)];
XYloc = CCS*[X Y]'; % Change to crack tip coordinates
r = sqrt(XYloc(1)^2+XYloc(2)^2); % Radius from crack tip to current gauss point
if r < 0.001*lXElem; r = 0.05*lXElem; end
theta = atan2(XYloc(2),XYloc(1)); % Angle from crack tip to current gauss point
elseif nCT == 2
X1 = Xgp-xCT(1);
Y1 = Ygp-yCT(1);
X2 = Xgp-xCT(2);
Y2 = Ygp-yCT(2);
CCS = [cos(omega(1)) sin(omega(1));-sin(omega(1)) cos(omega(1))];
XY1 = CCS*[X1 Y1]';
CCS = [cos(omega(2)) sin(omega(2));-sin(omega(2)) cos(omega(2))];
XY2 = CCS*[X2 Y2]';
r1 = sqrt(XY1(1)^2+XY1(2)^2); % Radius from crack tip to current gauss point
r2 = sqrt(XY2(1)^2+XY2(2)^2);
if r1 > r2
r = r2; theta = atan2(XY2(2),XY2(1));
CCS = [cos(omega(2)) sin(omega(2));-sin(omega(2)) cos(omega(2))];
elseif r2 > r1
r = r1; theta = atan2(XY1(2),XY1(1));
CCS = [cos(omega(1)) sin(omega(1));-sin(omega(1)) cos(omega(1))];
end
if r < 0.001*lXElem; r = 0.05*lXElem; end
end
c = 1/2/sqrt(r); ct = cos(theta+omega); st = sin(theta+omega); c11=E1^2/(E1-nuxy12^2*E2);c12=E1*E2*nuxy12/(E1-nuxy12^2*E2);c22= E1*E2/ (E1-nuxy12^2*E2);
c66=G12;A=1/2*(c66/c11+c22/c66-(c12+c66)^2/(c11*c66));P1=sqrt( A-sqrt(A^2-c22/c11)) ; P2=sqrt( A+sqrt(A^2-c22/c11)) ; d1=P1*sec(theta)^2/(P1^2+tan(theta)^2);
d2=P2*sec(theta)^2/(P2^2+tan(theta)^2);e1=(cos(theta)^2+sin(theta)^2/P1^2)^(-3/4);e2=(cos(theta)^2+sin(theta)^2/P2^2)^(-3/4);f1=(1-P1^2)/(2*P1^2);f2=(1-P2^2)/(2*P2^2);
g1= sqrt(cos(theta)^2+sin(theta)^2/P1^2); g2= sqrt(cos(theta)^2+sin(theta)^2/P2^2); theta1=atan(tan(theta)/P1); theta2=atan(tan(theta)/P2); % Constants
if NN(iN,12) == 0 % Crack tip enrichment
a1gp = sqrt(r)*sin(theta1/2)*sqrt(g1); % Node crack tip enrichment value1
a2gp = sqrt(r)*sin(theta2/2)*sqrt(g2); % Node crack tip enrichment value2
a1 = a1gp-NN(iN,5); % Shifted alpha 1 enrichment value
a2 = a2gp-NN(iN,7); % Shifted alpha 2 enrichment value
% Derivative of crack tip enrichment functions with respect to X
Px = c*[ct*sin(theta1/2)*sqrt(g1)-st*(d1*sqrt(g1)*cos(theta1/2)+f1*e1*sin(theta1/2)*sin(2*theta));...
ct*sin(theta2/2)*sqrt(g2)-st*(d2*sqrt(g2)*cos(theta2/2)+f2*e2*sin(theta2/2)*sin(2*theta))];
% Derivative of crack tip enrichment functions with respect to Y
Py = c*[st*sin(theta1/2)*sqrt(g1)+ct*(d1*sqrt(g1)*cos(theta1/2)+f1*e1*sin(theta1/2)*sin(2*theta));...
st*sin(theta2/2)*sqrt(g2)+ct*(d2*sqrt(g2)*cos(theta2/2)+f2*e2*sin(theta2/2)*sin(2*theta))];
B1 = [Nx(iN)*a1+N(iN)*Px(1);
Ny(iN)*a1+N(iN)*Py(1)];
B2 = [Nx(iN)*a2+N(iN)*Px(2) ;
Ny(iN)*a2+N(iN)*Py(2)];
Bb = [B1 B2 ];
Benr(:,indexindex+1)) = Bb;
index = index+2;
if (i == length(gp))
local(iLociLoc+1)) = [NN(iN,4) NN(iN,6)]; %Local is the vector which maps the local stiffness matrix into the global stiffness matrix
iLoc = iLoc+2; % enrElem=[ 46 47 51 52 53 56 57 58 59 61 62 63 64 67 68 69];
end
end
end
end
B = [Bt Benr];
localKT = localKT + W*B'*C*B*detJ;
end
end
if length(localKT) == 4 % Unenriched element
for row = 1:4
for col = 1:4
nIndexT = nIndexT+1;
allRowT(nIndexT) = local(row);
allColT(nIndexT) = local(col);
allValT(nIndexT) = localKT(row,col);
end
end
else
globalKT(local,local) = globalKT(local,local) + localKT;
% Assemble the global stiffness globalKT =spdiags((globalKT)+(localKT)); globalKT=zeros (4000,8000);
end
end
globalKT = globalKT + sparse(allRowT,allColT,allValT,globalTDOF,globalTDOF);
，
z770428
??? Error using ==> plus
Matrix dimensions must agree.
Error in ==> conductMatrix at 313
globalKT(local,local) = globalKT(local,local) + localKT;
1528
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
