X射線聚焦望遠鏡鏡片裝配一體化建模和仿真

郝朝洋，邱樹涵，潘文婷，李 芬，顏 遠，關曉飛

（1.同濟大學數學科學學院，上海200092；2.同濟大學交通運輸工程學院，上海200092；3.同濟大學經濟與管理學院，上海200092）

近20年來，隨著高能天文觀測的蓬勃發展，國際上先后發射了多顆X射線天文衛星。這些衛星多采用嵌套式Wolter-I型［1］或圓錐近似的類Wolter-I型［2］望遠鏡結構，實現了高角分辨率和大集光面積的觀測，并取得了豐碩的成果。而中國在X射線天文領域的起步較晚，時間投入和資金投資相對有限，成果相對較少，目前為止僅自主研制了一顆X射線天文衛星，并在硬X射線調制望遠鏡（HXMT）［3］領域初步獲得一些成果。此外，中國計劃于2025年發射X射線時變與偏振探測衛星（XTP衛星）［4］，采用類似NuSTAR望遠鏡的圓錐嵌套Wolter-Ⅰ型結構［5］，并且用圓柱玻璃鏡片進行壓制，設計的角分辨率在1'左右，這將會是世界上首個搭載大面積X射線偏振望遠鏡的天文觀測臺，擁有超5 000 cm2的總探測面積。XTP衛星的發射將極大有利地推動極端條件下基本物理規律的研究進展。但是現今對光學元件的加工精度以及裝配檢測手段達不到XTP衛星的X射線聚焦能力的要求，導致X射線成像和探測技術發展困難重重。

圓錐嵌套Wolter-I型結構既相對Wolter-I型結構有著性能的提升，又沒有嵌套Wolter-I型結構那樣高的制作成本和難度，目前被廣泛采用。考慮在裝配初期進行柱面玻璃熱彎效果的檢查和控制，可以利用目前的計算機全息檢測技術［6］，來減少每層鏡片的幾何誤差。而影響裝配質量的另一個關鍵因素是裝配過程產生的偏差，對此則很難找到合適途徑在裝配過程中進行實時的檢測和控制，并且實驗時間長、材料成本高，人工、設施等成本高，故計算機仿真模擬是一個很好的解決手段。

目前利用有限元軟件模擬望遠鏡鏡片裝配過程研究的建模方式主要為手動建模，所建模型較為簡單，分析方式不夠精細，不能快速、自動、精確地得到裝配過程中所需的關鍵值以用于指導實際生產和制作。本文針對X射線聚焦望遠鏡鏡片裝配過程，采用有限元軟件Abaqus結合Python語言進行二次開發一體化地建立貼近實際的裝配模型，研究芯軸前后端半徑、石墨條間隔角度、鏡片厚度、鏡片半徑和裝配載荷對結構面形偏差的影響，節省建模及數據處理時間，方便進行大量模擬試驗并得到系列對比結果，以此指導實驗室的裝配實際流程。

1 鏡片裝配過程分析

圓錐嵌套Wolter-I型結構性能較好，又沒有嵌套Wolter-I型結構那樣高的制作成本和難度，是較為廣泛采用的結構。它的光路結構如圖1所示。圖中，Lp和Lh分別為主、次軸的軸向長度，g為主次軸間的間隙，Z為光軸方向，N為第N個主鏡入射X射線方向，d為相鄰兩個主鏡光路之間的距離，θN為主鏡與軸線之間的傾斜角，RN、rN分別為次鏡在中心和端點到光軸的距離，其大小與主鏡傾斜角相配合，使得射線經過每一層主次鏡都聚焦到O點，f為焦距。

若實際主次鏡面與理想錐面在入射位置和反射角度上有出入，則將出現成像誤差。成像誤差的具體機理如圖2所示。圖中，Δd為主鏡面在徑向上的面形誤差，θ1N、θ2N分別為主、次鏡面與軸線之間的傾斜角。令一束與軸平行的X射線射入，ΔL1為從主鏡發生一次反射之后的軸向誤差，ΔL2為再從次鏡發生一次反射之后的軸向誤差。本文選取的最大面型誤差范圍為3μm［7］。

