電泳過程車身結構變形仿真方法概述

李鼎 于保君 劉創舉 朱學武 曹正林

（1. 中國第一汽車股份有限公司 研發總院，長春130013；2. 汽車振動噪聲與安全控制綜合技術國家重點實驗室，長春130013）

主題詞：涂裝電泳 有限元法 有限體積法 流固耦合

縮略語

FVM Finite Volume Method

FSI Fluid-Structure Interaction

VOF Volume Of Fluid

FEM Finite Element Method

FSI Fluid-Structure Interaction

CFD Computational Fluid Dynamics

CSD Computational Structural Dynamics

1 前言

電泳涂裝工藝一般由涂裝前預處理、電泳涂裝、電泳后清洗、電泳涂膜的烘干等多道工藝組成，如圖1所示。為提高電泳質量和生產節拍，新型翻轉線、穿梭機等電泳線逐步被采用，而與之帶來的是白車身閉合件的變形、生產過程中車身傾斜和運動裝置的脫鉤等現象在電泳過程中時有發生。

圖1 電泳涂裝工藝示意

傳統結構強度設計主要考慮整車的碰撞、強度、模態和剛度參數，并未充分考慮白車身在各大工藝過程的應力應變累積情況。而白車身在沖壓、焊接工藝設計時，充分考慮了材料、模具、夾具因素的影響，白車身焊接總成整體變形量不大。而當焊接總成進過涂裝生產線，由沖壓、焊接引入的殘余應力，電泳泳池中的流體壓力積累的塑性變形，烘烤引入的材料應力釋放、熱應力變形，最終造成總裝時部分車身結構面差較大，無法順利完成裝配。

由于白車身經過涂裝后，變形影響因素復雜，因此，在研究此類系統問題時，首先需要將問題進行分解，然后分析不同工序對結構變形的影響，找出設計的薄弱環節，加以優化，再通過多學科仿真工具將全生命周期下的變形情況進行分析和驗證。

白車身在電泳槽中的受力主要來自其在泳池中改變姿態而引起的阻力、浸入液面以下受到的浮力、自身重力、工裝的約束力，因此可以將模型抽象為求解流體力學及結構力學問題。流體力學部分需要考慮氣-液-固3相的交互作用，目前考慮多相流的方法主要是有限體積法（Volume Of Fluid,VOF），由于白車身在泳池中是運動的，因此若采用網格算法進行仿真，則需要采用重疊網格技術。考慮到流體-固體間的交互作用并不顯著，因此計算時可以采用單向耦合法進行計算，即先計算不同位置、時間下白車身各零部件表面受到的流體壓力情況，再將其作為壓力邊界條件導入到結構力學計算中。由于白車身鈑金為標準固體模型，目前金屬彈塑性變形數值模擬中的一種有效數值計算方法為有限元法，將求解對象離散為由各種單元組成的計算模型，計算流體載荷激勵下各離散單元非線性位移響應，從而獲取鈑金殘余變形。

當前對電泳過程車身結構變形仿真方法研究尚處于起步階段。ESS斯太爾工程軟件公司提出了一種基于有限差分法對工業EPD 涂層過程進行準確有效的模擬方法[1]。埃爾朗根-紐倫堡大學的Dominik Bartuschat 介紹了一種用于微流控電泳大規模并行直接數值模擬的耦合算法，這種多物理算法采用了流體和離子的歐拉描述，結合了移動帶電粒子的拉格朗日表示[2]。卡爾斯魯厄理工學院用數值方法證明了分段式對電極在不同電位差下可以有效地減少涂層材料的用量，將現有的模型擴展到三維空間中的任意非平面幾何[3]。上汽通用五菱汽車股份有限公司應用電泳仿真CAE技術對車身電泳膜厚進行了仿真分析，提升了產品的強度和防腐性能[4]。Institut für Computerphysik利用分子動力學模擬方法研究了聚電解質在平行非帶電平板間的電泳拉伸，模擬結果顯示，流體動力相互作用的影響，減少了限制板之間的距離[5]。

2 流體側仿真

電泳過程車身結構流固耦合仿真主要采用計算流體力學計算流體部分。白車身在電泳槽中運動，會產生復雜多變的流體流動特性，電泳液的流動符合湍流流動特性，而車身運動中包含進出槽過程，當中會卷入空氣進入車身腔體內部，電泳液在車身表面存在浸潤、物質交換的物理化學過程，攜帶的電泳液在出槽過程會有流淌現象。因此白車身電泳槽的運動過程實質上是一個多相流問題。根據具體工程應用，可以將電泳液定義為流體的液相，空氣定義為氣相，整個流動過程符合氣液2相流計算模型。

