光纤陀螺仪光纤尾纤装配动力学建模仿真与应力计算方法

　　光纤陀螺仪(FOG)中光纤尾纤的装配主要依靠工人手工完成，装配应力控制凭借人工经验，一致性难以保障。针对光纤装配应力定量控制的难题，提出一种基于离散微分几何理论的光纤动力学建模仿真方法。首先，基于离散微分几何理论建立了光纤运动学框架，并由弹性势能推导了包含拉伸、弯曲和扭转变形的光纤动力学方程以及应力计算方法。其次，提出了基于Newmark隐式算法的光纤动力学方程数值求解方法，并设计了光纤在装配过程中与周围物体的碰撞检测与响应策略。随后，进行了光纤圆周布设和弯扭变形的仿真测试。结果表明，光纤控制端与固定端之间的自由长度越短，光纤圆周布设的位置越精确;光纤截面半径越大，在同样的扭转角下光纤产生的弯扭变形程度越大。最后，设计了一种光纤陀螺仪装配布局，并通过动力学仿真得到了光纤装配过程的应力分布及最大应力变化，为陀螺仪光路低应力装配与路径规划提供了理论支撑。

　　关键词

　　光纤陀螺仪;光纤尾纤;动力学模型;离散微分几何;应力计算

　　论文《光纤陀螺仪光纤尾纤装配动力学建模仿真与应力计算方法》发表在《计算机集成制造系统》，版权归《计算机集成制造系统》所有。本文来自网络平台，仅供参考。

　　0 引言

　　光纤陀螺仪是一种基于萨格纳克效应的角速度传感器，是惯性导航系统的核心部件之一，它以全固态、高精度、高稳定和低成本等优点在航空航天、船舶及军事等领域得到广泛应用[1]。光路装配是光纤陀螺仪装配的核心环节，光路装配对象包括5大光学器件、光纤、熔接点等，光路装配质量的好坏直接关系到光纤陀螺仪的功能与输出性能参数。其中光纤尾纤装配是光路装配的重要组成步骤，光纤尾纤的装配质量对于中高精度光纤陀螺仪最终性能的影响不容忽视。随着中高精度光纤陀螺仪应用范围越来越广泛，对光纤尾纤装配质量和一致性控制的研究显得日益重要。

　　光纤尾纤装配质量对于光纤陀螺仪精度和性能的影响主要体现在光纤弯曲应力、扭曲应力造成的损耗及光波偏振态变化上[2-3]。在目前光纤陀螺仪光路的应力研究中，学者们大多是针对光纤环中的光纤，通过应力分析仪来进行光纤应力的分析和测试[4-6]，较少涉及光路装配过程中光纤尾纤应力状态的评估。而光纤尾纤在光纤陀螺仪中所受到的扭曲应力和弯曲应力与光路装配路径有很大关系。

　　另一方面，光纤陀螺仪中光纤本身由纤芯和包层组成，涉及到复杂的材料和力学特性，目前研究者主要通过将光纤简化为弹性变形进行建模分析[7-8]。而近年来，在弹性线体的变形预测中，通过建立包含重力、拉压、弯曲和扭转等特性的物理模型来进行动态仿真已成为主要研究方法[9-11]。Cosserat弹性杆模型[12-13]包含了弹性细杆的轴向线应变及弯曲剪应变等物理参数，具有更精确的仿真计算，但涉及到求解非线性偏微分方程难以获得精确解、离散化后参数过多等问题。考虑到模型仿真的真实性和实时性，离散微分几何理论被用于处理连续线体的离散化建模中[10,14]，通过保留表征实际(光滑)物理系统的关键几何结构来对连续线体进行离散化，在保证连续系统的鲁棒性的同时极大地简化了计算量[15]。其中，BERGOU等[16]通过一种离散处理的自适应框架曲线，建立了弹性细杆的离散几何模型，并模拟了弹性细杆的弯曲、扭转以及打结现象。LV等[17]利用离散微分几何理论建立了柔性线缆的离散弹性细杆动力学模型，并采用中心线-角度的表示方法，提出了离散弹性细杆模型的动力学求解方法，从而对柔性线缆进行实时动态变形模拟。目前，离散微分几何理论已经广泛地应用于具有复杂的线体或类线体的场景建模中，特别是在计算机图形学中[18-20]。

　　因此，本文提出一种基于离散微分几何理论的光纤动力学建模仿真方法，研究光纤陀螺仪装配过程中光纤变形与应力变化，用于建立光纤尾纤应力与布设路径之间的关系，以规划最佳光纤布设路径。该方法采用自适应参考框架与中心线结合的方式，描述光纤装配过程中包括拉伸、弯曲和扭转变形的运动学位姿，并通过动力学模型快速计算光纤装配过程中的应力变化。

　　1 光纤动力学模型

　　1.1 运动学框架

　　在光纤陀螺仪中，光纤尾纤的几何表示是建立光纤动力学模型的关键问题。本文假设光纤由一条中心线和一系列均匀各向同性的圆形截面组成，其中，中心线穿过每个圆形截面的几何中心，中心线是一条空间光滑曲线且至少具有二阶可微性，中心线上的截面可以相对移动和旋转。