該結構主、次鏡都疊合有多層鏡片，每層6片鏡片，每片角度約有60°，每個玻璃鏡片都撐有5根石墨條。裝配時把熱彎而成的柱面玻璃放在涂有超薄環氧樹脂膠的石墨條上，垂直于玻璃表面對橡膠壓條施加位移進行加壓，在石墨條與玻璃充分接觸之后，等待環氧樹脂膠凝固，最后按一定次序卸載壓條。表1為鏡片裝配時其中一組常用的幾何參數。

把柱面的玻璃壓制到圓錐面的安裝面上，玻璃會變形產生偏差。會影響到玻璃鏡片安裝的面形偏差因素有：玻璃自身熱彎質量、玻璃半徑（假設芯軸尺寸已固定）和裝配工藝等。本文默認玻璃質量達標，僅考慮玻璃半徑和裝配工藝等的影響。由于要進行大量的實驗對比分析，通過有限元軟件多次手工建模的成本很高，本文將通過有限元軟件Abaqus結合Python語言進行二次開發，建立一體化建模流程，在此基礎上進行仿真實驗，并嘗試為裝配優化工藝提供技術支持和建議。

2 一體化建模和仿真流程

2.1 模型設置

不同于賀鵬飛等［7］的相關研究，本文采用三維模型進行仿真和分析，考慮到芯軸剛度較高，在建模中不考慮其形變，并將石墨條與芯軸的接觸面簡化為固定約束。考慮到6個鏡片在結構上具有對稱性，但單個鏡片在裝配過程中并不具有對稱性，本文將模型簡化為研究一個鏡片。如圖3所示，1號壓條為中間壓條，2、3號壓條為二級壓條，4、5號壓條為三級壓條，壓條對應的石墨命名順序也相同，這里壓條加載順序為1—2—3—4—5。

圖3 三維黏接分析有限元模型Fig.3 3-D finite element model

建模和仿真采用Abaqus有限元分析軟件，為避免剪切閉鎖問題使用C3D8R單元，也可使用C3D8或C3D8I單元。局部有限元網格演示見圖3，各材料的性能參數見表2。

表2 主要材料參數Tab.2 Parameters of materials

2.2 荷載加載和邊界條件

本文研究不同半徑鏡面的加載黏合過程，其與芯軸前后端的半徑有以下3種關系：鏡面半徑大于芯軸前后端半徑、鏡面半徑小于芯軸前后端半徑和鏡面半徑介于芯軸前后端半徑之間。這3種情況下力荷載的作用有所不同，第一種情況下中間壓條1上的力荷載是為了使鏡片接觸石墨1（自然狀態下分離），之后在2～5號壓條上施加的力荷載與壓條1的作用相同，這時壓條1所需的力荷載會略小；對于第二種情況，壓條1開始不需要荷載（自然狀態下鏡片接觸石墨1），隨著2～5壓條上荷載的施加，壓條1需要一定荷載來防止鏡片與石墨分離；情況三則是前端為第一種情況，后端為第二種情況。對于不同參數下所需荷載的最優數值，以前的研究均采用的是反復試算的方法，這無疑會大大增加工作量，并且實驗結果的精確性不夠高。

