如何设计并行计算，能使每个点的offset_angle取值范围为n，计算全部情况下的TE：clc
clear
close all
discrete_num=11;
n=-(0:0.001:0.1);
for i=1:size(n,2)
offset_angle=linspace(n(i),n(i), discrete_num);
TE(i)=solveEnergy(offset_angle,discrete_num);
end function TE= solveEnergy(offset_angle,discrete_num)
scale = 1;
angle_all = -pi* 10/180; % 总弧度
NumberOfCurve = 3; % 总曲线数,不能等于1
rotate_angle = 0; % 旋转角度
[X] = plotcurve4(scale,angle_all,discrete_num,NumberOfCurve,rotate_angle,offset_angle);
for strip = 1:NumberOfCurve-1
for i = 1:discrete_num-1
n1 = (strip-1)*discrete_num + i;
n2 = (strip-1)*discrete_num + i+1;
n3 = strip*discrete_num + i+1;
n4 = strip*discrete_num + i;
Panel{(strip-1)*(discrete_num-1)+i,1}= [n4, n3, n2, n1];
end
end
Node=X;
%% Set up boundary conditions
% 定义固定约束
n=discrete_num;
for j=1:2:NumberOfCurve
for i=n*(j-1)+2:j*n
st(i-1-(n+1)*(j-1)/2)=i;
end
end
stz=[st]';
m = size(Node,1);
Suppz=[stz,zeros(length(stz),1),zeros(length(stz),1),ones(length(stz),1)];
% Supp = [ 1, 0, 1, 1;
% (NumberOfCurve-1)*n+1,0,1,1;
% n+1,1,1,0;
% Suppz;
% ];
for j=2:2:NumberOfCurve
for i=n*(j-1)+2:j*n
st1(i-1-n*(j-1))=i;
end
end
stz1=st1';
Suppz1=[stz1,zeros(length(stz1),1),zeros(length(stz1),1),ones(length(stz1),1)];
Supp = [
n+1,1,1,1;
Suppz1;
];
% for j=3:2:NumberOfCurve
% for i=(j-2)*n+1:n*(j-1)
% inp(i-n*(j-1)/2)=i;
% end
% end
inp=[1,2*n+1];
% inp=[n+1];
%ff=2e3*ones(length(inp),1);
ff=1*ones(length(inp),1);
Load = [ inp',zeros(length(inp),1),zeros(length(inp),1),ff;];
%% Adopt generalized N5B8 model
nL=(n-1)*NumberOfCurve+(NumberOfCurve-1)*n+(NumberOfCurve-1)*(n-1)+4*(n-1)*(NumberOfCurve-1);
% 折痕标记
isCrease = false(nL,1);
I1=NumberOfCurve-2;
ver(1,:)=n:(I1+1)*(n-1);%竖直折痕
for i=1:NumberOfCurve-1
I=(n-1)*NumberOfCurve;
hor(1,(n-2)*(i-1)+1:i*(n-2))=(I+2)+n*(i-1):(I+2)+n*(i-1)+n-3;%水平折痕
end
crease_ids = [hor(1:n-2)];
for i=2:n-1
and(i,:)=[ver(i),hor(n-2+i-1)];
end
and=and(2:end,:);
and=reshape(and',1,[]);
and=[ver(1),and];
crease_ids=[crease_ids,and];
isCrease(crease_ids) = true;
% 定义折痕旋转刚度
K_theta=zeros(nL,1);
for i=ver(:)
K_theta(i)=0; % 竖直折痕刚度
end
for i=hor(:)
K_theta(i)=0.01; % 横向折痕刚度
end
K_theta=(K_theta(crease_ids));
%%
% AnalyInputOpt = struct(...
% 'ModelType','N5B8',...
% 'MaterCalib','manual',... % 'manual'
% 'LoadType','Force',...
% 'InitialLoadFactor', 0.00001,...
% 'MaxIcr', 800,...
% 'BarCM', @(Ex)Ogden(Ex, 1e10),...
% 'Abar', 1e-5,...
% 'Kf',K_theta,...
% 'Kb',1e4 ...
% );
AnalyInputOpt = struct(...
'ModelType','N5B8',...
'MaterCalib','manual',... % 'manual'
'LoadType','Displacement',...
'InitialLoadFactor', 0.00001,...
'DispStep', 100,...
'BarCM', @(Ex)Ogden(Ex, 1e10),...
'Abar', 1e-5,...
'Kf',K_theta,...
'Kb',1e4 ...
);
%% Adopt generalized N4B5 model
% AnalyInputOpt = struct(...
% 'ModelType','N4B5',...
% 'MaterCalib','manual',... % 'manual'
% 'BarCM', @(Ex)Ogden(Ex, 5e3),...
% 'Abar', 2.4652,...
% 'Kb',0.9726*([38.4187;38.4187;41.7612]/0.127).^(1/3)./[38.4187;38.4187;41.7612],...
% 'Kf',1,...
% 'RotSprBend', @SuperLinearBend,...
% 'RotSprFold', @(he,h0,Kf,L0)EnhancedLinear(he,h0,Kf,L0,15,345),...
% 'LoadType','Force',...
% 'InitialLoadFactor', 0.00001,...
% 'MaxIcr', 100,...
% 'StopCriterion',@(Node,U,icrm)(abs(U(5*3))>12));
%% Perform analysis
% Assemble input data
[truss, angles, AnalyInputOpt,Fold] = PrepareData(Node,Panel,Supp,Load,AnalyInputOpt);
% Specify initial deformation state
truss.U0 = zeros(3*size(truss.Node,1),1);
% Perform path-following analysis using MGDCM
[Uhis,Fhis] = PathAnalysis(truss,angles,AnalyInputOpt);
% Clean output data
Uhis = real(Uhis);
Fhis = real(Fhis);
STAT = PostProcess(Uhis,truss,angles);
TE=STAT.PE(end);
end
扫一扫
1666
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
