51 lines
1.7 KiB
Matlab
51 lines
1.7 KiB
Matlab
%%%%%%%%%%% 一阶六面体单元求解程序(求解位移矩阵) %%%%%%%%%%%
|
|
% U位移矩阵
|
|
% E弹性模量
|
|
% u泊松比
|
|
% Forces外力矩阵n*3 [节点 方向 大小]
|
|
% Constraints位移约束n*3 [节点 方向 大小]
|
|
% Nodes节点坐标信息
|
|
% Elements单元信息
|
|
function [U]=StaticsSolver(E,u,epsilon0,Constraints,Nodes,Elements)
|
|
Dof=3;
|
|
NodeCount = size(Nodes,1); % 节点个数
|
|
ElementCount= size(Elements,1); %单元个数
|
|
Dofs = Dof*NodeCount; %总自由度数
|
|
U=sparse(Dofs,1); % 初始化结构位移
|
|
K = sparse(Dofs,Dofs); %初始化总体刚度阵
|
|
Force = sparse(Dofs,1); %初始化外力向量
|
|
%计算应力-应变矩阵D
|
|
D=LinearIsotropicD(E,u);
|
|
for I=1:ElementCount
|
|
% 单元节点坐标
|
|
ElementNodeCoordinate=Nodes(Elements(I,:),:);
|
|
% 计算单刚
|
|
[ElementStiffnessMatrix,P_dT]=Ke(D,ElementNodeCoordinate,epsilon0);%返回单元等效温度载荷
|
|
% 计算单元节点自由度编号
|
|
ElementNodeDOF=zeros(1,24);
|
|
for J=1:8
|
|
II=(J-1)*Dof+1;
|
|
ElementNodeDOF(II:II+2)=(Elements(I,J)-1)*Dof+1:(Elements(I,J)-1)*Dof+3;
|
|
end
|
|
K(ElementNodeDOF,ElementNodeDOF)=K(ElementNodeDOF,ElementNodeDOF)+ElementStiffnessMatrix;%根据单元节点自由度 组装总刚
|
|
Force(ElementNodeDOF)=Force(ElementNodeDOF)+P_dT;%%等效温度载荷
|
|
end
|
|
% 施加外力
|
|
% if size(Forces,1)>0
|
|
% ForceDOF = Dof*(Forces(:,1)-1)+Forces(:,2); %计算外力自由度编号
|
|
% Force(ForceDOF) = Force(ForceDOF) + Forces(:,3);
|
|
% end
|
|
% 乘大数法施加位移约束
|
|
BigNumber=1e8;
|
|
ConstraintsNumber=size(Constraints,1);
|
|
if ConstraintsNumber~=0
|
|
FixedDof=Dof*(Constraints(:,1)-1)+Constraints(:,2); %被约束的自由度编号(列向量)
|
|
for i=1:ConstraintsNumber
|
|
K(FixedDof(i),FixedDof(i))=K(FixedDof(i),FixedDof(i))*BigNumber;
|
|
Force(FixedDof(i))=Constraints(i,3)*K(FixedDof(i),FixedDof(i));
|
|
end
|
|
end
|
|
%计算位移
|
|
U = K\Force;
|
|
save data Force
|
|
end |