本文將建立自適應的位移荷載求反力法減少計算量并提高計算精度。位移荷載求反力法的目的是使鏡片與石墨條接觸，使用位移荷載配合自由度約束可以取代力荷載。首先在建模過程中旋轉鏡面與壓條，使得鏡面中線、1號壓條與1號石墨條平行，再將2～5號壓條繞芯軸旋轉到1號壓條兩側，最后抬升鏡面與壓條使它們與石墨條分離；接著依次為1～5號壓條施加位移荷載，位移大小用幾何方法計算，并在每步對鏡面、壓條和石墨條施加相應的自由度約束。同時設置場輸出和歷程輸出時對壓條選擇反作用力（RF），這樣仿真結束后可以逐幀提取反作用力，其結果即為所需荷載，并可提取節點坐標計算面形誤差。位移荷載求反力法有其局限性，仿真實驗中發現多數情況下（實際實驗也證明）并非在剛好接觸的狀態下其面形偏差最小，而是在接觸的基礎上增加一點荷載（合力或均布壓強）可以使得面形偏差更小，于是對算法進一步改進建立自適應的位移荷載求反力法。

自適應的位移荷載求反力法的算法分為兩步：第一步為搜索階段——首先新設4個荷載參數，分別為1號壓條荷載、2、3號壓條荷載，4、5號壓條荷載（利用對稱性）和步長，接著對荷載參數在一定范圍內進行遍歷（此時選取的搜索步長是關鍵），按照之前的幾何參數進行幾何建模，施加位移，接觸后施加荷載，并取消位移約束，如果滿足接觸且計算得到的面形誤差小于位移荷載法計算得到的面形誤差，則記錄當前參數；第二步為優化階段——在搜索階段已經獲得了較優參數的大致值，接著再縮小步長，對參數附近的值進行遍歷，直到滿足需要。

以上面形誤差的計算根據鏡面幾何參數計算標準錐面方程，提取有限單元節點的三維坐標（x，y，z），根據y和錐面方程計算標準半徑R，再用計算徑向誤差（見圖1和圖2）。本文取所有節點最大誤差的絕對值作為整個鏡片的面形誤差。

2.3 Abaqus-GUI一體化開發

基于模型設置和載荷加載方法利用Python對Abaqus進行二次開發。Abaqus腳本接口可以基于Python語言的定制開發。這個接口與圖形界面無關，可與內核直接通信，可以執行重復性計算、創建模型、做優化分析、訪問輸出數據庫、創建Abaqus的插件等。使用Python對Abaqus進行二次開發的好處有：語言簡潔，易于理解，簡化手動繁瑣的操作；自動化，如進行后處理可以自己編寫專門的模塊；能夠實現參數化分析；可編寫獨立模塊，具有獨立性和可移植性；具有優秀的異常處理機制［8］。

Abaqus腳本與Abaqus的計算機輔助計算（CAE）模塊的通信關系如圖4所示。這里Python解釋器有3種輸入接口分別為GUI，命令行接口（CLI）和腳本，并保存為rpy文件。經過有限元分析計算以后結果保存并輸出數據庫中（ODB文件），其由模型數據和結果數據組成。輸出數據庫ODB文件的層次結構見圖5。

圖4 Abaqus腳本與Abaqus/CAE的關系Fig.4 Abaqus script versus Abaqus/CAE

圖5 輸出數據庫ODB的層級結構Fig.5 Hierarchy of the output database ODB

自適應建模分析算法利用Python編寫腳本實現，仿真生成ODB輸出數據庫文件后，則用Python腳本對仿真過程和結果中的關鍵數據進行提取，然后用數據分析工具（如Python、Matlab）進行分析。此外，對于傳統的手動建模過程中需要反復切換模塊、點擊按鈕、輸入數據等，而且難以撤銷，若使用純代碼進行建模，則每次都要修改部分代碼，并且用戶體驗不佳，操作比較復雜。本文同時設計了GUI界面進行仿真分析，有效地提高前后處理效率，節省建模及數據處理時間，利于用戶操作且還可以規避一些人為的錯誤。