2.1 有限體積法網格化理論

有限體積法（Finite Volume Method，FVM）是近年發展非常迅速的一種離散化方法，其特點是計算效率高，目前在流體領域得到了廣泛的應用。其基本思路是：將計算區域劃分為網格，并使每個網格點周圍有一個互不重復的控制體積；將待解的微分方程（控制方程）對每一個控制體積分，從而得到一組離散方程。其中的未知數是網格點上的因變量，為了求出控制體的積分，必須假定因變量值在網格點之間的變化規律。從積分區域的選取方法看來，有限體積法屬于加權余量法中的子域法，從未知解的近似方法看來，有限體積法屬于采用局部近似的離散方法。北京航空航天大學流體力學研究所的朱自強給出了一種生成分區對接網格和Euler 方程分區解，以二維多段翼型為對象分別討論了分區重疊結構網格方法、非結構網格方法和自適應的笛卡爾網格法[6]。楊建宏通過有限元空間和有限體積元空間的一種雙射投影，得到了不可壓縮流問題低次等階穩定有限體積元方法[7]。D.G.Roychowdhury 討論了計算流體力學中的各種控制方程及其推導，以及偏微分方程描述邊界條件的物理和數學意義[8]。

對于電泳過程車身高度復雜的幾何，捕捉此幾何可能比較困難。在這類情況下，包面可用于簡化和捕捉復雜幾何的水密表示。無論是否使用包面，通常都使用表面重構來重構幾何的初始表面。幾何的初始表面通常由許多三角表面組成。此三角形化也稱為網格化，通常包含高度偏斜的三角形，不適用于生成高質量體網格。重構通過生成尺寸更均勻的三角形提高這些表面的整體質量，此方法最適合網格模型。然后，使用重構表面作為體網格的基礎。

2.2 多相流仿真方法

目前在多相流計算中，如果采用網格算法，通常采用歐拉-歐拉法和歐拉-拉格朗日法2種方法，其中歐拉-歐拉法常用流體體積模型（VOF）、混合物模型及歐拉模型，歐拉-拉格朗日法常用離散相模型。

VOF 模型是一種在固定歐拉網格下表面追蹤的方法，它通過研究網格單元中流體和網格體積分數來確定自由面，追蹤流體的變化，而非追蹤自由液面上質點的運動。在計算過程中，利用各流體的體積分數為計算參數，在一個動量方程中進行求解，把每個時刻各流體所占的體積比都記錄下來。

體積追蹤/交界面捕捉法非常適用于模擬具有大尺度交界面的分層多相流體。必須注意的是，計算網格單元的尺寸在將交界面分類為大尺度還是小尺度時至關重要。由于體積追蹤/交界面捕捉法依賴于完全求解交界面，因此在離散流態建模精度同樣重要的多尺度流體模擬中使用這些方法會導致成本過高。流體體積(Volume Of Fluid,VOF) 多相模型實施屬于交界面捕捉方法系列，可預測不混溶相交界面的分布和移動。此建模方法假設網格分辨率足以求解相之間的交界面位置和形狀。

交界面的相分布和位置由相體積分數∝i的場來描述。相i的體積分數定義如式（1）：

式中，Vi為為網格單元中相i的體積，V為網格單元的體積。

網格單元中所有相的體積分數總和必須是式（2）：

式中，N為總相數；∝i=0表示網格單元完全沒有相i；∝i=1 表示網格單元完全由相i填充；0＜∝i＜1表示2個極限之間的值存在相間交界面i。

包含交界面的網格單元中計算的材料屬性依賴于組成流體的材料屬性。同一包含交界面的網格單元中的流體被視為混合物，其中混合物的密度、粘度及比熱如式（3）（4）（5）：

式中，ρi為相i密度，μi為相i動力粘度，(Cp)i相i比熱。

體積分數傳輸方程相i的分布由相質量守恒方程驅動，如式（6）：

式中，a為表面積矢量，v為混合（質量平均）速度，vd,t為擴散速度，Sai為相i的用戶自定義源項為相密度ρi的材料或拉格朗日導數。

計算采用如下所示計算相體積分數：