　　在全局坐标系中(Oξηζ)，光纤的中心线C可以用向量r(s, t)表示，其中s是以中心线的一端为原点建立的弧坐标系，t为时间。则中心线的运动学可由其运动速度$dot{r}(s, t)=partial r / partial t$与变形梯度即中心线的切向量$r'(s, t)=partial r / partial s$表示。为描述光纤的弯扭状态，定义一个与刚性横截面固连的材料框架{d₁(s, t), d₂(s, t), d₃(s, t)}，3个单位向量dₖ(k=1,2,3)两两互相垂直，且d₃=d₁×d₂是横截面的单位法向量。

　　对于细长光纤，通常忽略其剪切变形，这时光纤横截面始终与中心线正交，即横截面单位法向量d₃与中心线切线方向的单位矢量t重合：

　　[d_{3}=t=frac{r'}{left|r' ight|}]

　　则光纤的拉伸应变εᵣ可表示为：

　　[varepsilon_{r}=r' cdot d_{3}=left|r' ight|]

　　而根据无穷小旋转的定义，材料框架{d₁, d₂, d₃}在空间上沿弧长的变化dₖ'(k=1,2,3)和在时间上的变化$dot{d}_{k}$(k=1,2,3)为：

　　[d_{k}'=omega × d_{k}, dot{d}_{k}=Omega × d_{k}]

　　其中：ω表示光纤的弯扭应变即弯扭度矢量，Ω表示光纤的弯扭速度，即角速度矢量：

　　[omega=d_{3} × d_{3}'+left(d_{1}' cdot d_{2} ight) d_{3}]

　　[Omega=d_{3} × dot{d}_{3}+left(dot{d}_{1} cdot d_{2} ight) d_{3}]

　　在不考虑剪切应变的情况下，弯扭度矢量ω与角速度矢量Ω可表示为：

　　[omega=t × t'+left(d_{1}' cdot d_{2} ight) t=K b+omega_{3} t]

　　[Omega=t × dot{t}+left(dot{d}_{1} cdot d_{2} ight) t=t × dot{t}+Omega_{3} t]

　　其中：Kb=t×t'为中心线副法线曲率，K=|t'|为曲率大小，b=t×t'/|t'|为副法线方向的单位矢量，ω₃=d₁'·d₂为扭转应变即εₜ，Ω₃=$dot{d}_{1}·d_{2}$为扭转角速度。则光纤的弯扭应变εᵦ、εₜ可表示为：

　　[varepsilon_{b}=|omega-(omega cdot t) t|=K, varepsilon_{t}=|omega cdot t|=omega_{3}]

　　为了表示光纤材料框架的位置和方向，需要光纤的中心线坐标和材料框架的扭转度。对于材料框架的表示，有四元数[21-22]和自适应框架[14,16]两种典型的方案。中心线与四元数组合的方案需要7个变量，其中包含冗余变量和多余约束，而中心线与自适应框架组合的方案只需要4个变量。因此，采用中心线和自适应框架组合的方案来表示材料框架的位置和方向。中心线的自适应框架{h₁, h₂, h₃}通常可以用Bishop框架即空间平行框架B={b₁, b₂, t}和时间平行框架A={a₁, a₂, t}表示[14]，其中空间平行框架的弯扭度矢量ω_B和时间平行框架的角速度矢量Ω_A分别表示为：

　　[omega_{B}=K b, Omega_{A}=t × dot{t}]

　　进一步将光纤离散为N+2个节点(0,1,⋯,N+1)和N+1个圆柱段(0,1,⋯,N)。初始时刻离散节点均匀分布，节点坐标为xᵢ={xᵢ, yᵢ, zᵢ}，光纤第i条线段的切向量为sⁱ=xᵢ₊₁−xᵢ，对应单位切向量为tⁱ=sⁱ/|sⁱ|，每个光纤线段上材料框架{d₁ⁱ, d₂ⁱ, d₃ⁱ}与时间平行框架A={a₁ⁱ, a₂ⁱ, tⁱ}之间的夹角为θ。离散化后光纤位姿的广义坐标系表示为X，大小为4N+7：

　　[X=left{x_{0}, heta^{0}, x_{1}, heta^{1}, cdots, x_{N}, heta^{N}, x_{N+1} ight}^{T}]

　　离散化后光纤第i条线段的拉伸应变εₑⁱ可表示为：

　　[varepsilon_{e}^{i}=frac{left|s^{i} ight|}{left|ar{s}^{i} ight|}]

　　其中|sⁱ|和|$ar{s}$ⁱ|分别为第i条线段的当前长度和初始长度。

　　第i个节点的弯扭应变εᵦⁱ、εₜⁱ可表示为：

