求解具有規定測地曲率曲線的高效幾何方法

緱建杰，張 鵬，楊 揚，陳代鑫

（1.成都飛機工業（集團）有限責任公司檢驗檢測部，四川 成都 610091 2.長沙理工大學汽車與機械工程學院，湖南 長沙 410114）

1 引言

曲面上的曲線在曲面裁剪[2]、曲面過渡[2]、數控加工刀具路徑規劃[3]等眾多應用中發揮著重要的作用。在一些應用場景中，需要設計出的曲線上測地曲率滿足一定的要求，這種曲線被稱為具有規定測地曲率的曲線。

例如，在自動鋪帶工藝中，鋪放路徑上測地曲率的大小直接決定了復合材料預浸帶的面內變形[4]。由于預浸帶面內變形能力很小，因此帶料通常沿著模具曲面上的測地線鋪放[5]。

測地線是具有規定測地曲率曲線中的一個特例，其上各處測地曲率均為零。自動鋪帶軌跡規劃中經常要處理測地線求解的初值問題，即給出鋪放起始點位置以及初始鋪放方向，求解曲面上的唯一一條測地線。

到目前為止，該問題已經有很多經典的求解方法，具體可分為數值解法[6-7]、離散解法[8-9]以及幾何解法[10-11]。

相較于測地線，本研究主要關注參數曲面上測地曲率不為零的曲線，這類曲線在圖案設計、齒輪設計、機械結構設計等場合均有所應用。

文獻[12]提出了一種保測地曲率的曲面曲線設計方法，其目的是在光滑曲面上進行圖案設計。

在該項研究中，所需的曲面圖案首先被構造為一平面圖案。在曲面上指定起始點和方向后，該平面圖案被逐段映射到目標曲面上。

映射過程中，為保證視覺效果，要求曲面曲線的測地曲率與平面曲線的曲率相等。

文獻[13]基于保測地曲率的曲線設計方法對非圓錐齒輪傳動進行了設計與分析，利用球面上具有規定測地曲率的曲線生成了非圓錐齒輪的空間三維齒廓。

文獻[14]則基于球面上具有規定測地曲率的曲線設計了一種不等速行星輪系分插機構。

具有規定測地曲率的曲線也常應用于復合材料自動化加工領域中。

文獻[15]提出了一種針對圓錐面模具的自動鋪絲軌跡規劃方法，利用等測地曲率曲線實現了纖維變形與纖維方向之間的平衡。

該方法所依據的幾何學原理是可展曲面上的等測地曲率曲線實際上是其對應展開面上的一段圓弧。

文獻[16]給出了一種針對自由曲面模具的自動鋪帶軌跡規劃方法，利用具有微小測地曲率的曲線實現了預浸帶在模具曲面上的轉向，消除了相鄰帶料間過大的鋪放間隙。

本研究提供一種參數曲面上具有規定測地曲率曲線的求解方法，當給出目標曲面上的起始點、方向以及測地曲率的分布時，可追蹤出符合設計要求的空間曲線。

其主要由以下幾個部分組成：第二節對現有的計算方法進行了分析；第三節詳細介紹了提出方法的原理；第四節將所提出的方法與現有的方法在計算效率上進行了比較；第五節展示了本方法的兩種應用場景；第六節則進行總結并給出結論。

2 研究現狀與分析

在設計制造的一些應用場景中，需要設計出的曲面曲線上測地曲率滿足一定的要求，并且測地曲率大小不為零。然而，現有的研究主要處理測地線的計算問題，其中多數方法無法直接用于構造具有非零測地曲率的曲線，例如數值解法以及離散解法這兩類。

到目前為止，具有非零測地曲率曲線的計算方法主要可分為數值解法、分析解法以及幾何解法這三類。

文獻[16]基于微分幾何建立了具有規定測地曲率曲線的數學模型，該數學模型由一個微分方程組表達，需要利用數值積分器進行求解。

