等幾何拓撲優化非齊次邊界條件施加技術研究

摘要：非均勻有理B樣條（NURBS）基函數缺乏插值性，導致等幾何拓撲優化方法無法直接在控制頂點處施加非齊次邊界條件。為此，介紹了面向等幾何分析的兩類非齊次邊界條件處理方法，通過強施加方法和弱施加方法分別對邊界條件施加的精度進行提升。數值算例結果驗證了前述方法在等幾何分析非齊次邊界條件施加過程的可行性。將前述方法嵌入至等幾何拓撲優化流程，創新性地提出懲罰因子自適應的非齊次邊界條件處理方法，提高了等幾何拓撲優化求解的精度，進一步論證了文中方法的有效性。

關鍵詞：等幾何拓撲優化；非齊次邊界條件；邊界配點法；自適應罰函數法

中圖分類號：TH122

Research on Non-homogeneous Boundary Condition Imposing Techniques for Isogeometric Topology Optimization

WANG Shuting1 XIE Qingtian1 YANG Aodi1 LI Xiaobing2 XIONG Tifan1 XIE Xianda^2*

1.School of Mechanical Science and Engineering，Huazhong University of Science and Technology，Wuhan，430074

2.School of Advanced Manufacturing，Nanchang University，Nanchang，330031

Abstract： The lack of interpolation property of non-uniform rational B-splines（NURBS） basis functions led to the failure of imposing non-homogeneous boundary conditions directly at control points by isogeometric topology optimization method. Therefore， two types of non-homogeneous boundary condition imposing methods were introduced to isogeometric analysis herein， which improved the accuracy of imposing boundary conditions through two perspectives， namely strong imposition method and weak imposition method. Then， the numerical results verify the feasibility of the above methods in applying non-homogeneous boundary condition for isogeometric analysis. Finally， the above methods were embedded into isogeometric topology optimization， and an innovative non-homogeneous boundary condition imposing method was proposed with adaptive penalty factor， which improves the solution accuracy of isogeometric topology optimization， and further demonstrates the effectiveness of the proposed method.

Key words： isogeometric topology optimization; non-homogeneous boundary condition; boundary collocation method; adaptive penalty function method

0 引言

等幾何拓撲優化^［1-3］方法以等幾何分析代替傳統拓撲優化方法中的有限元分析，由于非均勻有理B樣條（NURBS）基函數的非插值性質，導致描述幾何的控制點通常不在幾何體上，這種情況下，控制點自由度不能直接控制幾何體的形狀，而需要通過調整權重來實現。在處理非齊次本質邊界條件時，等幾何分析無法直接準確施加邊界條件^［4］，影響了分析計算精度。同時，狄利克雷邊界條件對非齊次較為敏感，使得在等幾何拓撲優化過程中施加邊界條件的精度受到較大限制。HUGHES等^［5］將邊界上控制點處的場變量值直接施加到對應的控制點上，但這種方式在一般情況下只適用于常數邊界等齊次的邊界條件，對于非齊次邊界條件將會導致精度下降以及收斂速度變緩。

有研究表明，NURBS基函數和無網格法中采用的移動最小二乘基函數都具有高階連續性且都有非插值性^［6-7］，它們的基函數可以通過控制點位置進行修改，以產生形狀變化或實現幾何約束；同時，它們的基函數也具有高度的連續性，可以較為方便地計算導數。

國內外學者從無網格法^［8-10］中借鑒相關方法，應用至等幾何分析的非齊次邊界條件施加。參照無網格法，等幾何分析將邊界條件修正方法分為兩類^［11］：第一類是對解空間進行修正的“強施加”處理方法；第二類是對變分形式進行修正的“弱施加”處理方法。

關于第一類方法，WANG等^［12］采用配點法思想施加等幾何分析的邊界條件，提出了一種改進的強制邊界條件處理方法。隨后，軒軍廠^［13］將上述方法和應變光滑子域積分結合，該算法將單元分割為若干子域， 通過子域的應變平均定義應變光滑，進一步提高了分析精度和計算效率。陳濤等^［14］證明了配點選取方法對計算精度的影響。

第二類方法主要包含罰函數法、Nitsche’s法。王東東等^［15］將罰函數法用于等幾何分析的邊界條件處理，與直接施加法對比，得到了更好的收斂效果和精度，并將其應用在三維薄梁板結構中取得了較好的效果^［16］；HOANG等^［17］將罰函數懲罰理念應用至流體力學問題；LEI等^［18］在C0/G1多面片的應力約束條件下研究了基于罰函數法的等幾何分析技術，在Kirchhoff-Love殼等結構中具有較好的有效性。EMBAR等^［19］討論了基于Ｂ樣條的有限元分析中施加狄利克雷邊界條件的問題，并提出了將Nitsche’s應用于等幾何分析；CHEN等^［20］在采用 Nitsche’s變分邊界的基礎上推導出穩定系數的計算公式；NGUYEN等^［21］采用 Nitsche’s法處理二維和三維多片結構之間的拼接問題；HU等^［22］在彈性力學問題中引入斜對稱Nitsche’s公式處理邊界問題；DU等^［23］將Nitsche’s應用至橡膠、生物組織等多材料的復雜幾何模型。在接觸力學領域，胡清元^［24］將Nitsche’s方法用于等幾何分析處理位移、轉動和耦合邊界條件，得到了較高精度的分析結果。此外，其他學者也貢獻了邊界處理的解決方案^［25-28］。

