摩阻扭矩-轨迹参数计算#
井眼曲率可以改
traj_para_linear#
% 轨迹参数计算函数 从钻头往上计算
function [alpha,k,kal,kph,tau,dk_ds,tz,nz,bz,dalpha_ds,dpha_ds]=traj_para_linear(trajectory_parameters)
dL = 1; % 计算间距
depth_m = trajectory_parameters(:,1); % 导入井深
incl = trajectory_parameters(:,2)*pi/180; % 导入井斜角
azim = trajectory_parameters(:,3)*pi/180; % 导入井斜方位角
depth_traj_para = [0:dL:floor(max(depth_m))]'; % 计算井深
alpha = interp1(depth_m,incl,depth_traj_para,'linear'); %井斜角的三次样条插值，转换为弧度制
pha = interp1(depth_m,azim,depth_traj_para,'linear'); %方位角的三次样条插值，转换为弧度制
%% 计算轨迹参数，从井底往井口计算，故向前差分为 [(i-1) - (i)]/dL
dalpha_ds2 = zeros(length(depth_traj_para)); % 线性插值，井斜二次导数为零，井段断点处一阶导数可能不连续，也取零
dpha_ds2 = zeros(length(depth_traj_para)); % 线性插值，方位二次导数为零，井段断点处一阶导数可能不连续，也取零
for i=1:length(depth_traj_para)
if i == 1
dalpha_ds(i,:) = 0; % 井斜变化率
dpha_ds(i,:) = 0; % 方位变化率
dk_ds(i,:) = 0;
tz(i,:) = 1;
nz(i,:) = 0;
bz(i,:) = 0;
continue;
end
dalpha_ds(i,:) = (alpha(i-1)-alpha(i))/dL; % 井斜变化率
dpha_ds(i,:) = (pha(i-1)-pha(i))/dL; % 方位变化率
k(i,:) = ((dalpha_ds(i))^2+(sin(alpha(i)))^2*(dpha_ds(i))^2)^0.5; % 井眼曲率
dk_ds(i,:) = (k(i-1) - k(i))/dL; % 井眼曲率变化率
kal(i,:) = dalpha_ds(i); % 井斜平面井眼曲率
kph(i,:) = sin(alpha(i))*dpha_ds(i); % 方位平面井眼曲率
tz(i,:) = cos(alpha(i)); % 切线方向投影分量
if dalpha_ds(i) == 0 && dpha_ds(i) ~= 0
nz(i,:) = 0;
bz(i,:) = sin(alpha(i));
elseif dpha_ds(i) == 0 || (dalpha_ds(i) == 0 && dpha_ds(i) == 0)
nz(i,:) = sin(alpha(i));
bz(i,:) = 0;
else
nz(i,:) = dalpha_ds(i)*sin(alpha(i))/k(i); % 主法线方向投影分量
bz(i,:) = -dpha_ds(i)*(sin(alpha(i)))^2/k(i); % 次法线方向投影分量
end
if k(i)<1e-6
k(i,:)=1e-6; % 如果井眼曲率k小于1e-6，为计算准确，令k=1e-6
end
tau(i,:)=(sin(alpha(i))*(dalpha_ds(i)*dpha_ds2(i)-dpha_ds(i)*dalpha_ds2(i))+dpha_ds(i)*((dalpha_ds(i))^2+(k(i))^2)*cos(alpha(i)))/(k(i))^2; % 井眼挠度
end
end
traj_para_pchip#
% 轨迹参数计算函数 从钻头往上计算
function [alpha,k,kal,kph,tau,dk_ds,tz,nz,bz,dalpha_ds,dpha_ds]=traj_para_pchip(trajectory_parameters)
dL = 1;
depth_m = trajectory_parameters(:,1); % 导入井深
incl = trajectory_parameters(:,2)*pi/180; % 导入井斜角
azim = trajectory_parameters(:,3)*pi/180; % 导入井斜方位角
depth_traj_para = [0:dL:floor(max(depth_m))]'; % 计算井深
alpha = interp1(depth_m,incl,depth_traj_para,'pchip'); %井斜角的三次样条插值，转换为弧度制
pha = interp1(depth_m,azim,depth_traj_para,'pchip'); %方位角的三次样条插值，转换为弧度制
%% 计算轨迹参数，从井底往井口计算，故向前差分为 [(i-1) - (i)]/dL
dalpha_ds = zeros(length(depth_traj_para),1);
dpha_ds = zeros(length(depth_traj_para),1);
dalpha_ds2 = zeros(length(depth_traj_para),1);
dpha_ds2 = zeros(length(depth_traj_para),1);
k = zeros(length(depth_traj_para),1);
dk_ds = zeros(length(depth_traj_para),1);
kal = zeros(length(depth_traj_para),1);
kph = zeros(length(depth_traj_para),1);
tz = zeros(length(depth_traj_para),1);
nz = zeros(length(depth_traj_para),1);
bz = zeros(length(depth_traj_para),1);
tau = zeros(length(depth_traj_para),1);
for i=1:length(depth_traj_para)
if i == 1
dalpha_ds(i,:) = 0; % 井斜变化率
dpha_ds(i,:) = 0; % 方位变化率
dalpha_ds2(i,:) = 0;
dpha_ds2(i,:) = 0;
k(i,:) = 0;
dk_ds(i,:) = 0;
tz(i,:) = 1;
nz(i,:) = 0;
bz(i,:) = 0;
tau(i,:) = 0;
continue;
end
dalpha_ds(i,:)=(alpha(i-1)-alpha(i))/dL; % 井斜变化率
dalpha_ds2(i,:) = (dalpha_ds(i-1)-dalpha_ds(i))/dL; % 井斜二阶倒数
dpha_ds(i,:) = (pha(i-1)-pha(i))/dL; % 方位变化率
dpha_ds2(i,:) = (dpha_ds(i-1)-dpha_ds(i))/dL; % 方位二阶倒数
k(i,:) = ((dalpha_ds(i))^2+(sin(alpha(i)))^2*(dpha_ds(i))^2)^0.5; % 井眼曲率
dk_ds(i,:) = (k(i-1) - k(i))/dL; % 井眼曲率变化
kal(i,:) = dalpha_ds(i); % 井斜平面井眼曲率
kph(i,:) = sin(alpha(i))*dpha_ds(i); % 方位平面井眼曲率
tz(i,:) = cos(alpha(i)); % 切线方向投影分量
if dalpha_ds(i) == 0 && dpha_ds(i) ~= 0
nz(i,:) = 0;
bz(i,:) = sin(alpha(i));
elseif dpha_ds(i) == 0 || (dalpha_ds(i) == 0 && dpha_ds(i) == 0)
nz(i,:) = sin(alpha(i));
bz(i,:) = 0;
else
nz(i,:) = dalpha_ds(i)*sin(alpha(i))/k(i); % 主法线方向投影分量
bz(i,:) = -dpha_ds(i)*(sin(alpha(i)))^2/k(i); % 次法线方向投影分量
end
if k(i)<1e-6
k(i,:)=1e-6; % 如果井眼曲率k小于1e-6，为计算准确，令k=1e-6
end
tau(i,:) = (sin(alpha(i))*(dalpha_ds(i)*dpha_ds2(i)-dpha_ds(i)*dalpha_ds2(i))+dpha_ds(i)*((dalpha_ds(i))^2+(k(i))^2)*cos(alpha(i)))/(k(i))^2; % 井眼挠度
end
end