其中曲線的設計精度與數值積分器的階數直接相關。文獻[15]提出的分析解法可以構造出準確的目標曲線，然而該方法只適用于可展曲面，應用范圍較窄。

文獻[12]提出的幾何解法在每一步計算中需要將點投影回曲面，降低了計算效率。

文獻[11]提出了一種構造測地線的幾何方法，經過調整，該方法可用于求解具有非零測地曲率的曲線。該方法則直接在參數空間中進行計算，無需額外的投影操作。

現有的幾何方法主要基于伏雷內公式與泰勒展開設計目標曲線。設Q(s)—目標曲線上一點，Q'(s)—曲線在該點的延伸方向，Q''(s)—曲線在該點的曲率。根據現有的幾何方法可知，Q(s)在目標曲線上臨近點的位置以及方向可表達為：

式中：h—計算步長。

在得到臨近點位置及延伸方向后，可采用公式追蹤其它離散點從而得到整段曲線。

需要注意的是，文獻[12]中提出幾何方法是在三維空間中追蹤曲線，而文獻[11]中的幾何方法則在參數空間中追蹤曲線。

由方程可知，現有的幾何方法可以保證二階位置計算精度與一階方向計算精度。

設延伸方向Q'(s)計算誤差為Ah2，其中A—常數。若忽略曲率Q''(s)的計算誤差，則位置計算誤差主要來源于以下兩個部分：

（1）由于方向計算誤差引起的誤差：

（2）泰勒展開的余項：

顯然，在每一步計算中，由于方向計算誤差引起的誤差ε1相對于泰勒展開的余項ε2不能忽略。

因此，現有的幾何方法雖然具有二階位置計算精度，但計算過程中由于方向計算誤差引起的誤差ε1不斷累積，降低了最終精度。

由前述分析可知，若想獲得較好的曲線設計精度，需保證方向計算精度為高階，而這正是現有的幾何方法難以保證的，其主要原因在于這些方法難以得到比曲率Q''(s)更高階的幾何信息。

對現有的規定測地曲率曲線計算方法的主要特點進行了總結與展示，如表1所示。

表1 現有計算方法的主要特點Tab.1 Main Features of Existing Methods

3 求解方法

本小節介紹文中所提出方法的數學原理。設r(u，v)=[x(u，v)，y(u，v)，z(u，v)]為模具曲面S的參數方程，其中x、y、z是參數u和v的可微函數。

令C是曲面S上一條具有規定測地曲率的弧長參數化曲線，如圖1所示。

圖1 測地曲率與法曲率的定義Fig.1 Definition of Geodesic Curvature and Normal Curvature

圖中：P—曲線C上任意一點；N—曲面S在點P的單位法向量；t—曲線C在P點的單位切向量。它可以用如下方程描述：

式中：s—弧長參數。

令b=N×t，則N，t，b三者彼此正交并構成一個右手坐標系。根據微分幾何，曲線C在P點處的曲率向量k可表示為：

式中：kn—法曲率向量；

kg—測地曲率向量；

kn—曲面S在P點處沿切向量t的法曲率；

kg—曲線C在P點處的測地曲率。

由公式可知，曲線C在P點處的測地曲率kg可通過如下公式求解：

設P0為曲線C的初始點；

C'與P0'為曲線C與點P0在參數空間中的逆映射。

根據微分幾何，曲線C'在P0'處的延伸方向l0=與曲線C在P0處的單位切向量t0之間的關系可通過下式表達：

在方程左右兩側同時點乘列向量[ru(P0')rv(P0')]T得到：

3.1 構建求解位置與方向信息的F算子

令P0為曲面上的正則點，則由公式可得，已知曲線C在三維空間中的單位切向量t0，可以求得C'在參數空間中的延伸方向l0為：