界面的設計如圖6所示，有常規計算和確定最佳荷載組合計算兩種窗口。在常規計算窗口中，有幾何參數欄（前后端半徑、鏡筒高、石墨條厚度與旋轉角度、鏡面半徑與厚度、壓條厚度等）、材料參數欄、仿真參數欄（石墨條、鏡面、壓條網格大小）以及位移荷載求反力法和自定義荷載兩種模式，勾選自定義荷載后可輸入想要的壓條荷載；在確定最佳荷載組合計算窗口中，有幾何參數欄、仿真參數欄以及自定義搜索區間和在求反力法結果附近搜索兩種模式，勾選求反力法結果附近搜索后將先進行求反力法計算再搜索包含改計算結果的區間，此外加入參數輸入約束如各數值為正數、前端半徑需大于后端半徑等。計算的過程輸出將放在Abaqus安裝路徑下的Temp文件夾中。

圖6 GUI界面Fig.6 GUI interface

2.4 收斂性及精度的影響因素

在仿真過程中，有限元分析比較困難的是計算達到收斂，大多數情況下不收斂的原因難以找到，需要進行大量的試算。不收斂原因主要有：①網格劃分不合理。如網格劃分過粗或者接觸時主從面的網格粗細沒有掌握好，有時候主面和從面的網格大小為倍數關系會造成穿透現象而不準確且無法收斂。點對面離散時，若從面網格比主面網格細，則不會發生穿透，從面和主面都發生了正常的變形；若從面網格比主面網格粗，則會出現穿透。面對面離散時，情況類似。②針對接觸分析問題。可以設置微小的過盈量，以保證在分析一開始就已經建立起接觸關系。另外不能在接觸面上使用C3D20、C3D20和C3D10等單元類型，避免過約束。此外，對本問題使用減縮積分單元收斂效果會比完全積分單元和非協調單元要好。③分析步設置不合理，兩個復雜步驟應該放在不同的分析步里，分析步的初始步長、最小步長、步長縮減系數、最大嘗試次數都關系到能否收斂。在接觸分析步中，位移主要有兩個階段，一個是鏡面的兩邊還沒有接觸到最外側的石墨條，這時的運算很簡單，分析步長也很大，但一旦鏡面的兩邊接觸到了最外側的石墨條，也就是進入了第二階段，這時候運算開始變得復雜，分析步步長將急劇減小，若減小的嘗試次數超過5，則Abaqus會以為算不出來而報錯，這時候需要在分析步編輯器和通用求解控制編輯器中調整初始步長、最小步長和最大嘗試次數。④邊界條件約束不當，缺乏約束會產生數值奇異或零主元。⑤使用不恰當的位移荷載，由于精度問題，位移荷載可能會造成較大的局部應力，產生不合理的形變，使得準確度較低、收斂很困難。⑥使用不恰當的力荷載，如靜力分析模塊中在沒有其他約束的時候對物體施加純力，會導致分析步長越來越小直至不能收斂。

針對計算的精度，考慮到Abaqus軟件默認精度及計算機存儲精度的限制，可以利用Abaqus仿真計算的無單位性，對仿真參數進行放大或縮小以提高計算結果的有效數字，使得仿真計算能最大程度保留有效數字，提高計算精度。此外，還需要對X射線聚焦望遠鏡鏡面壓彎有關材料的彈性形變特性與塑性形變特性、鏡面與環氧樹脂、石墨間隔條的接觸摩擦進行細節上的處理，另外也可使用Fortran語言用戶子程序的開發，以上改進會使得數值仿真結果更為準確。

3 仿真結果分析

根據上述一體化建模和仿真流程，本文將依次研究前后端半徑、鏡片半徑、石墨條間隔角度等對裝配結果的影響，并對數據進行進一步的分析。下文中所有演示圖默認為前端在下，1號石墨條指向1號壓條為x軸正向，鏡片軸線向上為y軸正向，4號石墨條指向5號石墨條為z軸正向（圖3）。

3.1 前后端半徑、鏡片半徑對裝配結果的影響

根據鏡片半徑與芯軸前后端半徑的3種大小關系，本文用多組數據進行實驗，得出這3種關系應力分布、面型殘余的分布等基本特征。這里芯軸參數見表1和表2，鏡片內表面半徑見表3，鏡面厚度為0.2 mm，鏡筒高為100.0 mm，石墨條之間間隔角度為14.1°，石墨條與壓條厚度為1.0 mm，石墨條、鏡片、壓條網格大小分別為2、3、6 mm，計算結果見表3（僅展示其中3組）。

