close all

clear all

clc

%% 离散曲率值矩阵k，弧长矩阵ss，(包含边界点k=0,s=0)求出对应弧段的弧长-曲率线性方程

%Strain=[0 1.75044E-04 1.53376E-04 1.32575E-04 1.12096E-04 9.16930E-05 7.12955E-05 5.08464E-05 3.01552E-05]%应变

%Dspl=[0 -8.57915E-04 -2.36228E-03 -4.51358E-03 -7.22375E-03 -1.04064E-02 -1.39754E-02 -1.78449E-02 -2.19286E-02]%位移

Data= load('F:\ANSYS working directory\modal2\simulation.txt'); %从数据采集文件中提取

Datas=Data(:,1);

Strain=[zeros(1,4);reshape(Datas,8,4)];

s(1)=0;i=2:9;s(i)=0.15+0.1*(i-2); %弧长：单位m

ss=repmat(s',1,4);

h=0.005;%厚度

k=2*Strain/h;%曲率

for i=1:1:8

for j=1:4

a(i,j)=(k(i+1,j)-k(i,j))/(ss(i+1,j)-ss(i,j)); %线性方程的参数

b(i,j)=(k(i,j)*ss(i+1,j)-k(i+1,j)*ss(i,j))/(ss(i+1,j)-ss(i,j));

end

end

%% 每一个弧段上每个点的坐标差

syms p

step=0.01; %分段步长

u=1/step;

for j=1:4

X(1,j)=0;Z(1,j)=0; %先计算起始点

k1=(k(2,j)/ss(2,j))*step;

X(2,j)= sin(k1)/k1;

Z(2,j)=(1-cos(k1))/k1;

for i=1:1:8

for q=(ss(i,j)+step):step:ss(i+1,j) %对第i段分段

m=q; %求出微元弧段端点的弧长位置和曲率

n=q+step;%微元弧段末点

theta=int(a(i,j)*p+b(i,j),m,n); %坐标系旋转的角度

id = round(q*u+1); %索引取正整数

X(id+1,j)=X(id,j)*cos(theta)-Z(id,j)*sin(theta); %坐标增量

Z(id+1,j)=X(id,j)*sin(theta)+Z(id,j)*cos(theta);

end

end

end

%% 计算每个点的坐标

for j=1:4

XX(1,j)=X(1,j);

ZZ(1,j)=Z(1,j);

for i=1:1:ss(9,j)/step

XX(i+1,j)=XX(i,j)+X(i+1,j);

ZZ(i+1,j)=ZZ(i,j)+Z(i+1,j);

end

end

XXX = XX*step; %X坐标矩阵

ZZZ = -ZZ*step;%Z坐标矩阵

mm=size(XXX,1); %获得坐标矩阵的行数

j=1:4;YY(j)=0.02+0.06*(j-1);

YYY=repmat(YY,mm,1);%Y坐标矩阵

%% 对比校核

Datax=Data(:,2);

Datay=Data(:,3);

Dataz=Data(:,4);

xb=[zeros(1,4);reshape(Datax,8,4)];

yb=[zeros(1,4);reshape(Datay,8,4)];

zb=[zeros(1,4);reshape(Dataz,8,4)];

%% 曲面插值

[X1,Y1]=meshgrid(0:0.001:1,0:0.001:0.22);

Z1=griddata(XXX,YYY,ZZZ,X1,Y1);

[X2,Y2]=meshgrid(0:0.0001:1,0:0.0005:0.22);

Z2=griddata(xb,yb,zb,X2,Y2);

%% 画图

figure

mesh(X1,Y1,Z1)

hold on

mesh(X2,Y2,Z2)

title('曲面重建');

xlabel('X'); ylabel('Y');zlabel('Z');

set(gca,'XTick',[0,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,1]);

set(gca,'YTick',[0,0.02,0.08,0.14,0.2,0.22]);

set(gca,'ZTick',[-0.025:0.005:0.005]);

axis([0,1,-0.05,0.25,-0.025,0.005])