　　[varepsilon_{i}^{b}=K_{i}=left|t_{i}' ight|=frac{kappa_{i}}{ar{l}_{i}}, varepsilon_{i}^{t}=frac{alpha_{i}}{ar{l}_{i}}]

　　其中：$ar{l}_{i}$=(|$ar{s}$ⁱ⁻¹|+|$ar{s}$ⁱ|)/2为第i个节点的Voronoi区域长度，即第i-1条线段和第i条线段长度和的一半;αᵢ为第i个节点的离散扭转角;Kᵢ为第i个节点的副法线曲率;κᵢ为第i个节点的离散积分曲率，定义为：

　　[kappa_{i}=2 tan frac{varphi_{i}}{2},(kappa b)_{i}=frac{2 t^{i-1} × t^{i}}{1+t^{i-1} cdot t^{i}}]

　　其中φᵢ为相邻线段i-1和i之间的夹角。

　　为了计算离散扭转角αᵢ，引入了平行传递的概念，即向量c₁到c₂的最小旋转为：

　　[P_{c_{1}}^{c_{2}}=Rleft(c_{1} × c_{2}, angleleft(c_{1}, c_{2} ight) ight)]

　　具体而言，平行传递指沿着c₁×c₂的方向使向量c₁旋转角度∠(c₁, c₂)。空间平行框架只取决于中心线位置，这也意味着空间平行框架之间的平行传递为：

　　[b_{k}^{i}left(t_{n} ight)=P_{t^{i}}^{t^{i-1}}left(t_{n} ight) b_{k}^{i-1}left(t_{n} ight), k=1,2]

　　将材料框架{d₁, d₂, d₃}和空间平行框架{b₁, b₂, t}从第i-1条线段平行传递到第i条线段得到离散扭转角αᵢ=θ_Bⁱ−θ_Bⁱ⁻¹，其中θ_Bⁱ⁻¹和θ_Bⁱ分别为相邻的第i-1条线段和第i条线段上的材料框架与空间平行框架之间的夹角。类似地，同一线段上的时间平行框架也具有平行传递关系：

　　[a_{k}^{i}left(t_{n} ight)=P_{tleft(t_{n-1} ight)}^{tleft(t_{n} ight)} a_{k}^{i}left(t_{n-1} ight), k=1,2]

　　将时间平行框架的向量a₁ⁱ⁻¹进行空间上的平行传递得到向量Pᵢⁱ⁻¹a₁ⁱ⁻¹，使得向量Pᵢⁱ⁻¹a₁ⁱ⁻¹和b₁ⁱ夹角为βⁱ⁻¹，即与第i-1条线段上的空间平行框架和时间平行框架之间的夹角保持一致。最后第i个节点上的离散扭角计算为：

　　[alpha_{i}= heta_{B}^{i}- heta_{B}^{i-1}= heta^{i}- heta^{i-1}+alpha_{r j}^{i}]

　　其中αᵣⱼⁱ为向量Pᵢⁱ⁻¹a₁ⁱ⁻¹和a₁ⁱ之间的夹角。

　　1.2 动力学模型

　　光纤陀螺仪中的光纤在装配过程中变形大但应变小，可近似视作纯弹性线。基于1.1节中的光纤运动学框架，在不考虑剪切变形时，光纤的弹性势能Eₙ包括拉压弹性势能Eₑ(d)、弯曲弹性势能Eᵦ(d)和扭转弹性势能Eₜ(d)，离散化形式为：

　　[

　　left{

　　egin{array}{c}

　　E_{e(d)}=frac{1}{2} sum_{i=0}^{N} frac{E_{e} A^{i}}{left|ar{s}^{i} ight|}left(left|s^{i} ight|-left|ar{s}^{i} ight| ight)^{2}

　　E_{b(d)}=frac{1}{2} sum_{i=1}^{N} frac{E_{b} I_{i}}{ar{l}_{i}}left(kappa_{i}-ar{kappa}_{i} ight)^{2}

　　E_{t(d)}=frac{1}{2} sum_{i=1}^{N} frac{G I_{p, i}}{ar{l}_{i}}left(alpha_{i}-ar{alpha}_{i} ight)^{2}

　　end{array}

　　 ight.

　　]

　　其中：Eₑ、Eᵦ和G分别为光纤的拉伸模量、弯曲模量和剪切模量(本研究中Eₑ=Eᵦ=E，E为杨氏模量，G=E/2(1+ν)，ν为泊松比);Aⁱ为第i条线段上的横截面积;Iᵢ和Iₚ,ᵢ分别为第i个节点上的截面惯性矩和极惯性矩;$ar{varepsilon}_{i}^{b}$、$ar{varepsilon}_{i}^{t}$、$ar{kappa}_{i}$和$ar{alpha}_{i}$分别为第i个节点上Voronoi区域在初始时刻的弯曲应变、扭转应变、离散积分曲率和离散扭转角。

　　节点力是总弹性势能对位置矢量的偏导[23]，而线段上的扭矩是弹性势能对角度的偏导。第i个节点上的弹性力Fᵈⁱ(Fᵈⁱ=Fₑⁱ+Fᵦⁱ+Fₜⁱ)包括拉伸力Fₑⁱ、弯曲力Fᵦⁱ和扭转力Fₜⁱ，而扭矩Mᵢ可以由扭转弹性势能Eₜ(d)对角度θⁱ的偏导得到。

　　因此，光纤在第i个节点受到的轴向应力σₙⁱ和切向应力σₜⁱ可由节点弹性力Fᵈⁱ得到：

　　[sigma_{i}^{n}=frac{left(F_{i}^{d} ight)^{n}}{A}, sigma_{i}^{t}=frac{left(F_{i}^{d} ight)^{t}}{A}]

　　其中(Fᵈⁱ)ⁿ和(Fᵈⁱ)ᵗ分别为第i个节点的轴向弹性力与切向弹性力，定义第i个节点的轴向矢量(即沿中心线的切向矢量)为相邻线段矢量的角平分线矢量t₁₂=tⁱ⁻¹+tⁱ，节点处的截面面积A=πr²。

　　广义弹性力Fₑ的列向量可表示为：

　　[F_{e}=left{cdots, F_{i, x}^{d}, F_{i, y}^{d}, F_{i, z}^{d}, M_{i}, cdots ight}^{T}]

　　此外，光纤还受到广义外力Fₑₓₜ的作用(例如重力和来自周围环境的接触力及接触力矩等)，以及广义速度阻尼力Fₙ。根据牛顿定律，总的动力学方程为：

　　[M cdot a=F_{e}+F_{n}+F_{e x t}]

　　其中：M为节点质量和转动惯量组成的广义质量矩阵，a为广义加速度的列向量。

　　2 动力学算法及碰撞策略

　　2.1 动力学算法

　　为保证动力学模型求解的高效性与良好的鲁棒性，本文采用Newmark多步时间积分隐式算法[24]求解动力学系统。首先，将动力学方程改写为二阶矩阵形式：

　　[M ddot{U}left(t_{n+1} ight)+D dot{U}left(t_{n+1} ight)+K Uleft(t_{n+1} ight)=F_{e x t}left(t_{n} ight)]

　　其中：U(tₙ₊₁)、$dot{U}(t_{n+1})$和$ddot{U}(t_{n+1})$分别为下一时刻tₙ₊₁的未知广义位移、广义速度和广义加速度;Fₑₓₜ(tₙ)为当前时刻tₙ的已知广义外力;广义位移U包括节点位移和线段扭转角：

　　[U=left{u_{1}, eta^{1}, u_{2}, eta^{2}, cdots, u_{N}, eta^{N}, u_{N+1} ight}]

　　其中uᵢ={uᵢ, vᵢ, wᵢ}为第i个节点的位移，βⁱ为第i条线段的扭转角。因此可以得到广义位移U与广义坐标X之间的关系。D为速度阻尼刚度矩阵(即对角线上元素为阻尼系数的对角矩阵)，K为弹性刚度矩阵，可由广义弹性力Fₑ线性化得到：

　　[K=frac{partial F_{e}}{partial U}, F_{e}=K U+b_{e}]

　　其中bₑ(tₙ)=Fₑ(tₙ)−K U(tₙ)为常数项，最后需要移动到Fₑₓₜ(tₙ)中。

　　在确定好动力学方程的系数矩阵后，进行Newmark算法。首先，计算tₙ₊₁时刻的广义位移和广义速度的预测值：

　　[ ilde{U}left(t_{n+1} ight)=Uleft(t_{n} ight)+dot{U}left(t_{n} ight) Delta t_{n+1}+ddot{U}left(t_{n} ight)left(frac{1}{2}-lambda_{2} ight) Delta t_{n+1}^{2}]

　　[ ilde{dot{U}}left(t_{n+1} ight)=dot{U}left(t_{n} ight)+ddot{U}left(t_{n} ight)left(1-lambda_{1} ight) Delta t_{n+1}]

　　其中取权重系数λ₁=1/2和λ₂=1/4，得到梯形规则下时间步长的无条件稳定性。最后得到tₙ₊₁时刻的广义加速度$ddot{U}(t_{n+1})$、广义速度$dot{U}(t_{n+1})$、广义位移U(tₙ₊₁)为：

　　[

　　egin{gathered}

　　ddot{U}left(t_{n+1} ight)=left(M+D lambda_{1} Delta t_{n+1}+K lambda_{2} Delta t_{n+1}^{2} ight)^{-1} cdot

　　left(F_{e x t}left(t_{n} ight)-D ilde{dot{U}}left(t_{n+1} ight)-K ilde{U}left(t_{n+1} ight) ight)

　　end{gathered}

　　]

　　[dot{U}left(t_{n+1} ight)= ilde{dot{U}}left(t_{n+1} ight)+ddot{U}left(t_{n+1} ight) lambda_{1} Delta t_{n+1}]

　　[Uleft(t_{n+1} ight)= ilde{U}left(t_{n+1} ight)+ddot{U}left(t_{n+1} ight) lambda_{2} Delta t_{n+1}^{2}]

　　2.2 碰撞检测与响应

　　在光纤陀螺仪的光纤装配过程中存在着光纤与其他物体之间的接触碰撞，为了保证动力学模型仿真的真实性，需要进行光纤的碰撞检测与响应策略研究。

　　在本研究中，采用杨啸东等[25]提出的分层碰撞检测方法，即第一层为包围光纤整体的最小圆柱体包围盒，第二层为包围光纤线段的最小球体包围盒，第三层为包围光纤线段的最小圆柱体包围盒。当检测到第一层发生碰撞时，执行第二层碰撞检测;当检测到第二层发生碰撞时，执行第三层碰撞检测。此外，对于相邻的光纤线段认为不存在碰撞情况。

　　对于光纤与周围物体之间的碰撞响应(如光纤装配过程与陀螺仪结构面的接触碰撞)，采用一种接触力的响应策略。当光纤第i个节点发生碰撞时，节点会受到碰撞物体的接触力Fᵢ，包括法向支撑力Fₙⁱ和切向摩擦力Fբⁱ：

　　[F_{i}=F_{i}^{n}+F_{i}^{f}]

　　其中支撑力Fₙⁱ使得节点沿接触面法向的受力和速度变为零(假设碰撞过程为完全非弹性碰撞)，则碰撞后的法向支撑力Fₙⁱ和节点切向速度vₜₐₙⁱ为：

　　[F_{i}^{n}=-left(F_{i}^{g}+F_{i}^{d} ight)^{n}]

　　[v_{i}^{t a n}=v_{i}-v_{i}^{n}]

　　其中：F₉ⁱ为节点重力，Fᵈⁱ为节点弹性力，vᵢ为节点碰撞前的速度，vₙⁱ为节点法向速度。对于切向摩擦力Fբⁱ，若节点碰撞后速度vₜₐₙⁱ不为零或碰撞后节点所受切向力Fₜₐₙⁱ大于最大静摩擦力Fₘₐₓʲˢ(其中：Fₘₐₓʲˢ=kₛFₙⁱ，kₛ为静摩擦系数)，则节点受到切向的滑动摩擦力：

　　[F_{i}^{f}=F_{i}^{f d}=k_{d}left|F_{i}^{n} ight| frac{v_{i}^{t a n}}{left|v_{i}^{t a n} ight|}]

　　其中k_d为动摩擦系数。若节点碰撞后速度vₜₐₙⁱ为零且碰撞后节点所受切向力Fₜₐₙⁱ小于最大静摩擦力Fₘₐₓʲˢ，则节点将保持静止，对应的切向摩擦力为：

　　[F_{i}^{f}=F_{i}^{f s}=-F_{i}^{t a n}]

　　3 光纤动力学仿真

　　本章主要基于光纤动力学模型的MATLAB代码进行数值仿真分析，首先检验了光纤动力学模型的收敛性，然后基于该模型进行了光纤圆周布设及弯扭的仿真分析，最后设计了一种光纤陀螺仪装配结构，并通过光纤的动力学仿真来分析光纤沿轴向与切向的应力分布及最大应力变化。光纤这类单元尺寸较小、刚度较大的柔性体在有实时性要求的动力学仿真情况下容易出现稳定性问题，而装配仿真属于实时性要求较高的场合。考虑到本章主要探究装配过程中布设路径对光纤应力的影响，且弹性范围内材料的弹性模量只影响光纤应力的大小，不影响光纤的应力与布设路径之间的相对关系，因此本节采用弹性模量较低的塑料光纤的模量进行仿真分析，本章所有仿真案例中涉及的主要参数如表1所示。

　　表1 光纤动力学仿真参数

　　| 参数名称 | 参数值 |

　　| 密度ρ_f/(kg/m³) | 0.32×10³ |

　　| 移动速度v_f/(m/s) | 0.01 |

　　| 截面半径r_f/m | 1×10⁻³ |

　　| 杨氏模量E_f/Pa | 1×10⁵ |

　　| 泊松比λ_f | 0.3 |

　　| 阻力系数 | 1×10⁻⁴ |

　　| 静摩擦系数 | 0.3 |

　　| 动摩擦系数 | 0.2 |

　　3.1 仿真收敛性检验

　　在进行光纤动力学仿真分析前，需要对动力学模型的网格尺寸(光纤离散化的每条线段长度)与时间步长进行收敛性检验。以纯弯曲仿真为例，研究不同网格尺寸与时间步长下仿真结果的形状变化与应力分布。具体过程为：一根长度L_f=0.4m的光纤悬空并水平放置，然后控制光纤两端以固定的移动速度相向移动，仿真时间为10s，其他相关的仿真参数采用表1中的参数值。

　　对于网格尺寸的收敛性研究，固定时间步长Δt为0.001s，网格尺寸分别取0.04m、0.02m和0.01m，其他仿真参数保持一致。结果显示，随着网格尺寸的等比减小，中心线形状差异越来越小，且应力分布的变化逐渐收敛。考虑到计算精度与效率，选取0.01m的网格尺寸作为后续仿真参数。

　　对于时间步长的收敛性研究，固定网格尺寸为0.01m，时间步长Δt分别取0.01s、0.001s和0.0001s，其他仿真参数保持一致。结果表明，时间步长对中心线形状及应力分布的影响很小，考虑到计算精度与效率，选取0.001s的时间步长作为后续的仿真参数。

　　3.2 仿真案例测试

　　光纤装配过程中光纤自由部分的长度(即光纤控制端与光纤固定端之间的距离)会影响装配过程中光纤尾纤布设的位置精度。为探究不同自由长度下光纤装配精度的差异，进行光纤圆周布设的仿真分析。XY平面表示光纤陀螺仪中的一个用于固定光纤的结构平面，光纤固定端定义为固定在XY平面上的光纤端部，光纤控制端定义为沿着半圆周轨迹移动的光纤端部，布设过程中光纤不断下落至平面上并保持光纤控制端距XY平面的高度H不变。分别取高度H=0.04m、0.008m和0.002m，进行3种不同高度下的圆周布设仿真。其中，光纤控制端的水平移动速度v_xy=0.01m/s，竖直下落速度v_z=0.01m/s，圆周半径R=0.1m，其它相关的仿真参数采用表1中的参数值。结果显示，圆周布设的高度H越小，光纤圆周布设越来越接近目标圆周曲线。因此，为保证光纤尾纤在装配过程中具有较好的布设位置精度，光纤控制端与固定端的距离应尽可能短，以便更好地控制光纤布设。

　　为研究光纤的弯扭特性，进行弯扭仿真测试分析。结果显示了光纤在不同的圆周切向扭转角度时的空间构型的变化，同一光纤在固定端和控制端相对扭转变形一周与四周的仿真模拟中，光纤长度L₁取为0.4m，光纤截面半径r_f取为1×10⁻³m，其他相关的仿真参数采用表1中的参数值，仿真过程分别为光纤两端相对扭转一周(总扭转角度为360°)和四周(总扭转角度为1440°)后再逐渐沿着水平直线分别靠近移动0.15m和0.06m(每端移动分别为0.15m和0.06m，两端总靠近位移分别为0.3m和0.12m)。

　　在同一扭转角度下，不同半径的光纤会产生不同的空间变形效果。在光纤两端相对扭转四周、相向位移为0.12m的情况下，光纤半径r_f=0.002m时可以实现光纤空间螺旋变形两周的效果，但通过仿真发现当光纤半径r_f=0.001m时不能实现光纤空间螺旋变形两周的效果。进一步分析不同光纤半径下的拉伸力Fₑ、弯曲力Fᵦ和扭转力Fₜ作用，发现在光纤半径r_f=0.001m时，沿光纤中心线分布的拉伸力Fₑ占主导作用，而弯曲力Fᵦ和扭转力Fₜ几乎为零，从而使得光纤无法产生较大的扭转变形;在光纤半径r_f=0.002m时，沿光纤中心线分布的扭转力Fₜ占主导作用，其次是弯曲力Fᵦ与拉伸力Fₑ，这也使得光纤受到扭转力的作用较大，从而产生较大的扭转变形。因此，在光纤长度和扭转角度一定的情况下，光纤半径越大，所产生的扭转力也越大，光纤越容易扭转变形。

　　3.3 光纤陀螺仪光路应力分析

　　在光纤陀螺仪装配过程中，光纤的盘纤过程是非常关键的环节。由于装配过程的复杂性，其装配工艺参数及盘纤路径的设计目前仍主要基于人工手动装配的经验积累，缺乏定量评估光纤尾纤应力和光纤布设路径的方法。本节基于光纤动力学模型对光纤装配过程进行动力学仿真，研究装配过程中光路的应力分布及变化。

　　一种光纤陀螺仪光路装配中探测器的光纤尾纤盘纤过程的切向应力云图，基于C++语言编程开发可视化操作界面进行仿真操作与显示，并采用三维造型引擎ACIS和三维显示交互工具包HOOPS建立三维环境。光纤装配过程的变形及光纤沿轴向和切向的应力分布首先由光纤动力学模型的MATLAB代码仿真计算得到，然后通过C++与MATLAB的数据交互实现光纤在可视化操作界面上的装配过程与应力云图显示。其中，探测器及周围结构为Solidworks中的三维模型，其凹槽内径D=0.12m，凹槽深度H₁=0.02m，主要仿真参数采用表1的参数值，对探测器的光纤尾纤盘纤过程进行装配仿真，盘纤过程从探测器管壳端的光纤开始，在盘纤过程中光纤控制端高度H₂=0.04m保持不变。

　　为进一步分析光纤陀螺仪中光路应力分布及变化，对探测器光纤尾纤作轴向应力和切向应力沿中心线长度分布的变化曲线。结果显示，光纤最大的轴向应力σⁿ和切向应力σᵗ均出现在光纤尾纤的端部。其中，光纤最大切向应力远远大于最大轴向应力，且光纤固定端的切向应力远大于控制端。这是由于盘纤过程中光纤主要发生弯扭变形，而光纤受到的拉伸作用很小，且相较于光纤固定端，控制端的光纤弯扭程度相对较小。此外，由于盘纤过程近似为圆周轨迹运动且每一圈的圆周半径近似相等，光纤中间部分的轴向和切向应力较小且变化较均匀。

　　在光纤装配过程中，光路的最大应力也在随时间变化。结果表明，最大轴向应力与最大切向应力在最开始盘纤过程中变化较大，随后逐渐稳定，但最大轴向应力相较于最大切向应力波动较大。这是因为光纤在刚开始盘纤时固定端附近的光纤弯扭变形较大，稳定布设后光纤弯扭变形较均匀，而光纤拉伸变形在光纤固定端与控制端附近变形均较大，所以在光纤稳定布设过程中最大轴向应力受光纤控制端影响较大。此外，整个盘纤过程中切向应力的最大值出现在盘纤最开始的时刻(即仿真时间t=0s)，对应于初始状态的切向应力云图，表明此时光纤弯扭变形程度最大;而轴向应力的最大值出现在盘纤开始一段时间后(即仿真时间t=1.2s)，对应于该时刻的轴向应力云图，表明此时光纤拉伸变形程度最大。这也为未来进一步定量设计最小应力的光纤盘纤路径提供了借鉴意义。

　　4 结束语

　　本文针对光纤陀螺仪中光纤尾纤装配应力定量控制的难题，以光纤陀螺仪的光纤尾纤为研究对象，对其装配过程中的变形及应力变化进行了动力学建模与仿真，主要研究结论如下：

　　1. 基于离散微分几何理论建立了光纤运动学框架，该方法能在简化光纤自由度、降低计算量的同时，保证计算结果具有良好的鲁棒性与精度。随后，基于该运动学框架由弹性势能推导了光纤动力学方程，得到了描述光纤装配过程中拉伸、弯曲和扭转变形以及光纤应力分布的动力学模型。

　　2. 基于Newmark隐式算法提出了一种求解光纤动力学模型的数值方法，并设计了光纤装配过程中的碰撞检测与响应策略，实现了光纤装配仿真的高效稳定计算与物理意义上的接触碰撞响应效果。

　　3. 通过对光纤圆周布设的仿真发现，在光纤装配过程中，光纤控制端与固定端之间自由长度越短，光纤布设的位置越精确;在光纤弯扭仿真过程中，同样的光纤长度与扭转角度下，光纤半径越大使得扭转力作用越大，光纤扭转变形的程度越大。

　　4. 基于光纤动力学模型对一种光纤陀螺仪的盘纤过程进行了动力学仿真，得到了光纤盘纤过程中光路应力分布与最大应力变化曲线，确定了光纤轴向应力与切向应力最大时的光纤布设路径，实现了光纤布设路径下的轴向应力与切向应力的定量分析。

　　本文通过建立光纤动力学模型，对光纤陀螺仪装配过程中光路应力分布与变化进行了研究，为光纤陀螺仪低应力装配的定量控制提供了理论支撑。未来也将进一步研究小尺寸大刚度下光纤装配过程的动力学仿真，并开发光纤装配仿真系统，实现虚拟环境下光纤低应力装配的实时仿真与路径规划。

　　参考文献

　　[1] WANG Wei. Interferometric fiber optic gyroscope technology[M]. Beijing: China Aerospace Press, 2010 (in Chinese). [王巍. 干涉型光纤陀螺仪技术[M]. 北京: 中国宇航出版社, 2010.]

　　[2] WU Changxin, YANG Mingwei, YANG Yuanhong, et al. Effect of fiber bend on interferometric fiber optical gyroscope performance[J]. Acta Optica Sinica, 2014, 34(3): 87-91 (in Chinese). [吴长莘, 杨明伟, 杨远洪, 等. 光纤弯曲对干涉式光纤陀螺性能的影响[J]. 光学学报, 2014, 34(3): 87-91.]

　　[3] WU Weifeng. Research on the bending loss of panda polarization maintaining fiber[D]. Nanjing: Nanjing University of Science and Technology, 2009 (in Chinese). [吴卫锋. 熊猫型保偏光纤弯曲损耗的研究[D]. 南京: 南京理工大学, 2009.]

　　[4] YANG Yuanhong, YI Xiaosu, MENG Zhaokui. Experimental study on strain distribution in fiber coil used in fiber optic gyroscope[J]. Piezoelectrics and Acoustooptics, 2005, 27(2): 98-101 (in Chinese). [杨远洪, 伊小素, 孟照魁. 光纤陀螺用光纤环的应力分布实验研究[J]. 压电与声光, 2005, 27(2): 98-101.]

　　[5] YIN Qiqi, LI Jing, ZHANG Zhihua. Experimental study on the stress in fiber sensing coil[J]. Navigation and Control, 2016, 15(2): 79, 98-103 (in Chinese). [尹其其, 李晶, 张智华. 光纤环内部应力试验研究[J]. 导航与控制, 2016, 15(2): 79, 98-103.]

　　[6] QIU Hongfang. Analysis on influence of inflection point in optical fiber loop on performance of optical fiber gyroscope[J]. Tianjin Science & Technology, 2023, 50(5): 23-25 (in Chinese). [邱红芳. 光纤环圈拐点对光纤陀螺性能影响的分析[J]. 天津科技, 2023, 50(5): 23-25.]

　　[7] YANG Xuewei, DU Xinzheng, LIN Heng, et al. Thermal-stress effect of package optical fiber coil adhesive on optical fiber[J]. Journal of Chinese Inertial Technology, 2003(6): 95-98 (in Chinese). [杨学围, 杜新政, 林恒, 等. 封装光纤线圈的胶粘剂对光纤产生的热应力影响[J]. 中国惯性技术学报, 2003(6): 95-98.]

　　[8] TANG Xu, ZHU Yazhou, SHEN Zhongxiang, et al. Research on application of optical fiber sensing material in ship structural health monitoring[J]. Ship Science and Technology, 2020, 42(23): 181-185 (in Chinese). [唐旭, 朱亚洲, 沈中祥, 等. 光纤传感材料在船舶结构健康监测中的应用[J]. 舰船科学技术, 2020, 42(23): 181-185.]

　　[9] LV N J, LIU J H, XIA H X, et al. A review of techniques for modeling flexible cables[J]. Computer-Aided Design, 2020, 122: 102826. DOI: 10.1016/j.cad.2020.102826.

　　[10] JAWED M K, NOVELIA A, O'REILLY O M. A primer on the kinematics of discrete elastic rods[M]. Berlin, Germany: Springer-Verlag, 2018.

　　[11] ZHANG Y, CHEN D S, QIAN H L. Computational method for the deformation mechanism of non-prestressed cable net structures based on the vector form intrinsic finite element method[J]. Engineering Structures, 2021, 231: 111788. DOI: 10.1016/j.engstruct.2020.111788.

　　[12] LIN Haili, LIU Jianhua, TANG Chengtong, et al. Physical characteristic oriented modeling for flexible cable harness based on Cosserat elastic rod theory[J]. Journal of Graphics, 2016, 37(1): 34-42 (in Chinese). [林海立, 刘检华, 唐承统, 等. 基于Cosserat弹性杆理论的柔性线缆物理建模方法[J]. 图学学报, 2016, 37(1): 34-42.]

　　[13] LIU L, YAN Z, HE Y M, et al. On the configuration evolution of soft filaments under combined tension and torsion[J]. Acta Mechanica Sinica, 2023, 39(11): 222498. DOI: 10.1007/s10409-023-22498-x.

　　[14] BERGOU M, AUDOLY B, VOUGA E, et al. Discrete viscous threads[J]. ACM Transactions on Graphics, 2010, 29(4): 116. DOI: 10.1145/1778765.1778853.

　　[15] DESBRUN M, POLTHIER K, SCHRÖDER P, et al. Discrete differential geometry: An applied introduction[EB/OL]. (2006-05-11) [2023-11-29]. http://geometry.caltech.edu/pubs/GSD06.pdf.

　　[16] BERGOU M, WARDETZKY M, ROBINSON S, et al. Discrete elastic rods[J]. ACM Transactions on Graphics, 2008, 27(3): 63. DOI: 10.1145/1399504.1360662.

　　[17] LV N J, LIU J H, JIA Y Y. Dynamic modeling and control of deformable linear objects for single-arm and dual-arm robot manipulations[J]. IEEE Transactions on Robotics, 2022, 38(4): 2341-2353.

　　[18] AUDOLY B, CLAUVELIN N, BRUN P T, et al. A discrete geometric approach for simulating the dynamics of thin viscous threads[J]. Journal of Computational Physics, 2013, 253: 18-49.

　　[19] SPERL G, SÁNCHEZ-BANDERAS R M, LI M, et al. Estimation of yarn-level simulation models for production fabrics[J]. ACM Transactions on Graphics (TOG), 2022, 41(4): 65. DOI: 10.1145/3528223.3530167.

　　[20] ARMANINI C, BOYER F, MATHEW A T, et al. Soft robots modeling: A structured overview[J]. IEEE Transactions on Robotics, 2023, 39(3): 1728-1748.

　　[21] GREGOIRE M, SCHÖMER E. Interactive simulation of one-dimensional flexible parts[J]. Computer-Aided Design, 2007, 39(8): 694-707.

　　[22] LANG H, LIN N J, ARNOLD M. Multi-body dynamics simulation of geometrically exact Cosserat rods[J]. Multibody System Dynamics, 2011, 25(3): 285-312.

　　[23] LOOCK A, SCHÖMER E, STADTWALDI. A virtual environment for interactive assembly simulation: From rigid bodies to deformable cables[EB/OL].[2023-11-29]. https://www.staff.uni-mainz.de/schoemer/publications/SCI01.pdf.

　　[24] NEWMARK N M. A method of computation for structural dynamics[J]. Journal of the Engineering Mechanics Division, 1959, 85(3): 67-94.

　　[25] YANG Xiaodong, LIU Jianhua, ZHAO Yingfeng, et al. Rigid and flexible hybrid assembly process simulation technology for complex mechatronic products[J]. Computer Integrated Manufacturing Systems, 2019, 25(2): 340-349 (in Chinese). [杨啸东, 刘检华, 赵瑛峰, 等. 面向复杂机电产品的刚柔混合装配过程仿真技术[J]. 计算机集成制造系统, 2019, 25(2): 340-349.]