不同尺寸鏡面應力分布、面形誤差分布見圖7到圖9。總體來看，壓條附近區域面形誤差都較小，這是因為壓條本身與石墨條的貼近對鏡面的約束強度較大，而其他區域則缺乏直接約束；而壓條附近區域鏡面應力都比較大，但計算而得的最大應力均低于玻璃的強度極限。對于前兩種情況的鏡片，前端應力普遍大于后端，而第三種情況前端與后端應力分布較為均勻。從應力的最大值看，最大應力與最大誤差隨鏡片內表面半徑增加而減少。第一種情況下徑向實際值往往偏大，而第二種情況則偏小，第三種情況則較為對稱。

表3 基于三種鏡片內表面半徑的計算結果Tab.3 Computational results based on three inner surface radius of mirror

圖7 鏡面內表面半徑為83.5 mm時的應力和面形誤差分布Fig.7 Distribution of stress and surface error when inner surface radius of mirror is 83.5 mm

注意到雖然在最大誤差上三維模型與二維模型類似，但面形誤差分布上兩者還是有所不同。尤其是當鏡片半徑介于前后端半徑之間時，其分布情況較為特殊，且綜合來看鏡片半徑略大于前后端半徑的時候，最大應力與最大誤差以及誤差分布的情況都較好，而過大的鏡面半徑也會造成較大誤差。

從所需荷載的角度看，可以觀察到當鏡片內表面半徑小于前后端半徑時，壓條1最終所需荷載幾乎為0；介于兩者之間時，所需的荷載較大（要把其中一端壓緊）；大于兩者時，需要一點力來抵消其他4根壓條的影響，與預想的結果相符。

圖8 鏡面內表面半徑為84.5 mm時的應力和面形誤差分布Fig.8 Distribution of stress and surface error when inner surface radius of mirror is 84.5 mm

圖9 鏡面內表面半徑為85.5 mm時的應力和面形誤差分布Fig.9 Distribution of stress and surface error when inner surface radius of mirror is 85.5 mm

3.2 鏡片厚度對裝配結果的影響

當前的鏡片厚度主要有0.2、0.3、0.4 mm 3種，由上一小節知，鏡片半徑略大于前后端半徑時效果較好，故本小節所使用的鏡片半徑為85.5 mm，其余參數相同，以此分別分析 0.20、0.25、0.30、0.35、0.40、0.45、0.50 mm這7種規格鏡片厚度對裝配結果的影響，計算結果如表4所示。圖10展示其中一種情況下整體應力和面形誤差分布。

表4 7種鏡片厚度下計算結果Tab.4 Computational results of 7 different thicknesses of mirror

由圖11可知，鏡片厚度改變，應力與面形誤差分布規律相同。可以發現，隨著鏡面厚度的增加，最大應力增加（但仍遠低于玻璃的強度極限），所需荷載無規律變化，而最大面形誤差減少，因為根據材料力學，厚度大的板殼結構不容易發生較大的形變；但當鏡面厚度增加到一定程度后，應力增加較快，使得鏡面出現更大的形變，最大面形誤差會突增到比原來還大。經過多組數據分析后鏡片厚度在0.3～0.4 mm之間為佳，厚度過大或過小均會使面形誤差增加。

3.3 石墨條間隔角度對裝配結果的影響

由上一小節知鏡面厚度為0.4 mm時面型最大誤差小，故本小節所使用的鏡片厚度為0.4 mm，其余參數與上一小節相同，以此分析石墨條間隔角度對裝配結果的影響。表5給出了石墨條間隔角度不同時的計算結果。石墨條間隔角度改變，應力與面形誤差分布規律相同，此處不再贅述。可以發現，隨著石墨條間隔角度的增加，最大應力先減后增，而最大誤差的改變則沒有明顯規律（總體上先減后增），故石墨條最佳間隔角度只能通過自適應算法求得。需要注意的是，石墨條間隔角度的影響要綜合考慮每一層鏡片的情況，要找到對每一層鏡片來說最大面形誤差最小的石墨條最佳間隔角度，而不能只針對某一層最小。這些都體現出了一體化建模和仿真的優勢。

