Q=Qbarra;
% Astiff(5,5)=0;Bstiff(5,5)=0;Fstiff(5,5)=0;Istiff(5,5)=0;
for k=1:size(alfas,2) for i=1:3
for j=1:3 Astiff(i,j)=Astiff(i,j)+Q(i,j,k)*(z(k+1)-z(k));
Bstiff(i,j)=Bstiff(i,j)+Q(i,j,k)*(z(k+1)^2-z(k)^2)/2;
Fstiff(i,j)=Fstiff(i,j)+Q(i,j,k)*(z(k+1)^3-z(k)^3)/3; end
end
for i=4:5 for j=4:5
Istiff(i,j)=Istiff(i,j)+Q(i,j,k)*(z(k+1)-z(k)); end
Có thể bạn quan tâm!
-
Kết Quả Phát Hiện Hư Hỏng Của Pso, Ann Và Annpso -
L. F. F. Miguel, R. H. Lopez, And L. F. F. Miguel. A Hybrid Approach For Damage Detection Of Structures Under Operational Conditions. Journal Of Sound And Vibration. 2013, 332, 4241-4260. -
N. M. Nawi, R. Ransing, And M. Ransing. An Improved Conjugate Gradient Based Learning Algorithm For Back Propagation Neural Networks. International Journal Of Computational Intelligence. 2007, -
Chẩn đoán dầm cầu bằng phương pháp phân tích dao động trên mô hình số hoá kết cấu được cập nhật sử dụng thuật toán tối ưu hoá bầy đàn kết hợp mạng nơ ron nhân tạo - 18 -
Chẩn đoán dầm cầu bằng phương pháp phân tích dao động trên mô hình số hoá kết cấu được cập nhật sử dụng thuật toán tối ưu hoá bầy đàn kết hợp mạng nơ ron nhân tạo - 19
Xem toàn bộ 154 trang tài liệu này.
end
end
pi=double(pi); % come back to numeric computation
% constitutive matrices CMembranaMembrana=Astiff(1:3,1:3); CMembranaFlexao0=Bstiff(1:3,1:3); CFlexao0Flexao0=Fstiff(1:3,1:3); CCorte0Corte0=Istiff(4:5,4:5);
%Mesh generation L = 1;
numberElementsX=10; numberElementsY=10;
numberElements=numberElementsX*numberElementsY; [nodeCoordinates, elementNodes] = ...
MeshRectangular(L,L,numberElementsX,numberElementsY); xx=nodeCoordinates(:,1);
yy=nodeCoordinates(:,2);
% figure
% PlotMesh(nodeCoordinates,elementNodes);
% axis off
numberNodes=size(xx,1); % the number of Nodes
% drawingEigenmodes3d_ele(nodeCoordinates,elementNodes)
% figure
% drawingMesh(nodeCoordinates,elementNodes);
% GDof: global number of degrees of freedom GDof=5*numberNodes;
%%
AMatrix=CMembranaMembrana; BMatrix=CMembranaFlexao0; DMatrix=CFlexao0Flexao0; SMatrix=CCorte0Corte0;
% Gauss quadrature for bending part [~,activeDof,fixedNodeW]=...
EssentialBC5dof('cccc',GDof,xx,yy,nodeCoordinates,numberNodes);
%% Mass matrix M=zeros(GDof); nb_mode=length(activeDof);
mt_e=zeros(nb_mode,nb_mode,numberElements);
% Gauss quadrature for bending part [gaussWeights,gaussLocations]=gaussQuadrature('complete');
% cycle for element
for e=1:numberElements
% indice : nodal condofectivities for each element indice=elementNodes(e,:);
% cycle for Gauss point for q=1:size(gaussWeights,1)
GaussPoint=gaussLocations(q,:); xi=GaussPoint(1); eta=GaussPoint(2);
% shape functions and derivatives [shapeFunction,naturalDerivatives]=shapeFunctionQ4(xi,eta);
% Jacobian matrix, inverse of Jacobian, derivatives w.r.t. x,y [Jacob,~,~]=...
Jacobian(nodeCoordinates(indice,:),naturalDerivatives);
% mass matrix M(indice,indice)=M(indice,indice)+...
shapeFunction*shapeFunction'*thickness*rho*gaussWeights(q)*det(Jacob); M(indice+numberNodes,indice+numberNodes)=...
M(indice+numberNodes,indice+numberNodes)+... shapeFunction*shapeFunction'*I*rho*gaussWeights(q)*det(Jacob);
M(indice+2*numberNodes,indice+2*numberNodes)=...
M(indice+2*numberNodes,indice+2*numberNodes)+... shapeFunction*shapeFunction'*I*rho*gaussWeights(q)*det(Jacob);
M(indice+3*numberNodes,indice+3*numberNodes)=...
M(indice+3*numberNodes,indice+3*numberNodes)+... shapeFunction*shapeFunction'*thickness*rho*gaussWeights(q)*det(Jacob);
M(indice+4*numberNodes,indice+4*numberNodes)=...
M(indice+4*numberNodes,indice+4*numberNodes)+... shapeFunction*shapeFunction'*thickness*rho*gaussWeights(q)*det(Jacob);
end % Gauss point
mt_e(:,:,e) =M(activeDof,activeDof); end % element
Mg=M;
M=M(activeDof,activeDof);
%% Stiffness matrix
% Les Matrices de Rigidité Saines K=zeros(GDof); kt_e=zeros(nb_mode,nb_mode,numberElements); Kd=zeros(GDof); kt_d=zeros(nb_mode,nb_mode,numberElements);
%
Element_endo=zeros(1,numberElements); Pourcentage=zeros(1,numberElements);
% Introduction an damage
alpha=zeros(1,numberElements ); alpha(10)=.9;
dim=4; nument=5;
K_eu=zeros(dim*nument,dim*nument,numberElements); K_ed=zeros(dim*nument,dim*nument,numberElements);
% Gauss quadrature for bending part [gaussWeights,gaussLocations]=gaussQuadrature('complete');
% cycle for element
for e=1:numberElements Element_endo(e)=e; Pourcentage(e)=alpha(1,e);
indice : nodal condofectivities for each element elementDof: element degrees of freedom indice=elementNodes(e,:);
elementDof=[ indice indice+numberNodes indice+2*numberNodes ... indice+3*numberNodes indice+4*numberNodes];
ndof=length(indice);
% cycle for Gauss point
for q=1:size(gaussWeights,1) GaussPoint=gaussLocations(q,:); xi=GaussPoint(1); eta=GaussPoint(2);
% shape functions and derivatives [~,naturalDerivatives]=shapeFunctionQ4(xi,eta);
% Jacobian matrix, inverse of Jacobian,derivatives w.r.t. x,y [Jacob,~,XYderivatives]=Jacobian(nodeCoordinates(indice,:),naturalDerivatives);
% [B] matrix bending B_b=zeros(3,5*ndof);
B_b(1,ndof+1:2*ndof) = XYderivatives(:,1)'; B_b(2,2*ndof+1:3*ndof) = XYderivatives(:,2)'; B_b(3,ndof+1:2*ndof) = XYderivatives(:,2)'; B_b(3,2*ndof+1:3*ndof) = XYderivatives(:,1)';
% [B] matrix membrane B_m=zeros(3,5*ndof);
B_m(1,3*ndof+1:4*ndof) = XYderivatives(:,1)'; B_m(2,4*ndof+1:5*ndof) = XYderivatives(:,2)'; B_m(3,3*ndof+1:4*ndof) = XYderivatives(:,2)'; B_m(3,4*ndof+1:5*ndof) = XYderivatives(:,1)';
% stiffness matrix Healthy
% ... bending-bending K(elementDof,elementDof)=K(elementDof,elementDof)+...
B_b'*DMatrix*B_b*gaussWeights(q)*det(Jacob);
% ... membrane-membrane K(elementDof,elementDof)=K(elementDof,elementDof)+...
B_m'*AMatrix*B_m*gaussWeights(q)*det(Jacob);
% ... membrane-bending K(elementDof,elementDof)=K(elementDof,elementDof)+...
B_m'*BMatrix*B_b*gaussWeights(q)*det(Jacob);
% ... bending-membrane K(elementDof,elementDof)=K(elementDof,elementDof)+...
B_b'*BMatrix*B_m*gaussWeights(q)*det(Jacob);
% stiffness matrix Damaged
% ... bending-bending Kd(elementDof,elementDof)=Kd(elementDof,elementDof)+...
B_b'*DMatrix*B_b*gaussWeights(q)*det(Jacob);
% ... membrane-membrane Kd(elementDof,elementDof)=Kd(elementDof,elementDof)+...
B_m'*AMatrix*B_m*gaussWeights(q)*det(Jacob);
% ... membrane-bending Kd(elementDof,elementDof)=Kd(elementDof,elementDof)+...
B_m'*BMatrix*B_b*gaussWeights(q)*det(Jacob);
% ... bending-membrane Kd(elementDof,elementDof)=Kd(elementDof,elementDof)+...
B_b'*BMatrix*B_m*gaussWeights(q)*det(Jacob); end % Gauss point
end % element
% Gauss quadrature for shear part [gaussWeights,gaussLocations]=gaussQuadrature('reduced');
% cycle for element
for e=1:numberElements Element_endo(e)=e; Pourcentage(e)=alpha(1,e);
% % Element_endo(e)=e;
% % Pourcentage(e)=alpha(1,e);
% indice : nodal condofectivities for each element
% elementDof: element degrees of freedom indice=elementNodes(e,:);
elementDof=[ indice indice+numberNodes indice+2*numberNodes ... indice+3*numberNodes indice+4*numberNodes];
ndof=length(indice); nndof(:,e)=elementDof;
% cycle for Gauss point
for q=1:size(gaussWeights,1) GaussPoint=gaussLocations(q,:); xi=GaussPoint(1); eta=GaussPoint(2);
% shape functions and derivatives [shapeFunction,naturalDerivatives]=shapeFunctionQ4(xi,eta);
% Jacobian matrix, inverse of Jacobian, derivatives w.r.t. x,y [Jacob,~,XYderivatives]=Jacobian(nodeCoordinates(indice,:),naturalDerivatives);
% [B] matrix shear B_s=zeros(2,5*ndof);
B_s(1,1:ndof) = XYderivatives(:,1)'; B_s(2,1:ndof) = XYderivatives(:,2)'; B_s(1,ndof+1:2*ndof) = shapeFunction; B_s(2,2*ndof+1:3*ndof)= shapeFunction;
% Stifness Matrix Healthy K(elementDof,elementDof)=K(elementDof,elementDof)+...
B_s'*SMatrix*B_s*gaussWeights(q)*det(Jacob); localstiff_U = zeros(dim*nument, dim*nument);
localstiff_U=localstiff_U+B_s'*SMatrix*B_s*gaussWeights(q)*det(Jacob); K_eu(:,:,e)= localstiff_U;
% Stifness Matrix Damaged SS=1-Pourcentage(e);
Kd(elementDof,elementDof)=Kd(elementDof,elementDof)+...
SS*B_s'*SMatrix*B_s*gaussWeights(q)*det(Jacob); localstiff_D=localstiff_U+SS*B_s'*SMatrix*B_s*gaussWeights(q)*det(Jacob); K_ed(:,:,e)= localstiff_D;
end % gauss point
% kt_e(:,:,e)=K(activeDof,activeDof);
% KT(:,:,e)=K;
% kt_d(:,:,e)=Kd(activeDof,activeDof); end % element
% Kg=K;
% Kdg=Kd; K=K(activeDof,activeDof); Kd=Kd(activeDof,activeDof);
% Les matrices élémentaires de dimension de activdof k_t=zeros(nb_mode,nb_mode,numberElements); m_t=zeros(nb_mode,nb_mode,numberElements); k_d=zeros(nb_mode,nb_mode,numberElements);
for tj=1:numberElements if tj==1
k_t(:,:,1)=kt_e(:,:,1);
m_t(:,:,1)=mt_e(:,:,1);
k_d(:,:,1)=kt_d(:,:,1);
K_T(:,:,1)=KT(:,:,1);
elseif tj==2
k_t(:,:,2)=kt_e(:,:,2)-kt_e(:,:,1);
m_t(:,:,2)=mt_e(:,:,2)-mt_e(:,:,1);
k_d(:,:,2)=kt_d(:,:,2)-kt_d(:,:,1);
K_T(:,:,2)=KT(:,:,2)-KT(:,:,1);
else
k_t(:,:,tj)=kt_e(:,:,tj)-kt_e(:,:,tj-1);
m_t(:,:,tj)=mt_e(:,:,tj)-mt_e(:,:,tj-1);
k_d(:,:,tj)=kt_d(:,:,tj)-kt_d(:,:,tj-1);
K_T(:,:,tj)=KT(:,:,tj)-KT(:,:,tj-1);
end end
[V,D] = eig(K,M); % V is to calculate mode shape, D is to calculate natural frequency.
% eigenproblem: free vibrations numberOfModes=15;
% Liew, p-Ritz
D0=e2*h^3/12/(1-miu12*miu21);
D = diag(sqrt(D)*L*L/pi/pi*sqrt(rho*h/D0)); [D,iid] = sort(D); ii = iid(1:numberOfModes); V_ef = V(:,iid);
VV = V(:,ii);
activeDofW=setdiff((1:numberNodes)',fixedNodeW); NNN=size(activeDofW); %NNN is the number of DOF
VVV(1:numberNodes,1:numberOfModes)=0; for i=1:numberOfModes
VVV(activeDofW,i)=VV(1:NNN,i); end
%%
drawingEigenmodes3d(numberNodes,numberOfModes,NNN,... numberElementsX,numberElementsY,...
zz=zeros(1,numberElements); for i=1:numberElements
zz(i)=norm(K_eu(:,:,i)-K_ed(:,:,i)); end
[V_D,D_D] = eig(Kd,M);
% eigenproblem: free vibrations numberOfModes=15;