目前，等幾何分析中對非齊次邊界條件的處理基本圍繞借鑒無網格方法展開研究，并取得了一定進展，但在等幾何拓撲優化方法中對邊界條件施加的相關研究較少。實際上拓撲優化過程中，每個迭代步數都會影響單元材料屬性，而現有非齊次邊界條件等幾何拓撲優化方法未考慮單元材料屬性變化對邊界條件施加技術精度的影響，這也會導致在非齊次邊界條件下等幾何拓撲優化精度降低。

本文針對上述問題，基于非齊次邊界條件的強施加與弱施加方法，研究基于等幾何分析的拓撲優化在非齊次邊界條件下的精確等效施加方法，主要基于格雷維爾標記配點法的強施加方法，以及基于罰函數法的弱施加方法。結合優化過程中目標函數及剛度矩陣特征值的變化，對懲罰因子的取值進行優化設計，研究自適應罰函數法的等幾何拓撲優化方法。

1 理論基礎

本節首先對邊界條件的分類和定義進行闡述，然后基于NURBS非插值特性對等幾何方法的具體影響及內在機制進行了分析。

1.1 非齊次邊界條件

邊界條件是給出具體物理現象在邊界上所處的物理情況。根據邊界條件數學表達方式的不同，在偏微分方程求解中將邊界條件分為狄利克雷邊界條件、紐曼邊界條件以及柯西邊界條件三類。

其中，在等幾何分析過程中，后兩種邊界條件一般直接構造在方程的等效積分形式中就能滿足條件，下式展示了三種邊界條件的形式：

