| 查看: 2685 | 回复: 12 | ||
z770428
|
[求助]
MATLAB FEM 刚度矩阵
|
| 各位高手，我用MATLAB编程得到的热传导矩阵中为什么有些对角线元是0？应该全大于0呀，是可能哪儿错了 |
回复此楼» 猜你喜欢
求标准粉末衍射卡号 ICDD 01-076-1802
已经有0人回复
-
新西兰Robinson研究所招收全奖PhD
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有286人回复
-
石墨烯转移--二氧化硅衬底石墨烯
已经有0人回复
-
笼目材料中量子自旋液体基态的证据
已经有0人回复
-
数学教学论硕士可以读数学物理博士吗？
已经有0人回复
-
德国亥姆霍兹Hereon中心汉堡分部招镁合金腐蚀裂变SCC课题方向2026公派博士生
已经有4人回复
-
澳门大学 应用物理及材料工程研究院 潘晖教授课题组诚招博士后
已经有11人回复
-
求助NH4V4O10晶体的CIF文件
已经有0人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- matlab 轴的强度与刚度校核软件
已经有3人回复
2楼2012-09-14 15:54:01
z770428 
金虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 1508.1
- 红花: 2
- 帖子: 416
- 在线: 178.3小时
- 虫号: 1360748
- 注册: 2011-08-04
- 性别: GG
- 专业: 固体力学
3楼2012-09-17 07:36:32
zhuomec
新虫 (正式写手)
- 应助: 21 (小学生)
- 金币: 3927.6
- 散金: 39
- 红花: 8
- 帖子: 805
- 在线: 249.9小时
- 虫号: 1723083
- 注册: 2012-03-28
- 专业: 固体力学
4楼2012-09-17 08:31:57
z770428
金虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 1508.1
- 红花: 2
- 帖子: 416
- 在线: 178.3小时
- 虫号: 1360748
- 注册: 2011-08-04
- 性别: GG
- 专业: 固体力学
5楼2012-09-17 12:35:37
6楼2012-09-17 13:21:44
z770428
金虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 1508.1
- 红花: 2
- 帖子: 416
- 在线: 178.3小时
- 虫号: 1360748
- 注册: 2011-08-04
- 性别: GG
- 专业: 固体力学
|
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 = 1  nYElem*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 elementCTN = 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(:,index index+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(:,index index+1)) = Bb;index = index+2; if (i == length(gp)) local(iLoc iLoc+1)) = [NN(iN,4) NN(iN,6)]; %Local is the vector which maps the local stiffness matrix into the global stiffness matrixiLoc = 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); |
7楼2012-09-17 15:45:41
z770428
金虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 1508.1
- 红花: 2
- 帖子: 416
- 在线: 178.3小时
- 虫号: 1360748
- 注册: 2011-08-04
- 性别: GG
- 专业: 固体力学
8楼2012-09-17 15:49:16
z770428
金虫 (正式写手)
- 应助: 0 (幼儿园)
- 金币: 1508.1
- 红花: 2
- 帖子: 416
- 在线: 178.3小时
- 虫号: 1360748
- 注册: 2011-08-04
- 性别: GG
- 专业: 固体力学
9楼2012-09-17 15:51:19
10楼2012-09-17 17:09:35