Trang chủ / Phân Tích Dao Động

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 - 17


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!

Xem toàn bộ 154 trang tài liệu này.

end

end

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 - 17

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;

..... Xem trang tiếp theo?
⇦ Trang trước - Trang tiếp theo ⇨

Ngày đăng: 14/07/2022