 |
 |
 |
| 24小時(shí)熱門版塊排行榜 |
- >論壇更新日志 (3670)
- >導(dǎo)師招生 (416)
- >考研 (415)
- >文獻(xiàn)求助 (177)
- >蟲(chóng)友互識(shí) (155)
- >基金申請(qǐng) (129)
- >休閑灌水 (85)
- >論文投稿 (63)
- >考博 (55)
- >教師之家 (54)
- >招聘信息布告欄 (51)
- >碩博家園 (41)
- >博后之家 (31)
- >論文道賀祈福 (29)
- >留學(xué)生活 (27)
- >無(wú)機(jī)/物化 (22)
 |
 |
 |
||
 |
 |
 |
 |
 |
| 查看: 2798 | 回復(fù): 12 | |||
z770428
|
[求助]
MATLAB FEM 剛度矩陣
|
| 各位高手，我用MATLAB編程得到的熱傳導(dǎo)矩陣中為什么有些對(duì)角線元是0？應(yīng)該全大于0呀，是可能哪兒錯(cuò)了 |
回復(fù)此樓» 猜你喜歡
2026年循環(huán)經(jīng)濟(jì)功能材料國(guó)際會(huì)議（ICFMCE 2026）
已經(jīng)有0人回復(fù)
-
2026年第五屆電氣、電子與信息工程國(guó)際會(huì)議（ISEEIE 2026）
已經(jīng)有0人回復(fù)
-
物理學(xué)I論文潤(rùn)色/翻譯怎么收費(fèi)?
已經(jīng)有84人回復(fù)
-
哈爾濱理工大學(xué)物理系招收物理學(xué)考研調(diào)劑
已經(jīng)有5人回復(fù)
-
0702一志愿吉大B區(qū)求調(diào)劑
已經(jīng)有5人回復(fù)
-
求調(diào)劑
已經(jīng)有0人回復(fù)
-
0702一志愿吉大B區(qū)求調(diào)劑有論文
已經(jīng)有0人回復(fù)
-
請(qǐng)問(wèn)還有沒(méi)有用Latex寫(xiě)文章的小伙伴們？
已經(jīng)有0人回復(fù)
-
光學(xué)工程學(xué)碩調(diào)劑信息
已經(jīng)有27人回復(fù)
-
歡迎加入課題組
已經(jīng)有0人回復(fù)
-
散金幣，求好運(yùn)，祝面上順利！
已經(jīng)有44人回復(fù)
» 本主題相關(guān)價(jià)值貼推薦，對(duì)您同樣有幫助:
- matlab 軸的強(qiáng)度與剛度校核軟件
已經(jīng)有3人回復(fù)
csgt0
榮譽(yù)版主 (著名寫(xiě)手)
彩色掛圖
-
專家經(jīng)驗(yàn): +2 - 計(jì)算強(qiáng)帖: 2
- 應(yīng)助: 367 (碩士)
- 貴賓: 0.457
- 金幣: 7807
- 散金: 8897
- 紅花: 20
- 沙發(fā): 7
- 帖子: 2809
- 在線: 515.8小時(shí)
- 蟲(chóng)號(hào): 449046
- 注冊(cè): 2007-11-02
- 專業(yè): 金屬非晶態(tài)、準(zhǔn)晶和納米晶
- 管轄: 計(jì)算模擬
2樓2012-09-14 15:54:01
z770428 
金蟲(chóng) (正式寫(xiě)手)
- 應(yīng)助: 0 (幼兒園)
- 金幣: 1508.1
- 紅花: 2
- 帖子: 416
- 在線: 178.3小時(shí)
- 蟲(chóng)號(hào): 1360748
- 注冊(cè): 2011-08-04
- 性別: GG
- 專業(yè): 固體力學(xué)
3樓2012-09-17 07:36:32
zhuomec
新蟲(chóng) (正式寫(xiě)手)
- 應(yīng)助: 21 (小學(xué)生)
- 金幣: 3927.6
- 散金: 39
- 紅花: 8
- 帖子: 805
- 在線: 249.9小時(shí)
- 蟲(chóng)號(hào): 1723083
- 注冊(cè): 2012-03-28
- 專業(yè): 固體力學(xué)
4樓2012-09-17 08:31:57
z770428
金蟲(chóng) (正式寫(xiě)手)
- 應(yīng)助: 0 (幼兒園)
- 金幣: 1508.1
- 紅花: 2
- 帖子: 416
- 在線: 178.3小時(shí)
- 蟲(chóng)號(hào): 1360748
- 注冊(cè): 2011-08-04
- 性別: GG
- 專業(yè): 固體力學(xué)
5樓2012-09-17 12:35:37
csgt0
榮譽(yù)版主 (著名寫(xiě)手)
彩色掛圖
-
專家經(jīng)驗(yàn): +2 - 計(jì)算強(qiáng)帖: 2
- 應(yīng)助: 367 (碩士)
- 貴賓: 0.457
- 金幣: 7807
- 散金: 8897
- 紅花: 20
- 沙發(fā): 7
- 帖子: 2809
- 在線: 515.8小時(shí)
- 蟲(chóng)號(hào): 449046
- 注冊(cè): 2007-11-02
- 專業(yè): 金屬非晶態(tài)、準(zhǔn)晶和納米晶
- 管轄: 計(jì)算模擬
6樓2012-09-17 13:21:44
z770428
金蟲(chóng) (正式寫(xiě)手)
- 應(yīng)助: 0 (幼兒園)
- 金幣: 1508.1
- 紅花: 2
- 帖子: 416
- 在線: 178.3小時(shí)
- 蟲(chóng)號(hào): 1360748
- 注冊(cè): 2011-08-04
- 性別: GG
- 專業(yè): 固體力學(xué)
|
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) ;... % 溫度場(chǎng)一樣，名稱改為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 ;... % 溫度場(chǎng)分析同 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
金蟲(chóng) (正式寫(xiě)手)
- 應(yīng)助: 0 (幼兒園)
- 金幣: 1508.1
- 紅花: 2
- 帖子: 416
- 在線: 178.3小時(shí)
- 蟲(chóng)號(hào): 1360748
- 注冊(cè): 2011-08-04
- 性別: GG
- 專業(yè): 固體力學(xué)
8樓2012-09-17 15:49:16
z770428
金蟲(chóng) (正式寫(xiě)手)
- 應(yīng)助: 0 (幼兒園)
- 金幣: 1508.1
- 紅花: 2
- 帖子: 416
- 在線: 178.3小時(shí)
- 蟲(chóng)號(hào): 1360748
- 注冊(cè): 2011-08-04
- 性別: GG
- 專業(yè): 固體力學(xué)
9樓2012-09-17 15:51:19
csgt0
榮譽(yù)版主 (著名寫(xiě)手)
彩色掛圖
-
專家經(jīng)驗(yàn): +2 - 計(jì)算強(qiáng)帖: 2
- 應(yīng)助: 367 (碩士)
- 貴賓: 0.457
- 金幣: 7807
- 散金: 8897
- 紅花: 20
- 沙發(fā): 7
- 帖子: 2809
- 在線: 515.8小時(shí)
- 蟲(chóng)號(hào): 449046
- 注冊(cè): 2007-11-02
- 專業(yè): 金屬非晶態(tài)、準(zhǔn)晶和納米晶
- 管轄: 計(jì)算模擬
10樓2012-09-17 17:09:35
相關(guān)版塊跳轉(zhuǎn)
我要訂閱樓主 z770428 的主題更新
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 085601材料工程專碩求調(diào)劑 +5 | 慕寒mio 2026-03-16 | 5/250 |
|
|---|---|---|---|---|
|
|
[考研] 268求調(diào)劑 +7 | 好運(yùn)連綿不絕 2026-03-12 | 8/400 |
|
|
|
[考研] 085601求調(diào)劑 +4 | Du.11 2026-03-16 | 4/200 |
|
|
|
[考博] 26申博 +4 | 八旬速覽 2026-03-16 | 4/200 |
|
|
[考研] 275求調(diào)劑 +4 | 太陽(yáng)花天天開(kāi)心 2026-03-16 | 4/200 |
|
|
|
[考研] 289求調(diào)劑 +6 | 步川酷紫123 2026-03-11 | 6/300 |
|
|
|
[考研] 285化工學(xué)碩求調(diào)劑（081700） +9 | 柴郡貓_ 2026-03-12 | 9/450 |
|
|
|
[考研] 267一志愿南京工業(yè)大學(xué)0817化工求調(diào)劑 +6 | SUICHILD 2026-03-12 | 6/300 |
|
|
|
[考研] 一志愿985，本科211，0817化學(xué)工程與技術(shù)319求調(diào)劑 +5 | Liwangman 2026-03-15 | 5/250 |
|
|
|
[考研] 304求調(diào)劑 +5 | 素年祭語(yǔ) 2026-03-15 | 5/250 |
|
|
|
[考研] 304求調(diào)劑 +4 | ahbd 2026-03-14 | 4/200 |
|
|
|
[考研] 一志愿211 0703方向310分求調(diào)劑 +3 | 努力奮斗112 2026-03-15 | 3/150 |
|
|
|
[考博] 東華理工大學(xué)化材專業(yè)26屆碩士博士申請(qǐng) +6 | zlingli 2026-03-13 | 6/300 |
|
|
|
[考研] 22408總分284求調(diào)劑 +3 | InAspic 2026-03-13 | 3/150 |
|
|
|
[考研] 266求調(diào)劑 +4 | 學(xué)員97LZgn 2026-03-13 | 4/200 |
|
|
|
[考研] 290求調(diào)劑 +9 | ADT 2026-03-11 | 9/450 |
|
|
|
[考研] 281求調(diào)劑 +9 | Koxui 2026-03-12 | 11/550 |
|
|
|
[考研] 311求調(diào)劑 +3 | 冬十三 2026-03-13 | 3/150 |
|
|
|
[考研] 求b區(qū)學(xué)校調(diào)劑 +3 | 周56 2026-03-11 | 3/150 |
|
|
|
[考研] 工科材料085601 279求調(diào)劑 +8 | 困于星晨 2026-03-12 | 10/500 |
|