1）存在2 個VOF 相時，僅對第1 個相求解體積分數傳輸。在每個網格單元中，調整第2個相的體積分數，使這2個相的體積分數總和等于1。

2）存在3個或更多VOF相時，對所有相求解體積分數傳輸。然后，根據每個網格單元中的所有相的體積分數總和標準化每個相的體積分數。

連續性方程所有相的總質量守恒方程由式（7）給出：

式中，S為與相源項相關的質量源項，如式（8）所示：

通過由式（3）給出的密度來解釋對流體混合物組成相的體積分數的依賴性。

動量方程及能量方程如式（9）式（10）所示：

式中，p為壓力，I為單位張量，T為應力張量，fb為體積力矢量。

工程中關注的大多數流體流都具有不規則的波動流量。此類波動通常尺度小但頻率高，要在時間和空間上對其求解，需要高昂的計算成本。可以不求解湍流的精確控制方程（直接數值模擬），改為求解平均量或濾波量來近似小波動結構的影響，這樣成本較低。各種湍流模型為此類結構的建模提供了不同的方法。

2.3 重疊網格技術

由于白車身在電泳槽中是運動的，如果采用網格重構技術，網格計算量很大而且會造成后期結構有限元計算網格節點不統一。因此需要使用重疊網格技術。重疊網格也稱作“嵌套”網格，用于離散具有以任意方式相互重疊的多個不同網格的計算域。這些網格在處理涉及多個體或移動體的問題時以及優化研究中最為有用。Tapan K. Sengupta 通過計算二維Navier-Stokes 方程，使用單塊結構和重疊網格來求解圓柱繞流，利用相對較少的網格點捕獲了該布置的已知流動特性[9]。

中國船舶科學研究中心的趙發明采用RANS方法和重疊網格計算了帶自由液面的船舶繞流問題，描述了重疊網格和單相Level Set 自由液面模擬方法的數學模型及求解過程[10]。上海交通大學船舶海洋與建筑工程學院的王建華采用重疊網格技術開發的多功能水動力學求解器naoe-FOAM-SJTU，對標準船模DTMB5512裸船體在平面運動機構(PMM)控制作用下對純搖首運動進行了數值模擬[11]。海軍工程大學的周廣禮以Suboff標模為研究對象，基于RANS方程、VOF模型及重疊網格技術，建立了潛艇應急上浮運動的數值預報方法[12]。

在大多數情況下，使用重疊網格不需要在生成初始網格之后進行任何網格修改，比標準網格化方法更靈活。任何涉及重疊網格的研究都具有封閉整個求解域的背景區域，以及包含域中的體的一個或更多較小區域（圖2）。

圖2 重疊網格示意

在重疊網格中，網格單元分組為活動、不活動或受體網格單元。在活動網格單元中，將對離散控制方程進行求解。在不活動網格單元中，不會對任何方程進行求解，但是如果重疊區域在移動，這些網格單元可以變為活動狀態。受體網格單元可分離背景區域中的活動和不活動網格單元，并與重疊區域中的重疊邊界相連接。受體網格單元用于耦合2個重疊網格上的求解。一個網格的供體網格單元處的變量值，通過插值表示其他網格中受體網格單元處的變量值。供體網格單元是其他網格中最靠近受體網格單元的活動網格單元。

將同時為所有區域中的所有活動網格單元計算求解，即網格為隱式耦合。當引用離散方程中一個網格的受體網格單元中的變量值時，將使用其他網格中供體網格。

3 固體側仿真

電泳過程中，車身變形是由液體沖擊及車身自身運動導致，因此需將流體側計算獲得的不同時刻下的流體壓力值，映射到結構計算的有限元網格上，作為壓力邊界條件計算整個過程中白車身的彈塑性變形情況。

3.1 有限元法網格化理論

有限元法（Finite Element Method,FEM）是廣泛應用的固體力學數值計算方法。有限元法是將一個連續的求解域任意分成適當形狀的許多微小單元，并于各小單元分片構造插值函數，然后根據極值原理（變分或加權余量法），將問題的控制方程轉化為所有單元上的有限元方程，把總體的極值作為個單元極值之和，即將局部單元總體合成，形成嵌入了指定邊界條件的代數方程組，求解該方程組就得到各節點上待求的函數值。有限元法的基礎是極值原理和劃分插值，它吸收了有限差分法中離散處理的內核，又采用了變分計算中選擇逼近函數并對區域積分的合理方法，是這2 類方法相互結合，取長補短發展的結果。它具有廣泛的適應性，特別適用于幾何及物理條件比較復雜的問題，而且便于程序的標準化。Nathan Mendes 在Basics in Practical Finite-Element Method 一書中詳細介紹了有限元法及其應用[13]。莊茁在基于ABAQUS 的有限元分析和應用中詳細介紹了有限元的基本理論和數值計算方法，并對相關例題進行了討論研究[14]。杜平安詳細介紹了劃分有限元網格時應考慮的一些基本原則，包括網格的數量、疏密、階次、質量、分界點、布局、編號等[15]。