u=h""""" on ΓD

n·σ=ton ΓN

αu+n·σ=ron ΓR（1）

式中：u、h分別為實際邊界位移和給定邊界位移;n為邊界外法線方向;σ為應力張量；t為邊界處牽引力；α為用于確定邊界類型的常數；r為外載荷向量；ΓD、ΓN和ΓR分別為狄利克雷邊界、紐曼邊界和柯西邊界。

當h=0時，被稱為齊次邊界條件；反之，為非齊次邊界條件。在非齊次條件下，等幾何拓撲優化過程求解結果不穩定，精度下降。由于非齊次邊界條件可能導致解的奇異性或多解性，在彈性力學中，這可能會導致應力和應變分布的不連續或奇異點的出現，從而影響結構的穩定性和安全性。

1.2 NURBS非插值特性

在有限元分析中，由于其近似函數具有插值特性：基函數在其控制點處取值為1而在其他控制點處取值為0，故相對容易處理邊界條件。但NURBS基函數由于是通過遞歸地應用Cox-de Boor遞推公式來計算的，這個過程中權重因子使得基函數在其控制點處取值不為1，導致其不具有插值特性，因此對邊界條件的處理較為困難。如圖1所示，節點向量為Ξ=（0，0，0，0.25，0.5，0.75，1，1，1）的B樣條基函數，可觀察到：除了包含函數端點位置的基函數（N1和N6）外其他都不具備插值特性。

將B樣條推廣至多變量范圍，設兩個參數方向都采用上述節點矢量，可以得到圖2所示的雙變量NURBS基函數。

對雙變量 NURBS基函數也只在參數域的四個拐角點處插值，如圖2a所示，這類控制點被稱為角控制點。圖2b和圖2c分別表示其他邊界和內部節點，被稱為邊界控制點和內部控制點。可以看到，在二維空間中NURBS基函數也是非插值的。同時，在邊界位置處內部控制點的形函數為0，因此邊界值的大小只由角控制點和邊界控制點控制。

在幾何特性范疇，NURBS樣條的凸包性質限制了實際樣條幾何只能包含在其凸包內部。圖3以二維1/4圓環控制網格為例，分別給出了等幾何與有限元的控制網格，圖中實心圓點為控制頂點（或節點），連接控制頂點的虛線構成控制網格。可以看到由于NURBS 樣條缺乏插值性，控制頂點可能在幾何區域的內部、邊界以及外部（其中，紅色、黃色和藍色控制點分別表示角控制點、邊界控制點和內部控制點），而有限元的二次拉格朗日等參單元的節點都落在幾何區域內部或者邊界上。因此，需要對控制頂點進行適當調整，以確保實際的幾何區域被包含在調整后的凸包內部。若直接將非齊次本質邊界條件施加在控制點上，則可能會破壞NURBS樣條的凸包性質，導致計算精度降低甚至收斂退化。

1.3 非齊次邊界條件施加方法分類

為了能夠提高非齊次邊界條件下等幾何拓撲優化的計算精度，提高收斂特性，本文對邊界條件進行處理，并將處理方法分為強施加和弱施加兩類方法，如圖4所示。

第一類是強施加方法，該方法直接對解空間進行強修正，要求邊界條件必須滿足，直接將其作為方程的一部分來求解。第二類是弱施加方法，通過將邊界條件加入問題的弱形式中，將其轉化為一個變分問題，在變分問題中邊界條件以積分形式出現。這種方法通常用于邊界條件較為復雜或難以處理的情況。

目前主流的強施加方法包括邊界配點法、位移約束方程法、最小二乘擴展法、等價函數構造法；弱施加法包含罰函數法、拉格朗日乘子法、修正變分原理法、Nitsche’s法。為了能夠更好地處理非齊次邊界條件，本文分別研究了強施加方法及弱施加方法對邊界條件進行修正，提高計算精度。

2 基于強弱施加法的非齊次邊界條件施加技術

本節首先介紹了非齊次邊界條件下的等幾何拓撲優化的數學模型，然后介紹了基于B樣條的隱式過濾器。提出兩類非齊次邊界條件處理方法，通過強施加方法和弱施加方法兩個角度來提高邊界條件施加的精度。

2.1 非齊次邊界條件下等幾何拓撲優化模型

在非齊次的紐曼狄利克雷混合邊界中，考慮一個彈性力學模型：

·σ+f=0" in Ω

u=hon ΓD

n·σ=ton ΓN（2）

式中：為梯度算子，Ω為計算域； f為體積力向量；ΓN∪ΓD=Γ，ΓN∩ΓD=。

得到靜載荷作用下二維彈性問題的代數方程：

KU=F（3）

K∈R^Ndof×Ndof

U∈R^Ndof" F∈R^Ndof

式中：K為總體剛度矩陣；U、F分別為整體結構控制點位移和外力矢量；Ndof為控制點總自由度數。

根據有限元分析邊界定義^［29-30］，可以將F可分解為內部力向量、紐曼邊界條件以及狄利克雷邊界條件：

F=∫ΩRfdΩ+∫ΓNRtdΓN+∫ΓDRhdΓD（4）

式中：R為NURBS基函數矩陣。

為便于表達，將控制點分為邊界控制點和內部控制點。圖5以二維平面中的四分之一圓環模型為例展示了邊界控制點（圖5中的藍色控制點）和內部控制點（圖5中的紅色控制點）。

為了更準確地對邊界條件進行控制，借鑒了無網格法^［31］的變換方法，將控制變量與它們相應的物理邊界值聯系起來，將邊界單元與內部單元進行區分。因此，式（3）可以轉化為

KIIKIBKIBKBBUIUB=FIFB（5）

式中：KII、KIB和KBB分別為內部控制點剛度子矩陣、內部控制點與邊界控制點耦合剛度子矩陣和外部控制點剛度子矩陣；UI、UB分別為內部位移向量和外部強施加位移向量；FI、FB分別為內部控制點力向量和外部控制點力向量。

本文拓撲優化的目標是最小化廣義柔度，廣義柔度參考總勢能（total potential energy， TPE）概念^［32-33］，計算公式為

Π（a1，a2）=12UTKU－FUI（6）

a1=（ρ00，ρ01，…，ρn1n2）

a2=（ω00，ω01，…，ωn1n2）

式中：a1、a2分別為全部控制點密度值和權重值。

由TPE可以推導出非齊次邊界條件下的廣義柔度：

C（a1，a2）=－2Π（a1，a2）=FUI－UTBFB（7）

以柔度最小化為目標，優化問題以約束非線性規劃問題的形式表示為

min C（a1，a2）Cref

s.t. KU=F

V（a1，a2）V0≤Volfrac

a1τ∈［ρmin，ρmax］，a2τ∈［ωmin，ωmax］

τ={1，2，…，ncp}（8）

式中：Cref為給定的結構柔度的標準值；V（a1，a2）為優化結構的實體材料體積；V0為初始設計域體積；Volfrac為實體材料體積分數上限；a1τ、a2τ分別為每個控制點的密度值和權重；ρmin、ρmax分別為密度值下限及上限（一般情況下，ρmin≥0，ρmax≤1）；ωmin、ωmax分別為權重下限和上限；ncp為控制變量總數。

2.2 基于配點法的強施加法及實現

本節首先對邊界配點法的基本原理和實現過程進行介紹。隨后對直接配點法進行改進，設計了更為精準的配點選取的準則以及配點方案。

2.2.1 配點法

考慮對狄利克雷邊界條件進行處理，將式（2）中ΓD表示為非齊次邊界條件與齊次邊界條件之和：

u=uhom+unon=∑ni=1Ri（a）ui（9）

式中：u為狄利克雷邊界；uhom、unon分別為齊次部分和非齊次部分；Ri（a）為NURBS基函數；a為其參數坐標；ui為對應的控制變量。

進一步，對式（9）進行排序重組，根據圖3a的控制點分類思想，將其重寫為

uhom=∑wi=1R^homi（a）u^homi=0

unon=∑mj=1R^noni（a）u^noni

uhom∪unon=u（10）

式中：R^homi（a）為內部基函數，由于其支撐域與邊界交集為，因此在邊界處為0；R^noni（a）為邊界基函數；u^homi、u^noni分別為內部和邊界控制變量。

因此，由式（9）可得到

u=∑mj=1R^noni（a）u^noni=R（a）bCb=TTCb（11）

式中：R（a）b、Cb分別為非齊次邊界條件的基函數和控制變量集合，T={R^noni（a）}，可以得到Cb=T^－Tub。

由上述可知，對非齊次邊界條件進行配點法強施加需要首先確定狄利克雷邊界中的非齊次部分unon，然后通過等效積分對載荷向量和剛度矩陣進行修正，最終轉化為線性方程組求解內部的控制變量。由于NUBRS是開節點曲線，基函數在節點矢量的端點處插值，只有一個基函數在端點處有非零值，因此可通過尋找支撐域與邊界條件非空交集的基函數確定邊界和內部基函數。

在邊界曲線上引入一組插值配點{bw}gw=1代入式（9）中，可得到配點方程組：

t（bw）=∑gi=1Ri（b）ci=R（b）pCp（12）

其中，線性方程組形式中的R（b）p矩陣由Ri（b）構成。

得到在非齊次邊界條件下的配點法精確等效施加邊界的整體流程，如圖6所示。在定義非齊次邊界條件后，識別邊界控制點并獲取其所在單元索引；在組裝整體剛度矩陣及載荷向量后判斷是否為邊界單元，并更新邊界單元剛度矩陣及載荷控制點。

2.2.2 配點法算法實現

在等幾何拓撲優化過程中，配點法主要通過對迭代過程中的剛度矩陣和載荷向量優化實現，將實現原理代入式（5），可得到優化后的矩陣形式：

D=

=KIIKIBT^-T

T^-1（KIB）TT^-1KBBT^-T

D=CIUB=I00TTCIUB，

=IB

（13）

KII=∫Ω（RI）TRIdΩ

KIB=∫Ω（RI）TRBdΩ

KBB=∫Γ（RB）TRBdΓ

I=FI=∫Ω（RI）TlσdΩ

B=∫Ω（RB）TsdΩ+∫ΓN（RB）Tt-dΩ（14）

式中：KII、KIB和KBB分別為內部剛度子矩陣、內部與邊界耦合的剛度子矩陣和邊界剛度子矩陣；I、B分別為內部載荷分向量和邊界載荷分向量。

通過配點法進行強施加的過程中還需要考慮如何合理地選擇配點，以避免配點密度過小而導致誤差增大。選擇配點的基本準則：為保證配點方程組正常求解，系數矩陣必須是非奇異的，且盡可能使條件數變小。根據勛伯格惠特尼（Schoenberg-Whitney）定理，在選取配點插值時應盡可能保證避免配點超出或接近基函數支撐邊界，以保證求解精度。本節研究了兩種選取方式，第一類是均勻等分法，將樣條函數的定義域均勻等分；第二類是格雷維爾標記法，采用節點平均值作為配點，對于給定一組節點矢量（ζi）^n+p+1i=1，配點（ζ′i）ni=1標記為連續p個節點均值，計算方法如下：

ζ′i=1p∑pj=1ζi+j（15）

2.3 基于自適應罰函數法的弱施加方法及實現

本節首先對罰函數法的基本原理和實現過程進行介紹，并通過等幾何分析算例進行數值驗證，證明其可行性和對精確度的提升。隨后，對罰函數法進行改進，修正了懲罰因子取值標準，使其能夠更加精確等效施加邊界條件。

2.3.1 罰函數法及算法實現

設在空間R中的約束集合上尋求一個實值函數的極值。那么對目標函數f（x）而言，需要增加定義一個罰函數g（x），得到新函數f′（x）：

f′（x）=f（x）+g（x）（16）

將罰函數引入彈性力學問題中，可以得到罰函數的修正泛函，變分后得到

δΦ=Φ+α∫Γbδ（u（η））T·（u（η）－u-b（η））dΓ（17）

式中：α為懲罰因子。

將其引入二維彈性力學問題，得到

（K+Kp）D=F+Fp（18）

其中，Kp和Fp分別是剛度和載荷懲罰修正項：

Kp=α∫Γb（N（η））T·N（η）dΓ

Fp=α∫Γb（N（η））T·u-b（η）dΓ（19）

罰函數的系數矩陣是對稱正定矩陣，其結構與原始剛度矩陣相同，并且沒有增加自由度數量。

2.3.2 懲罰因子自適應修正方法

在等幾何分析中，罰函數法的計算精度主要受懲罰因子影響，而懲罰因子的選取并無直接現有規律可循。懲罰因子從理論上講應該趨近于無窮大，可以更好地滿足非齊次邊界條件，以提高收斂速度。然而懲罰因子過大容易導致稀疏矩陣趨向病態，因此在不同的工程情況需取不同量級的懲罰因子。

此外，在基于具有懲罰的實體各向同性材料（SIMP）法的等幾何拓撲優化過程中，由于每一個迭代步數都需要更新全局剛度矩陣以及所對應的單元密度值，因此也需要懲罰因子能夠隨著密度的更新而不斷自適應調整。為了使剛度矩陣修正項與全局剛度矩陣實時匹配，本文對懲罰因子的自適應修正進行研究。

在等幾何分析中，懲罰因子一般為α∈［103E，107E］（E為彈性模量）。這種取值方法無法實現懲罰因子自適應，在等幾何拓撲優化中無法有效匹配迭代過程的變化，且無法與模型幾何尺寸構建關聯。本文提出懲罰因子自適應調節技術，在優化迭代開始前設定初始懲罰因子：

α=K^maxK^maxp（20）

式中：K^max、K^maxp分別為剛度矩陣K和Kp特征值的最大值。

當進入優化迭代過程后，設計懲罰因子自動遞增方法，在連續兩次迭代之間目標函數值減小較慢的情況下更新懲罰參數值。

αn+1=

αn+0.1K^maxK^maxp" 45|On|lt;|On+1|lt;|On|

αn其他（21）

式中：On、On+1為迭代過程中的目標函數值。

自適應罰函數法的實現流程如圖7所示。定義非齊次邊界條件后，識別邊界控制點及所在單元索引，并設置初始懲罰值，進入迭代優化過程后，判斷目標函數收斂情況，根據其變化情況修正懲罰因子值，并對非齊次邊界單元的剛度矩陣及載荷向量進行高斯線積分修正。

3 數值算例

本節基于強弱施加兩種非齊次邊界條件施加技術原理，首先通過等幾何分析方法進行算例驗證，證明該方法在等幾何方法的可行性，然后將其嵌入等幾何拓撲優化方法內，進一步驗證等幾何拓撲優化非齊次邊界條件施加技術。硬件測試環境為Intel（R） Xeon（R） Gold 5120 CPU @ 2.20GHz 2.19GHz內核和32GB RAM內存的工作站；軟件測試環境采用Windows10 64位操作系統和MATLAB_R2021a。

3.1 強施加法算例驗證

為驗證配點法的有效性，利用等幾何分析中二維1/4圓環區域求解作先行驗證。幾何區域采用二次NURBS樣條曲面定義，該二維圓環問題描述如圖8所示。該算例內徑a=0.3 m，外徑b=0.5 m，受內壓p=3×104 kN/m，彈性模量E=30 GPa，泊松比ν=0.25。施加非齊次邊界條件左邊界ux=I，下邊界uy=I（I為全1方陣）。

2D annulus

存在徑向和切向兩個方向的應力和位移解析值：

σr（r）=a2pb2－a2（1－b2r2）

σθ（r）=a2pb2－a2（1+b2r2）

ur（r）=a2prE（b2－a2）［1+b2r2（1+ν）－ν］

uθ（r）=0（22）

根據已知條件可得到應力最大理論值σmax=63.750 MPa，位移最大值umax=7.1250×10^－4 m。取4組不同單元數量（4×4單元、8×8單元、16×16單元及32×32單元）的結果誤差作對比，比較直接施加法、均勻配點法和格雷維爾標記法的計算結果。圖9展示了應變和應力誤差值。

兩類邊界配點法在應變和應力方面都表現出明顯優于直接施加在控制頂點上的方法。在應變方面，相較于直接施加法，均勻配點法的誤差減小了8.87%～69.68%，格雷維爾標記法的誤差減小了19.25%～80.16%；在應力方面，均勻配點法的誤差減小了91.53%～94.53%，格雷維爾標記法誤差減小了92.41%～98.18%。同時也應注意到，兩類配點方法的分析結果稍有差異，其中格雷維爾標記法效果更好。

進一步，通過對經典的拓撲優化二維1/4圓環算例進行優化計算，以驗證非齊次邊界條件下的配點法邊界施加方法的有效性。二維1/4圓環的初始設計域及約束情況見圖10，最大迭代步數設置為301步。

本文設置了3組對比參照組，如表1所示，在各類的設計網格分辨率（nelx×nely）下對傳統等幾何拓撲優化方法、格雷維爾標記法進行比較。將上述的每一組參數作為等幾何拓撲優化方法的輸入參數，同時，為了收斂結果保持穩定，將模型最下方邊界單元的網格密度固定。圖11以組別分類，分別展示了不同網格尺寸下的優化結果，圖12展示了兩種方法優化過程中的KKT（Karush-Kuhn-Tucker）條件殘差的二范數變化趨勢。

圖11中的優化結果表明，傳統方法與本文設計的格雷維爾標記法可以得到相似的優化結果：兩者在網格密度較小（100×100）的情況下都出現了棋盤格現象，隨著網格密度不斷增大，優化結果也更加理想，得到了更多的結構細節和更高的清晰度。圖12表明了格雷維爾標記法的有效性，該方法能夠在非齊次邊界條件下的等幾何拓撲優化中提高計算精度，減小設計誤差。相較于傳統方法，格雷維爾標記法將KKT條件殘差二范數降低了6.33%～8.50%（為避免優化初始時目標函數的劇烈振蕩對結果造成的誤差，本文選取第51步作為計算起始步）。

3.2 弱施加法算例驗證

為驗證罰函數法的有效性，利用3.1節中的二維1/4圓環等幾何分析算例做先行驗證。取4組單元數量的結果誤差作對比，比較直接施加法、格雷維爾標記法和自適應罰函數法的計算結果。圖13展示了應變和應力誤差值。圖中顯示，罰函數法在等幾何分析中的表現明顯優于直接施加法和格雷維爾標記配點法，其內在原因是：基于配點法的強施加方法對復雜邊界問題可能會導致數值不穩定或收斂困難，而基于罰函數法的弱施加方法可以更加準確地處理復雜邊界條件，具有良好的數值穩定性，詳情可參考文獻［4］；對邊界裁剪單元的懲罰因子進行自適應更新能提高求解結果的精度^［34］。由對比圖可以看出，在應變方面，自適應罰函數法與直接施加法相比，誤差減小了96.48%～99.90%，與格雷維爾標記法相比，誤差減小了95.65%～99.50%；在應力方面，自適應罰函數法能夠實現與精確值相同結果，消除了等幾何分析誤差。

進一步，使用二維圓環算例對自適應罰函數法在等幾何拓撲優化中的有效性進行驗證，初始條件和對比組參數設置同3.1節，對比了自適應罰函數法、格雷維爾標記法和傳統方法在不同設計網格分辨率情況下的優化結果、收斂過程、優化過程中KKT條件殘差的二范數。自適應罰函數法的優化結果如圖14所示，收斂過程KKT條件殘差范數變化如圖15所示。

自適應罰函數法可以得到與傳統方法及配點法相似的優化結構，并且展示了更加細致的結構細節，同時可以提高收斂速度，在100×100的網格尺寸下，自適應罰函數法在256步就實現了收斂，相較于其余兩種觸發最大迭代步數，更快達到收斂效果。圖15中的KKT條件殘差二范數也證明了自適應罰函數法的有效性，相較于傳統方法和格雷維爾標記法，將KKT條件殘差二范數分別減小了8.71%～9.34%和1.51%～3.03%。

4 結論

本文針對現有等幾何拓撲優化的邊界條件施加方法存在計算精度不高、收斂效果不理想的問題，設計了基于強施加方法和弱施加方法的等幾何拓撲優化非齊次邊界條件施加技術。基于格雷維爾標記法設計了配點法施加邊界條件，并基于罰函數法設計了自適應罰函數法，

利用優化過程中目標函數的更新實時修正懲罰因子，實現準確施加懲罰。通過等幾何分析數值算例對兩種方法的有效性進行驗證，通過非齊次邊界條件等幾何拓撲優化算例進一步驗證。相較于傳統方法，格雷維爾標記配點法可以減小KKT條件殘差二范數6.33%～8.50%，自適應罰函數法可減小8.71%～9.34%，同時可以提高收斂速度14.95%。

參考文獻：

［1］ 謝賢達. 基于等幾何分析的移動可變形組件拓撲優化方法及應用研究［D］. 武漢：華中科技大學， 2021.

XIE Xianda. Research on Topology Optimization Method and Application of Mobile Deformable Components Based on Isogeometric Analysis［D］.Wuhan：Huazhong University of Science and Technology， 2021.

［2］ 丁延冬， 羅年猛， 楊奧迪， 等. Bézier單元剛度映射下的高效多重網格等幾何拓撲優化方法［J］. 中國機械工程， 2022， 33（23）：2801-2810.

DING Yandong， LUO Nianmeng， YANG Aodi， et al. Efficient Multigrid Isogeometric Topology Optimization under Bézier Element Stiffness Mapping［J］. China Mechanical Engineering， 2022， 33（23）：2801-2810.

［3］ 楊雨豪， 鄭偉， 王英俊. 一種自由度縮減和收斂加速的高效等幾何拓撲優化方法［J］. 中國機械工程， 2022， 33（23）：2811-2821.

YANG Yuhao， ZHENG Wei， WANG Yingjun. An Efficient Isogeometric Topology Optimization Method Using DOF Reduction and Convergence Acceleration［J］. China Mechanical Engineering， 2022， 33（23）：2811-2821.

［4］ 陳濤. 等幾何分析方法的本質邊界條件處理研究［D］. 西安：西北工業大學， 2016.

CHEN Tao. Study on the Treatment of Essential Boundary Conditions of Equal Geometry Analysis Method［D］.Xi’an：Northwestern Polytechnical University， 2016.

［5］ HUGHES T J R， COTTRELL J A， BAZILEVS Y. Isogeometric Analysis：CAD， Finite Elements， NURBS， Exact Geometry and Mesh Refinement［J］. Computer Methods in Applied Mechanics and Engineering， 2005， 194（39/40/41）：4135-4195.

［6］ WANG Huiping， WANG Dongdong. Efficient Meshfree Computation with Fast Treatment of Essential Boundary Conditions for Industrial Applications［J］. Journal of Engineering Mechanics， 2009， 135（10）：1147-1154.

［7］ CHEN J S， WANG Huiping. New Boundary Condition Treatments in Meshfree Computation of Contact Problems［J］. Computer Methods in Applied Mechanics and Engineering， 2000， 187（3/4）：441-468.

［8］ GINGOLD R A， MONAGHAN J J. Smoothed Particle Hydrodynamics：Theory and Application to Non-spherical Stars［J］. Monthly Notices of the Royal Astronomical Society， 1977， 181（3）：375-389.

［9］ WANG Lihua， HU Minghao， ZHONG Zheng， et al. Stabilized Lagrange Interpolation Collocation Method：a Meshfree Method Incorporating the Advantages of Finite Element Method［J］. Computer Methods in Applied Mechanics and Engineering， 2023， 404：115780.

［10］ LUCY L B. A Numerical Approach to the Testing of the Fission Hypothesis［J］. The Astronomical Journal， 1977， 82：1013.

［11］ 黃志強. 無網格法中本質邊界條件實施研究［D］. 西安：西北工業大學， 2007.

HUANG Zhiqiang. Research on the Implementation of Essential Boundary Conditions in Meshless Method［D］.Xi’an：Northwestern Polytechnical University， 2007.

［12］ WANG Dongdong， XUAN Junchang. AnImproved NURBS-based Isogeometric Analysis with Enhanced Treatment of Essential Boundary Conditions［J］. Computer Methods in Applied Mechanics and Engineering， 2010， 199（37/38/39/40）：2425-2436.

［13］ 軒軍廠.基于改進邊界條件施加方式和應變光滑子域積分的幾何精確NURBS有限元分析［D］. 廈門：廈門大學，2010.

XUAN Junchang. Geometrically Accurate NURBS Finite Element Analysis Based on Improved Boundary Condition Application and Strain-smooth Subdomain Integral［D］. Xiamen：Xiamen University，2010.

［14］ 陳濤， 莫蓉， 萬能. 等幾何分析中Dirichlet邊界條件的配點施加方法［J］. 機械工程學報， 2012， 48（5）：157-164.

CHEN Tao， MO Rong， WAN Neng. Imposing Dirichlet Boundary Conditions with Point Collocation Method in Isogeometric Analysis［J］. Journal of Mechanical Engineering， 2012， 48（5）：157-164.

［15］ 王東東， 軒軍廠， 張燦輝. 幾何精確NURBS有限元中邊界條件施加方式對精度影響的三維計算分析［J］. 計算力學學報， 2012， 29（1）：31-37.

WANG Dongdong， XUAN Junchang， ZHANG Canhui. A Three Dimensional Computational Investigation on the Influence of Essential Boundary Condition Imposition in NURBS Isogeometric Finite Element Analysis［J］. Chinese Journal of Computational Mechanics， 2012， 29（1）：31-37.

［16］ 張漢杰， 王東東， 軒軍廠. 薄梁板結構NURBS幾何精確有限元分析［J］. 力學季刊， 2010， 31（4）：469-477.

ZHANG Hanjie， WANG Dongdong， XUAN Junchang. Non-uniform Rational B Spline-based Isogeometric Finite Element Analysis of Thin Beams and Plates［J］. Chinese Quarterly of Mechanics， 2010， 31（4）：469-477.

［17］ HOANG T， VERHOOSEL C V， AURICCHIO F， et al. Skeleton-stabilized Iso Geometric Analysis：High-regularity Interior-penalty Methods for Incompressible Viscous Flow Problems［J］. Computer Methods in Applied Mechanics and Engineering， 2018， 337：324-351.

［18］ LEI Zhen， GILLOT F， JEZEQUEL L.AC0/G1 Multiple Patches Connection Method in Isogeometric Analysis［J］. Applied Mathematical Modelling， 2015， 39（15）：4405-4420.

［19］ EMBAR A， DOLBOW J， HARARI I. Imposing Dirichlet Boundary Conditions with Nitsche’s Method and Spline-based Finite Elements［J］. International Journal for Numerical Methods in Engineering， 2010， 83（7）：877-898.

［20］ CHEN Tao， MO Rong， GONG Zhongwei. Imposing Essential Boundary Conditions in Isogeometric Analysis with Nitsche’s Method［J］. Applied Mechanics and Materials， 2011， 121/122/123/124/125/126：2779-2783.

［21］ NGUYEN V P， KERFRIDEN P， BRINO M， et al. Nitsche’s Method for Two and Three Dimensional NURBS Patch Coupling［J］. Computational Mechanics， 2014， 53（6）：1163-1182.

［22］ HU Qingyuan， CHOULY F， HU Ping， et al. Skew-symmetric Nitsche’s Formulation in Isogeometric Analysis：Dirichlet and Symmetry Conditions， Patch Coupling and Frictionless Contact［J］. Computer Methods in Applied Mechanics and Engineering， 2018， 341：188-220.

［23］ DU Xiaoxiao， ZHAO Gang， WANG Wei， et al. Nitsche’s Method for Non-conforming Multipatch Coupling in Hyperelastic Isogeometric Analysis［J］. Computational Mechanics， 2020， 65（3）：687-710.

［24］ 胡清元. 等幾何分析中的閉鎖問題與Nitsche方法研究［D］. 大連：大連理工大學， 2019.

HU Qingyuan. On the Locking Problem and the Nitsche’s Method in Isogeometric Analysis［D］.Dalian：Dalian University of Technology， 2019.

［25］ DORNISCH W， KLINKEL S. Boundary Conditions and Multi-patch Connections in Isogeometric Analysis［J］. PAMM， 2011， 11（1）：207-208.

［26］ TEMIZER I， WRIGGERS P， HUGHES T J R. Contact Treatment in Isogeometric Analysis with NURBS［J］. Computer Methods in Applied Mechanics and Engineering， 2011， 200（9/10/11/12）：1100-1112.

［27］ MITCHELL T J， GOVINDJEE S， TAYLOR R L. A Method for Enforcement of Dirichlet Boundary Conditions in Isogeometric Analysis［M］∥MUELLER-HOEPPE D， LOEHNERT S， REESE S. Recent Developments and Innovative Applications in Computational Mechanics. Heidelberg：Springer， 2011：283-293.

［28］ JIANG Kai， ZHU Xuefeng， HU Changzhi， et al. An Enhanced Extended Isogeometric Analysis with Strong Imposition of Essential Boundary Conditions for Crack Problems Using B++ Splines［J］. Applied Mathematical Modelling， 2023， 116：393-414.

［29］ 王勖成. 有限單元法［M］. 北京：清華大學出版社， 2003.

WANG Xucheng. Finite Element Method［M］. Beijing：Tsinghua University Press， 2003.

［30］ 黃艾香， 周天孝. 有限元理論與方法—第一分冊［M］. 北京：科學出版社， 2009.

HUANG Aixiang， ZHOU Tianxiao. Finite Element Theory and Methods-Division 1［M］. Beijing：Science Press， 2009.

［31］ ATLURI S N，SHEN Shengping.The Meshless Method［M］. Tech. Science Press Encino， 2002.

［32］ MONTEMURRO M. On theStructural Stiffness Maximisation of Anisotropic Continua under Inhomogeneous Neumann-Dirichlet Boundary Conditions［J］. Composite Structures， 2022， 287：115289.

［33］ MONTEMURRO M， RODRIGUEZ T， PAILHS J， et al. On Multi-material Topology Optimisation Problems under Inhomogeneous Neumann-Dirichlet Boundary Conditions［J］. Finite Elements in Analysis and Design， 2023， 214：103867.

［34］ De PRENTER F， VERHOOSEL C V， Van BRUMMELEN E H， et al. Stability and Conditioning of Immersed Finite Element Methods：Analysis and Remedies［J］. Archives of Computational Methods in Engineering， 2023， 30（6）：3617-3656.

（編輯 王旻玥）

基金項目：國家重點研發計劃（2023YFB2504605）；國家自然科學基金（52205267）

作者簡介：

王書亭，男，1968年生，博士、教授。研究方向為復雜機電產品集成設計理論與技術。E-mail：wangst@hust.edu.cn。

謝賢達*（通信作者），男，1994年生，博士、副研究員。研究方向為等幾何分析與拓撲優化。E-mail：xiexd020@ncu.edu.cn。

本文引用格式：

王書亭，謝晴天，楊奧迪，等.等幾何拓撲優化非齊次邊界條件施加技術研究［J］. 中國機械工程，2025，36（3）：525-535.

WANG Shuting， XIE Qingtian， YANG Aodi， et al. Research on Non-homogeneous Boundary Condition Imposing Techniques for Isogeometric Topology Optimization［J］. China Mechanical Engineering， 2025， 36（3）：525-535.