式中：E(P0')，F(P0')，G(P0')—曲面在P0'處的第一類基本量。

由點P0'位置以及延伸方向l0，可以通過一階線性逼近得到P0'臨近點P1'的位置：

在獲得P1'的位置之后，為保證追蹤過程能夠持續進行，需計算曲線C在P1處的單位切向量t1以及C'在P1'處的延伸方向l1。由公式可知，單位切向量t0在P0點相對弧長的變化率可表達為：

在得到單位切向量t1后，可以利用與類似的公式計算l1。

顯然，通過以上步驟獲得的P1'與t1只具有一階精度，在多數場合中直接應用會有較大的計算誤差。為提高計算精度，首先將前述具有一階精度的計算過程整體封裝為一個F算子：

如上式所示，F算子輸入為C'上某一點P'的位置以及C在P處的單位切向量t，輸出為C'在P'的延伸方向l以及單位切向量t在P點相對弧長的變化率

具體對于P0'來說，F算子可表示為：

3.2 通過多點平均提高位置、方向計算精度

為提高P1'與t1的計算精度，可采用與龍格-庫塔法相似的原理，在P0'與P1'之間設置多個中間點，然后對中間點的幾何信息進行加權平均以精確預測曲線C'在P0'附近的幾何形狀。

若需保證P1'與t1具有二階計算精度，則依據二階龍格-庫塔法在P0'與P1'之間設置一個中間點

基于公式計算具有二階精度的P1'與t1，得到：

式中：2P1'與2t1的左上標2表示這些變量具有二階計算精度。

若需保證P1'與t1具有四階計算精度，則依據經典四階龍格-庫塔法在P0'與P1'之間設置三個中間點它們采用如下公式進行計算：

式中：4P1'與4t1的左上標4—這些變量具有四階計算精度。

需要注意的是，迭代過程中方程與方程無法保證2t1與4t1為單位向量，這可通過下式進行調整：

在求得P1'的位置以及曲線C在P1處的單位切向量t1后，可以采用與前述相似的步驟求得后續點P2'、P3'…..，直到曲線抵達曲面邊界或離散點的數量滿足規定要求。

顯然，這些離散點[P0'、P1'、P2'......]定義了參數空間中具有規定測地曲率的曲線C'。

4 驗證

為驗證所提出方法的有效性，本節將該方法與現有的方法在計算精度以及計算時間方面進行比較。

在接下來部分給出了三個算例。每個算例中采用現有的方法計算參數曲面上一系列具有規定測地曲率的曲線，直到曲線抵達曲面邊界。同時，記錄計算過程所消耗的時間。

另一方面，以計算曲線末端點與理想曲線末端點在三維空間中的距離Error作為評價計算精度的指標。

顯然，Error越小表示方法越精確。由于很難通過解析方法獲得一般曲面上理想的具有規定測地曲率的曲線，因此本節中的理想曲線均采用數值解法求得。

為保證理想曲線具有足夠高的精度，數值解法采用了RK4數值積分器以及極小的計算步長。

為便于比較計算精度，計算過程中令各方法的計算步長保持恒定。

文中所有的代碼均在MATLAB 2016b環境中實現，然后運行于個人電腦[Intel Core?i5-8400，16G]上來測試計算效率。

如圖2所示，算例1計算一個2×2B樣條曲面上的等測地曲率曲線。

圖2 2×2B樣條曲面上的等測地曲率曲線Fig.2 Curves with Constant Geodesic Curvature on a 2×2B Spline Surface

計算結果，如圖3、圖4以及表2所示。

圖4 算例1中計算時間對計算精度的影響Fig.4 Influence of Time Costs on Accuracy in Example 1

表2 算例1中各方法所消耗的計算時間Tab.2 Time Costs of Each Method in Example 1

圖3 算例1中計算步長對計算精度的影響Fig.3 Influence of Step Size on Accuracy in Example 1

其中，不同方法在相同計算步長下的計算精度以及計算步長對計算精度的影響，如表3所示。

不同方法在相同時間內所能獲得的精度，如表4所示。計算過程消耗的時間，如表2所示。

在以下圖示中，RK4、RK2表示采用文獻[16]中的數值方法以及RK4、RK2數值積分器進行求解，Hu’s表示采用文獻[12]中的幾何方法進行求解，Zhang’s表示采用文獻[11]中的幾何方法進行求解，Our4、Our2則表示采用提出的四階與二階方法進行求解。

另外，為減小文章篇幅，表2、表3、表4 只列出了部分實驗數據。

表4 算例3中各方法所消耗的計算時間Tab.4 Time Costs of Each Method in Example 3

算例2計算一個雙三次B樣條曲面上的等測地曲率曲線，如圖5所示。

圖5 雙三次B樣條曲面上的等測地曲率曲線Fig.5 Curves with Constant Geodesic Curvature on a 3×3B Spline Surface

各方法在相同計算步長下的精度以及計算步長對精度的影響，如圖6所示。

圖6 算例2中計算步長對計算精度的影響Fig.6 Influence of Step Size on Accuracy in Example 2

它們在相同時間內所能獲得的精度，如圖7所示。

圖7 算例2中計算時間對計算精度的影響Fig.7 Influence of Time Costs on Accuracy in Example 2

計算過程消耗的時間，如表3所示。

表3 算例2中各方法所消耗的計算時間Tab.3 Time Costs of Each Method in Example 2

算例3計算一個6×7B樣條曲面上的等測地曲率曲線，如圖8所示。

圖8 6×7B樣條曲面上的等測地曲率曲線Fig.8 Curves with Constant Geodesic Curvature on a 6×7b Spline Surface

計算結果，如圖9、圖10以及表4所示。

圖9 算例3中計算步長對計算精度的影響Fig.9 Influence of Step Size on Accuracy in Example 3

其中，計算精度隨步長以及計算時間的變化規律，如圖9、圖10所示。計算過程消耗的時間，如表4所示。

圖10 算例3中計算時間對計算精度的影響Fig.10 Influence of Time Costs on Accuracy in Example 3

由圖2、圖5、圖8可知，測地曲率的幅值決定了曲線偏離測地線的程度，測地曲率的正負決定了曲線所處的方位。

具體地說，當測地曲率為正時，曲線位于測地線延伸方向的左側；當測地曲率為負時，曲線位于測地線延伸方向的右側。

首先，考慮計算精度問題。

由圖3、圖6、圖9可看出，隨著計算步長的減小，各方法的計算精度均有所提高。

同時，方法階數越高，收斂速度越快。對于四階方法（RK4與Our4），當計算步長很小時，計算結果容易產生震蕩，此時繼續減小步長對改善計算精度的效果不顯著。由計算結果可看出，現有的幾何方法精度較差，這主要是由于這些方法無法保證較高的方向計算精度。

對于二階方法而言，文中提出的方法在精度上要優于文獻[16]中的數值解法；對于四階方法而言，文中提出的方法在多數情況下可以獲得更好的計算精度。

另外，計算結果表明四階方法的計算精度要顯著高于二階方法，這主要是由于所構造的曲線光滑連續，在解光滑的情況下，高階方法通常優勢明顯。其次，考慮計算時間。由表2、表3、表4可看出，隨著計算步長的減小，計算過程消耗的時間顯著增加。

由上表可看出，與四階方法相比，現有的幾何方法消耗的計算時間更少，這主要是由于這些方法在計算過程中無需進行多點平均、多次迭代。

文獻[12]中的幾何方法計算速度較慢，這主要是由于該方法在每一步計算中需要進行點到曲面的投射操作。通過比較發現，本方法在時間消耗上與現有的數值方法大致相當。

為綜合考慮計算精度與計算時間，圖4、圖7、圖10 展示了不同方法在相同時間內所能獲得的計算精度。結果表明，在相同時間內，本方法獲得的計算精度一般要優于現有的數值解法，即計算效率更佳。同時，現有的數值方法在計算效率上要優于幾何方法。

5 應用

5.1 自動鋪帶

自動鋪帶技術是一種增量制造技術，適用于制造翼面、壁板等大尺寸、小曲率復合材料構件，具有加工成本低、效率高等特點。

在自動鋪帶工藝中，復合材料預浸帶在鋪帶頭的作用下沿著特定軌跡鋪貼到模具表面，如圖11所示。

圖11 自動鋪帶工藝原理圖Fig.11 Sketch of Automated Tape Placement

鋪帶軌跡需保證單條預浸帶不產生過大變形進而發生鋪放褶皺，同時需保證相鄰兩條預浸帶滿足一定的鋪放間隙要求。通常情況下，相鄰兩條預浸帶間不允許產生覆蓋，且鋪放間隙大小不能超過2.5mm。

現有的商用鋪帶軟件一般基于測地線規劃預浸帶的鋪放路徑以保證帶料變形最小[19]。

由于模具曲面高斯曲率的影響，相鄰兩條預浸帶間難以保持完全平行，導致鋪放間隙問題的產生。

一自由曲面模具，如圖12所示。其尺寸為（1000×5000）mm。L1、L2…L7是自由曲面模具0°鋪層中所鋪放的預浸帶。預浸帶型號為5228/T700，寬度為150mm。

圖12 自由曲面模具上相鄰帶料間的鋪放間隙分布Fig.12 Distribution of Gap on Freeform Mould Surface

為保證帶料變形最小，L1、L2…L7的初始鋪放路徑沿著測地線。由于該自由曲面不可展，相鄰測地線間無法保證完全平行，導致鋪放間隙問題的產生。如圖所示，相鄰兩條帶料L6與L7之間末端的鋪放間隙大小接近10mm，這將嚴重影響最終產品的質量。鋪放間隙的大小采用了文獻[20]中的方法進行計算。

接下來部分將采用提出的四階方法構造自由曲面上的等測地曲率曲線，以實現預浸帶在模具曲面上的微小轉向，消除帶料間過大的鋪放間隙。

采用等測地曲率曲線的目的在于將整條預浸帶需要的轉向均勻賦予到軌跡上的每一點，可以防止預浸帶在軌跡某一點處出現變形集中進而產生鋪放褶皺，從而保證鋪帶產品的最終質量。

顯然，調整后Li的鋪放路徑應該在調整前Li鋪放路徑的右側才能消除Li與Li-1之間過大的鋪放間隙，由第四節內容可知，Li的鋪放路徑上測地曲率符號為負。為獲得鋪放路徑上測地曲率的幅值，可采用折半搜索法，具體，如圖13所示。

圖13 折半查找流程圖Fig.13 Flowchart of the Bi-Section Search

在圖13 中，Upper_bound 為折半查找法的上邊界，Lower_bound 為折半查找法的下邊界，滿足|Upper_bound|＞|Lower_bound|，Δ1、Δ2—初始上邊界與下邊界，它們的數值，如表5所示。

表5 折半查找的初始上邊界與下邊界Tab.5 The Initial Upper Bound and Lower Bound for the Bisection Search

折半搜索法的停止條件為1.5mm ≤End_Gap≤2.5mm，其中End_Gap表示Li與Li-1之間末端的鋪放間隙大小。當折半搜索完成時，Li的鋪放路徑上測地曲率記為kg（Li）。

折半查找完成后測地曲率及鋪放間隙的大小，如表6所示。

表6 折半查找完成后測地曲率及鋪放間隙的大小Tab.6 Value of Geodesic Curvature and Gap After the Search Procedure is Finished

調整鋪放路徑前、后相鄰預浸帶間鋪放間隙的分布情況，如圖14所示。

圖14 調整鋪放路徑前、后鋪放間隙的分布情況Fig.14 The Distribution of Gap on the Mould Surface

其中一條曲線表示優化前鋪放間隙的分布，另一條曲線表示優化后鋪放間隙的分布。

由該圖可看出，采用本方法所構造出的等測地曲率曲線成功實現了預浸帶在模具曲面上的轉向，消除了帶料間過大的鋪放間隙。

優化前后模具曲面上的鋪放路徑，如圖15所示。

圖15 優化前后自由曲面模具上的鋪放路徑Fig.15 The Lay-Up Trajectories on the Freeform Surface

其中，一條曲線表示原始的鋪放路徑；另一條曲線表示調整后的鋪放路徑。

采用第三節中的四階方法進行軌跡規劃時鋪放路徑L1、L2…L7的計算精度以及7條曲線總的計算時間消耗（未包含折半搜索的時間消耗），如圖16所示。

圖16 本方法與文獻[16]中方法的對比Fig.16 Comparison of Our Method with Reference[16]

由該圖可看出，與現有的數值方法相比，本方法在相當的計算時間內可以得到更好的計算精度，即計算效率更高。

具體的說，當計算步長為10與1時，本方法平均計算誤差為2.3E-7 與3.47E-11，而文獻[16]中的方法平均計算誤差為7.3E-7與7.92E-11。

5.2 纖維纏繞

纖維纏繞工藝同樣是一種復合材料自動化加工工藝，一般應用于壓力管道、壓力容器等回轉體的制造。

近年來，該工藝也逐漸被應用于制造葉片、飛機S型進氣道、三通管等非回轉體異型面[21]。

本方法可應用于纖維纏繞工藝中非測地線纏繞軌跡的設計。非測地線纏繞對測地曲率的要求為：纏繞軌跡上各點的測地曲率與軌跡在該點處法曲率的比值不能超過纖維與模具之間的最大靜摩擦系數[22]。

現有的研究在處理異型面纏繞時通?；诰W格化方法，即先將原始曲面離散為由許多微小面片組成的網格曲面，然后在網格曲面上設計纏繞軌跡。

網格化方法的計算精度與離散網格的精細程度直接相關，當網格劃分很精細時，網格曲面的生成需要消耗較多的計算時間。另外，該方法也難以保證設計出的軌跡完全位于原始曲面之上。

本方法不需要生成離散曲面，只要目標曲面能夠通過參數方程進行表達，就可以設計出高精度的非測地線纏繞軌跡。

采用本方法設計非測地線纏繞軌跡時，需對F算子中公式進行調整，此時單位切向量t0在P0點相對弧長的變化率需通過下式計算：

式中：λ—纖維纏繞的滑移系數。

圖17顯示了一個基于自由曲面表達的異型管道，以及采用第三節中四階方法所設計出的不同滑移系數的非測地線纏繞軌跡，這些軌跡具有很高的計算精度且完全位于目標曲面之上。

圖17 異型管道上不同滑移系數的非測地線纏繞軌跡Fig.17 Non-Geodesic Winding Paths on a Non-Revolution Surface

6 結束語

提出了一種參數曲面上具有規定測地曲率曲線的幾何方法，當給出目標曲面上的起始點、方向以及測地曲率的分布時，可追蹤出符合設計要求的空間曲線。

與現有數值方法相比，本方法無需建立以及求解反映具有規定測地曲率曲線的復雜數學模型，且對于不同的參數曲面具有統一的求解格式；與現有的幾何方法相比，本方法利用多點平均可以得到很高的計算精度。

實驗結果表明，本方法在計算效率上要顯著高于現有的幾何法，且一般情況下要優于數值解法。本方法可應用于自動鋪帶軌跡規劃以及纖維纏繞工藝中非測地線纏繞路徑的設計。

為便于比較各方法的計算效率，文中未考慮在計算過程中動態調整計算步長。

如何合理利用每一步計算所得到的幾何信息設計變步長方法將是下一步的研究重點。