3.2 多體動力學仿真

白車身在電泳槽中的運動模式與涂裝工藝相關，如圖3所示。在自身重力的作用下，不同的運動模式對車身變形有相應的影響，因此在模擬電泳過程固體側變形時，需考慮動力學影響。

圖3 車身電泳過程示意

應用拉格朗日待定乘子法，動力學方程如式（11）：

式中，Φq為約束方程雅克比矩陣；M為質量矩陣；λ為拉格朗日乘子；q?為加速度；F(q,q?,t)為作用力；γ為加速度公式中二階導數項。

其中，系統約束方程為式（12）：

求約束方程廣義坐標的偏微分，可得到約束方程雅克比矩陣Φq：

速度方程為式（14）

加速度方程為式（15）

由加速度方程得到的拉格朗日乘子γ為式（16）

3.3 結構力學非線性仿真

與流體方程不同，固體應力在拉格朗日架構中求解，網格依據材料中的材料點。這意味著特定網格單元中的質量將在變形時保持不變，因為它在物理上代表相同的材料。此網格單元質量由材料原始密度與原始體積的乘積來確定。因此，變形狀態下的密度是原始質量除以網格單元的當前體積。

固體力學用于描述固體連續體對應用負載的響應行為。應用負載包括體積力、表面負載、點力或固體溫度變化導致的熱負載。應用負載會在結構中產生應力場，并且可能導致結構位移從初始未變形配置到變形配置。

應用的負載可導致固體結構產生從初始配置到變形配置的位移。總位移由剛體運動和體中各點的相對位移（決定固體結構的變形）共同組成，如圖4所示。

圖4 各節點總位移示意[16]

如果未變形配置中材料點的位置為X，且此點相對于變形配置的位移為u(X,t)，則該材料點在變形配置中的位置為式（17）：

采用分量形式，位移可以表示為式（18）：

剛體的位移場完全由單個位移矢量定義。可變形體的位移場由其材料點的位移矢量集定義。

應力是作用于表面的單位面積力的度量。在固體結構中，外力或溫度變化會產生可能導致固體結構運動和變形的應力場。而阻礙結構變形的內力也會引起可將結構恢復到原始未變形狀態的內部彈性應力。一些材料還具有內置應力，這種應力在沒有作用力和變形時就已存在。

一般情況下，作用于體的平面截面上的應力由矢量（稱作應力矢量或拉力）定義如式（19）：

式中，F為作用于平面的力，A為平面的面積。

在任意點處，應力是單位面積上的力，應力狀態完全由與3個通過該點且相互垂直的平面關聯的應力矢量定義。因此，任何平面上的任何點處的應力狀態均通過式（20）形式的二階張量定義：固體側計算得到的離散網格節點坐標與原始坐標之差即為車身在電泳過程中的殘余變形。

4 流固耦合仿真方法研究

流固耦合(Fluid-Structure Interaction,FSI)是目前很多領域研究的熱點和難點之一，與國外相比，國內在該方面的研究相對滯后或片面。Lúcia Armiliato Sangalli 提出了一種用于顫振失穩主動控制橋梁數值模擬的流固耦合格式，采用2 步Taylor-Galerkin 顯式方法和任意拉格朗日-歐拉（ALE）描述求解了流動基本方程[17]。Andreas Hessenthaler 建立了16 個獨立的分析解決方案，能夠在二維和三維、靜態和瞬態情況以及線性和超彈性固體材料中實現法FSI 分析，證明了這些解決方案在分析收斂行為方面的實用性[18]。陳鋒紹了流固耦合的基本理論和算法，包括流固耦合問題的基本概念與耦合機理的分類，探討了強耦合和弱耦合算法的計算流程與適用范圍，闡述了任意拉格朗日歐拉法、無網格法、特征線分離算法及基于流動條件插值的有限元法等流固耦合分析方法[19]。王建軍詳細描述了自由液面大晃動的流固耦合問題的性質和特點，討論了交替求解方法及其對求解此類問題的特別適應性[20]。Yue Qian-bei 提出了一種基于子迭代回路的弱耦合算法，對彈性管的固有頻率與靜流頻率之比的影響進行了深入的數值研究[21]。José A.González提出了一種考慮沸騰和冷凝相變的水蒸汽能量系統瞬態模擬的流固耦合有限元方法，利用Crank-Nicolson 顯式隱式時間積分方法，提升了分析精度和穩定性[22]。

