光纖陀螺儀(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)和在時間上的變化$dottl7lrfpl_{k}$(k=1,2,3)為：

　　[d_{k}'=omega × d_{k}, dottl7lrfpl_{k}=Omega × d_{k}]

　　其中：ω表示光纖的彎扭應變即彎扭度矢量，Ω表示光纖的彎扭速度，即角速度矢量：

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

　　[Omega=d_{3} × dottl7lrfpl_{3}+left(dottl7lrfpl_{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(dottl7lrfpl_{1} cdot d_{2} ight) t=t × dot{t}+Omega_{3} t]

　　其中：Kb=t×t'為中心線副法線曲率，K=|t'|為曲率大小，b=t×t'/|t'|為副法線方向的單位矢量，ω?=d?'·d?為扭轉應變即ε?，Ω?=$dottl7lrfpl_{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}^tl7lrfpl ight)^{n}}{A}, sigma_{i}^{t}=frac{left(F_{i}^tl7lrfpl ight)^{t}}{A}]

　　其中(F??)?和(F??)?分別為第i個節點的軸向彈性力與切向彈性力，定義第i個節點的軸向矢量(即沿中心線的切向矢量)為相鄰線段矢量的角平分線矢量t??=t??¹+t?，節點處的截面面積A=πr²。

　　廣義彈性力F?的列向量可表示為：

　　[F_{e}=left{cdots, F_{i, x}^tl7lrfpl, F_{i, y}^tl7lrfpl, F_{i, z}^tl7lrfpl, 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}^tl7lrfpl 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_tl7lrfplleft|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.]