圖10 鏡片厚度為0.35 mm時的應力和面形誤差分布Fig.10 Distribution of stress and surface error when thickness of mirror is 0.35 mm

圖11 鏡片厚度與最大應力、最大面形誤差的近似關系Fig.11Mirror thickness versus maximum stress and surface error

3.4 最佳荷載組合研究

本文利用自適應算法在一定區間內對各鏡片半徑下的最佳荷載組合進行探究，裝配相關參數與鏡片厚度為0.3 mm的算例相同，得到結果見表6。從計算結果可看出，最佳荷載組合與施加位移求反力的結果是相近的，且由于所需荷載的數值較小，算法的遍歷搜索步驟可能要執行很長時間，故在施加位移求反力的結果附近搜索效率更高。

表5 不同石墨條間隔角度的計算結果Tab.5 Computational results of different angles of graphite strips

4.5 面形誤差的進一步表示

以前的研究較為重視徑向誤差對像差的影響［7］。由光路圖可知，影響鏡面成像的質量不只是徑向誤差Δd，鏡面傾斜角度的誤差Δθ對鏡面成像的質量也有一定影響，故面形誤差實際上應當用如下關系式來表示（此處以主鏡的誤差為例）：

對于Δθ的求算，本文考慮用曲面擬合仿真后鏡面得到的散點，裝配相關參數與鏡片厚度為0.4 mm的算例相同。具體算法如下：找到所求節點周圍的點（邊界點處理稍微特殊），對這些點進行二次曲面擬合，獲得曲面方程，求曲面在此處的法向量，并與理想錐面的法向量作比較。所得的傾斜角誤差分布圖如圖12所示。

表6 半徑82～86 mm鏡片最優荷載組合及最大應力Tab.6 Optimal loads and maximum stress of mirrors with a radius from 82 to 86 mm

從傾斜角誤差分布的情況來看，越往中間鏡面的角度誤差就越大，其計算精度與曲面擬合的程度有關。

圖12 傾斜角誤差分布圖Fig.12 Error of tilt angle

4 結論與展望

本文針對圓錐嵌套式Wolter-I型X射線聚焦望遠鏡鏡片裝配過程利用Abaqus建立了三維一體化建模和仿真算法流程，同時使用Python對Abaqus進行二次開發，提出了自適應建模計算方法，本方法可以有效提取各參數下的應力分布、面形偏差分布等結果。利用此算法流程對3種不同半徑鏡片進行加載，分析了前后端半徑、鏡片厚度、鏡片曲率和裝配載荷對結構面形偏差的影響，得到其面形誤差與二維模型有著較大的不同，在鏡片厚度偏小和偏大都會使得面形誤差偏大，理想的厚度應在0.3～0.4 mm之間，石墨條間隔角度增加往往會減小面形誤差，但要綜合考慮各層鏡面的情況等結論，并利用自適應算法自動地為不同半徑鏡片計算最佳裝配載荷。在面形偏差的評估方面，本文提出了需要綜合考慮面形徑向偏差與面形傾斜角度變化的思想，并用曲面擬合的方法來估算鏡面在某點處的法線方程，以此來對光路進行進一步的分析。此外，本文相應二次開發的算法和軟件能夠一體化地建模和仿真，可用于多層鏡片計算，節約可觀的實驗成本和人工，具有良好的可拓展性和廣泛的適用性。最后得到的最優鏡片半徑、最佳裝配載荷和考慮多種面型偏差影響因素等研究結果，為提高裝配工藝減少面形偏差提供了理論分析方法和裝配指導方向。