4.1 有限元法網流固耦合分類

流固耦合問題從物理上可分為2 類：第1 類是耦合僅作用在2相交界面上，由耦合面上的動力學平衡條件及運動學協調條件來引入方程上的耦合，如氣動彈性、水動彈性問題等；第2類是二相遇部分或全部重疊在一起，如滲流問題，其耦合效應通過描述問題的微分方程體現。流固耦合問題從控制方程解法上可分為直接求解(Direct)的強耦合(Strong coupling)和分區迭代求解(Partitioned)的弱耦合(Loose coupling)。所謂強耦合，是將流體域、固體域和耦合作用構造在同一控制方程中，在同一時間步內，同時求解所有變量。由于固體方程與流體方程存在很大的差異，聯立求解困難重重。目前還沒有一款商業軟件可以求解強流固耦合問題。弱耦合是在每一時間步內，分別依次對CFD方程和CSD方程求解,通過中介交換固體域和流體域的計算結果數據,從而實現耦合求解。目前，流固耦合基本上都是采用弱耦合。由于存在時間差，所以與現實情況存在一定的誤差。單向耦合與雙向耦合主要是針對弱耦合求解。

考慮到車身電泳過程流體-固體間的交互作用并不顯著，且車身結構較為復雜，因此，計算時采用單向弱耦合法進行計算，即先計算不同位置、時間下白車身各零部件表面受到的流體壓力情況，再將其作為壓力邊界條件導入到結構力學計算中，最終獲得白車身的整體彈塑性變形結果。

4.2 分類界面信息傳遞與網格匹配

流固耦合的關鍵是界面信息傳遞，其原理很簡單，即流體域和固體域在邊界上對應質點間滿足位移協調ds=df和作用力的平衡條件fs=ff。但進行數值模擬時，由于CFD 和CSD 對計算網格要求的不同，導致耦合界面上的網格不匹配，節點多數情況下并不重合，從而不能直接進行數據的交換，這就需要流固界面上節點采用一定的信息映射方法。

通常CFD計算要求的網格密度比CSD要密得多，由此產生耦合界面上2套非匹配網格之間的數據傳遞問題,在數學上是一個雙向插值問題。結構將在流體荷載作用下發生的變形通過耦合界面由CSD 網格傳遞給CFD 網格，CFD 求解器在新狀態下求解流場，再將流體荷載通過耦合界面傳遞給CSD 網格。可以通過采用一定的映射算子來進行計算，將流體節點應力首先映射到固體節點，然后通過積分(采用固體邊界的差值函數)至固體節點力，再將計算出固體模型的節點力，在有限元步驟里作為自然邊界條件。由于固體節點力是流體速度和壓強的函數，因此固體域的結構離散方程里包括了流域求解變量。

5 結束語

綜上所述，電泳過程車身結構變形仿真分析尚處于理論階段，未有全面系統的流程與試驗對比作為支撐，少數關于電泳過程車身的仿真仍主要以涂層本身為研究對象，且主要為靜態過程。本文所述的電泳過程車身結構變形預期仿真方法會在以下3點存在較大挑戰：

（1）車身電泳過程中邊界條件，包括工藝操作參數設定的整體移動速度、翻轉速度、噴嘴流場均會影響白車身水動力CFD分析結果，且車身在出水時會產生白車身集液問題，因此邊界條件的準確度及選用的算法模型將會對結果精度有較大影響。

（2）由于白車身在電泳槽中是運動的，因此與采用重疊網格技術，而白車身與開閉件結構較為復雜且有較多孔隙，結構的離散化規則會導致不同的結果響應，且會影響相應的計算時間。

（3）由于流體側與固體側分別采用有限體積法與有限元法進行求解，因此分類界面信息傳遞與網格匹配中，不同的映射算子、傳遞迭代步長、耦合策略等將會對固體域的結構離散方程求解有較大影響，如何優化上述參數需要進行深入研究。