1、板料成形數值模擬的關鍵技術,Engineering Technology Associates (China) , Inc,Contents,單元的選擇 材料的定義 模具與板料的相互作用接觸的定義 自適應網格細分 (Adaptive Meshing) 顯式算法和隱式算法 時間步長 一步算法，增量算法和2.5D 成形極限圖（FLD）,LSDYNA的單元類型,LS-DYNA的單元類型包括: 質量單元 慣性單元 梁單元 薄殼單元 (沖壓模擬分析所用的單元) 對于一個方向相比其他方向薄得多的部件來說,一般采用薄殼單元來劃分網格,采用的是平面應力條件,即在單元法向的應力為0 厚殼單元 體單元 彈簧阻尼單

2、元 SPH單元 安全帶單元,EQ.1: Hughes-Liu 殼單元 LSDYNA中最早的殼單元，面內單點積分，能對殼單元的參考面進行偏置，2倍于缺省殼單元時間開銷 EQ.2: Belytschko-Tsay (default) 殼單元 缺省殼單元公式，面內單點積分，計算速度很快，對于大變形問題是最穩(wěn)定有效的公式，建議在大多數的分析中使用 EQ.3: BCIZ triangular shell 三角形殼單元 根據Kirchhoff薄板理論導出，面內單點積分 EQ.4: C0 triangular shell 三角形殼單元 根據mindlin-Reissner薄板理論導出,面內單點積分,比BCI

3、Z三角形殼單元精度好 EQ.5: Belytschko-Tsay membrane 膜單元 沒有彎曲力矩和橫向剪切力,面內單點積分,建議在大多數薄膜分析中使用 EQ.6: S/R Hughes Liu 選擇性縮減積分殼單元 面內采用2X2積分點方式,沒有沙漏模式,比缺省公式慢10倍 EQ.7: S/R co-rotational Hughes Liu選擇性縮減積分殼單元 面內采用2X2積分點方式,沒有沙漏模式,比缺省公式慢5倍,Shell Elements 殼單元,單元方程,EQ.8: Belytschko-Leviathan 殼單元 面內單點積分,比缺省單元慢1.25倍,自動含有物理沙漏控制

4、 EQ.9: fully integrated membrane 全積分BT膜單元 面內采用2X2積分點,不需要沙漏控制,建議用于安全氣囊分析 EQ.10: Belytschko-Wong-Chiang 殼單元 面內單點積分,用于單元過渡翹曲情況,對于大變形問題沒有BT殼穩(wěn)定 EQ.11: Fast Hughes-Liu 快速殼單元 面內單點積分,有沙漏模式,比缺省公式慢1.5倍 EQ.12: Plane stress 2D element (x-y plane) 平面應力殼單元 EQ.13: Plane strain 2D element (x-y plane) 平面應變殼單元 EQ.14:

5、 Axisymmetric Petrov-Galerkin 2D solid 軸對稱體,面積加權 EQ.15: Axisymmetric Galerkin 2D solid 軸對稱體,體積加權 EQ.16: fully integrated 全積分殼單元 BT殼單元公式,采用2X2積分點方式,只比BT殼單元多23倍,可以處理翹曲的幾何體問題,此時需要激活第八種沙漏控制公式.,Shell Elements 殼單元,單元方程,為什么提供這么多單元方程? 全積分單元用于彈塑性，金屬成形，氣囊以及關心精確度的時候。 三角形單元(C0單元）剛度比較硬，不建議采用。但是在混合網格中采用，因為C0三角形殼單

6、元，比退化的四邊形網格算法好。 膜單元不能受彎曲和斷面剪應力，適用于非常薄的板料以及拉張為主的變形中 。 進化的板殼單元有較高的精確度，但是和BT單元相比穩(wěn)定性較差 BT單元速度非?？?！,單元方程,Characteristics of the Belytschko-Tsay SE (default) 基于隨動坐標系和速度位移方程. 隨動坐標系的單元方程, 避免了在單元中嵌入坐標系而導致的非線性動態(tài)的復雜性. 單元動量方程假設節(jié)點位于同一個平面. 如果是5個積分點, BT單元需要 725數學計算。但是一點HL單元需要4066次數學計算. HL單元的選擇減縮積分需要35,367次數學計算。 由于計

7、算速度效率非常高, BT單元方程通常是殼單元的單元方程的最佳選擇. 正因為如此, 它已經成為4節(jié)點板殼單元的缺省的單元方程。 BT單元簡單，速度非?？?。,單元方程,The cost of shell elements vary by formulation,單元方程,單元屬性,*SECTION_SHELL $PROPERTY NAME:blankpro $ SID ELFORM SHRF NIP PROPT QR/IRID ICOMP 5 2 0.833E+00 3.0 1.0 0.0 $ T1 T2 T3 T4 NLOC 1.000E+00 1.000E+00 1.000E+00 1.000

8、E+00,如果成形完后還要計算回彈，需選擇此項,參數選擇,為什么不用統一的參數? LSDYNA是通用求解器 不同的應用需要不同的參數 不當參數影響結果預測 取決于用戶的結果時常發(fā)生 對于板料成形，我們已經找到了合適的參數 做了非常多的分析實例 參數一致性得到長足改進 對不同的應用采用適當的參數，有利于分析效率 大多數參數，DYNAFORM都設定了合理的值 用戶的經驗仍然重要 缺省參數有其局限性 新的應用，需要新的參數 可以用來調試，,在鈑金成形仿真中通常使用的單元方程 #2單元方程（ Belytschko-Tsay） #16單元方程 (全積分單元方程) #2單元在平面內一個積分點 #16單元在

9、平面內四個積分點,#2 單元 #16單元,參數選擇 單元方程（1）,參數選擇 單元方程（2）,在成形分析中 平面內拉伸重要 對于減薄和拉伸，使用面內一個積分點不會導致明顯誤差 在回彈預測中 應力分布更為重要 彎曲和非彎曲, 在單元平面內應力是不一致的 在面內選擇更多的積分點是必要的,參數選擇 單元方程（3）,特點: Element #2 比 #16快 (35倍) Element #16 能夠捕獲合理的應力分布 建議: 對成形分析#2單元是比較好的選擇 節(jié)省CPU 降低內存需求 對回彈分析, element #16是更好的選擇 更加精確的應力分布 在關鍵字中選擇單元方程 *SECTION_SHE

10、LL ELFORM=2: Element #2 is selected ELFORM=16: Element #16 is selected (ELFORM is the second parameter under the card of *section_shell),厚度方向上的積分點(NIP) 在厚度方向的應力分布是復雜的 需要更多的積分點來捕獲在厚度方向不同的應力分布模式.,幾種在厚度方向的應力模式 Case I: 純彈性彎曲; Case II: 彈性和塑性彎曲; Case III: 彎曲 + 拉伸; Case IV: 回彈后的,參數選擇 厚度上的積分點（1）,NIP = 2 Cas

11、e I中沒有錯誤; Case II, III, IV中，導致錯誤 NIP = 3 Case I中沒有錯誤; case II, III中，比較小的偏差; Case IV中，較大偏差 NIP=5,7 Case I中沒有錯誤; Case II, III, and IV中，比較小的偏差.,NIP = 3 Real stress,參數選擇 厚度上的積分點（2）,NIP 7 精確度沒有明顯提高 增加不必要的CPU 需要更多的內存 不推薦 Suggestions: 重力載荷 (可以認為是純彈性變形) NIP = 2 可以增加計算速度 成形分析: (Case II and III) NIP=3 回彈分析: (

12、Case IV) NIP=5, 7 注意: 參數NIP 在 *SECTION_SHELL,參數選擇 厚度上的積分點（3）,單元使用提示,單元算法BT(2#)，BWC(10#)，BL(8#)采用面內單點積分，全積分采用四點積分 所有殼單元厚度方向可以用任意多個積分點 對彈性變形沿厚度方向用兩個積分點即可 塑性行為沿厚度方向用3到5個積分點，回彈可以用7個積分點 避免使用小單元，以免縮小時間步長。如果使用，請同時使用質量縮放。 減少使用三角形單元。 避免銳角單元和翹曲單元，否則會降低計算精度。,單向拉伸試驗驗證的塑性行為,材料模型/塑性行為,單軸拉伸試驗,在單軸拉伸下, 可以看到下列階段 彈性 屈

13、服 塑性硬化 頸縮 破裂,材料模型,Option 1: Bilinear Model 線性硬化模型 Option 2: Exponential Model 冪指數硬化模型 Option 3: Experimental Data 試驗曲線硬化模型,鈑金成形的材料塑性模型:,exponential,材料,Material 18: *Mat_Power_Law_Plasticity Material 24: *Mat_Piecewise_Linear_Plasticity Material 36: *Mat_3-Parameter_Barlet Material 37: *Mat_Transvers

14、ely_Anisotropic_Elastic_Plastic Material 1: *Mat_Elastic Material 64: *Mat_Rate_Sensitive_Power_Law_Plasticity Material 103: *Mat_Anisotropic_Viscoplastic Material 33: *Mat_Barlet_Yld96,材料的加工硬化, = K n,s,從實驗數據得到硬化曲線（決定彈塑性材料的屬性）,假設我們從拉伸實驗得到力的曲線，用下列步驟得到硬化曲線 產生工程應力/應變 產生真實應力/應變 去除彈性應變得到等效應力應變曲線,STEP 1:

15、Engineering Stress/Strain Curve,工程應力/應變曲線,工程應力=-,力,原始面積,工程應變=-,長度改變量,原始長度,STEP 2: True Stress/Strain Curve,真實應力/應變曲線,真實應力=-,力,當前面積,真實應變=-,長度改變量,當前長度,STEP 3: Shift to Remove Elastic Strains,有效應力/應變曲線,真實應力應變曲線去除彈性應變得到有效應力應變曲線,硬化曲線,所得到的硬化曲線是屈服應力和等效塑性應變之間的函數關系。 對金屬來說，真實應力是遞增的,*Mat_3-Parameter_Barlet,*MA

16、T_3-PARAMETER_BARLAT $MATERIAL NAME:DQ $ MID RO E PR HR P1 P2 1 7.850E-09 2.070E+05 2.800E-01 2.000E+00 5.320E+02 2.000E-01 $ M R00 R45 R90 LCID E0 SPI 6.000E+00 1.870E+00 1.270E+00 2.170E+00 0 0.000E+00 0.000E+00 $ AOPT 0.0 $ XP YP ZP A1 A2 A3 $ V1 V2 V3 D1 D2 D3,由Barlat和Lian開發(fā) 平面應力狀態(tài)下的各向異性材料 由Lank

17、ford參數定義各向異性. 有三種硬化準則可以選擇. 比37號材料計算慢,*Mat_Transversely_Anisotropic_Elastic_Plastic,*MAT_TRANSVERSELY_ANISOTROPIC_ELASTIC_PLASTIC $MATERIAL NAME:CQ $ MID RO E PR SIGY ETAN R HLCID 1 7.850E-09 2.070E+05 2.800E-01 1.760E+02 5.860E+02 1.400E+00 0.000E+00,s,p,使用Hill屈服函數. 只考慮橫斷面各向異性,平面內R=(R0+2R45+R90)/4 雙

18、線性或任意由定義的應力/等效塑性應變曲線. 不能模擬出凸耳現象,材料庫,材料參數曲線,通?？赡艿腻e誤 觀測結果: 在鈑金模擬過程中,過早失效(頸縮) 發(fā)現原因: 應力應變曲線不合理(高應變范圍太低),Bad stress-strain curve,Effective strain,Effective stress,材料參數曲線,產生原因: 通過實驗只能獲得一定范圍內的應力應變曲線,所以需要外插. 有幾種外插方式,Commonly extrapolated stress-strain curves,Constant slope (too high stress in large strain),

19、Power law (more reasonable fit),No hardening (result in earlier necking),使用材料提示,在定義材料屬性時，確保使用一致的單位制，不正確的單位制(密度,時間,應力)不僅影響材料響應，而且影響到接觸剛度的計算。 大多數高度非線性有限元分析的精確取決于所使用的材料參數。要得到好的分析結果，需要使用精確的材料參數。所以不要低估準確材料數據對結果的重要性，盡量花費額外的時間和金錢去獲得準確的材料數據。,接觸問題,在有限元分析中，接觸問題的處理往往是衡量有限元軟件分析能力的一個重要指標。LS_DYNA有40多種接觸類型供用戶選擇，具有

20、強大的接觸分析能力。之所以有這么多的接觸類型，一是由于有一些專門的接觸類型用于專門應用，另一個是由于有一些老的接觸類型一直保留，主要是為了使以前建立的有限元模型能一直使用。 選擇合適的接觸類型和定義接觸參數對于初學者來說可能有一定的困難。 主（Master）與從（Slave）的概念 在接觸定義中有主面與從面的概念，一旦定義了一個接觸對，程序檢查從面的節(jié)點是否與主面的段發(fā)生接觸。 從面考慮的是節(jié)點，所以對從面來說有限元模型是不是連續(xù)不是很重要，一般要求它的網格劃分比主面更密。 主面上考慮的是段，可以是連續(xù)的網格也可以是不連續(xù)的網格。所以一般來說： 粗網格表面定義為主面，細網格表面定義為從面；主從

21、面相關材料剛度相處懸殊，材料剛度大的一面為主面；平直或凹面為主面，凸面為從面。,接觸算法,Lsdyna處理接觸問題一般采用3種算法： 動力約束法 分配參數法 對稱罰函數法（缺省算法） 動力約束法是最早采用的接觸算法，1976年最先用于dyna2d程序。它的原理是：在每一時步修正構形前，檢查從節(jié)點是否貫穿主表面，并調整時間步大小，使那些貫穿從節(jié)點都不貫穿主表面，對所有已經和主表面接觸的從節(jié)點施加約束條件，保持從節(jié)點與主表面接觸；另外檢查與主表面接觸的從節(jié)點所屬單元是否存在受拉界面力，如有則用釋放條件使從節(jié)點脫離主表面。由于該算法比較復雜，目前僅用于固連接觸。 分配參數法用于有相對滑動沒有分離的滑

22、動處理，如炸藥爆炸的氣體對結構的壓力作用。其原理是：將每一個正在接觸的從單元一半質量分配到被接觸的主表面面積上，同時由每個從單元的內應力確定作用在接受質量的主表面面積上的分布壓力。在完成質量和壓力的分配后，程序修正主表面的加速度，然后對從節(jié)點的加速度和速度施加約束，保證從節(jié)點沿主表面運動。程序不允許從節(jié)點穿透主表面。,對稱罰函數法是LS-DYNA的缺省算法，其原理是：每一時步先檢查各從節(jié)點是否穿透主表面，沒有穿透則不對從節(jié)點做任何處理。如果穿透，則在該從節(jié)點與主表面間、主節(jié)點和從表面間引入一個較大的界面接觸力，大小與穿透深度、接觸剛度成正比，成為罰函數值。其物理意義相當于在其中放置一系列法向彈

23、簧，限制穿透。 接觸力由下面公式計算： FK 公式中K為接觸界面剛度（接觸界面剛度由單元尺寸和材料特性確定） 為穿透量,接觸的分類： LS_DYNA有40多種接觸類型，可以大致分類如下： 單向 雙向 單面 實體 固連,基于罰函數法，基于約束法(若有_CONSTRAINT),固連約束法,單向與雙向接觸的區(qū)別 單向的概念：僅檢查從節(jié)點對主面的穿透，計算效率高，使用單向接觸的接觸類型如下.,NODES_TO_SURFACE AUTOMATIC _NODES_TO_SURFACE FORMING_NODES_TO_SURFACE CONSTRAINT_NODES_TO_SURFACE ERODING_

24、NODES_TO_SURFACE 定義一些單個節(jié)點作為從接觸面 不對稱的 不是任意選擇主從接觸面 通常來講, 網格較粗的面作為主接觸面 Nodes-to-surface 接觸計算效率高 尤其適合點接觸剛體 使用單向接觸有時可能捕捉不到接觸行為,單向接觸,與單向接觸不同的是,雙向接觸既檢查從節(jié)點對主面的穿透,又檢查主面 節(jié)點對從面的穿透,這樣對于主從表面的定義是任意的,代價是處理接觸方面的計算量將增加2倍左右.使用雙向接觸的類型有: SURFACE_TO_SURFACE AUTOMATIC_SURFACE_TO_SURFACE FORMIMG_SURFACE_TO_SURFACE CONSTRA

25、INT_SURFACE_TO_SURFACE ERODING_SURFACE_TO_SURFACE,雙向接觸,厚度偏置,“厚度偏置”是兩個接觸面之間的偏置距離, 接觸面即單元中面 所有的 “automatic”接觸包含厚度偏置 對非自動類型的接觸，厚度偏置是可選的。,兩種厚度,單元厚度 在 *section_shell or *element_shell_thickness 影響單元剛度和單元質量 接觸厚度 在接觸中考慮單元厚度等距 不影響單元剛度和單元質量 缺省的接觸厚度 = 單元厚度 可以在*contact or *part_contact中設置接觸厚度,projected contact

26、 surface,初始接觸穿透 (Interference) 導致主要問題 計算不穩(wěn)定 LS-DYNA 盡力更正幾何形狀 在第一次計算的時候，把穿透的節(jié)點移到接觸面上。,接觸/穿透,成形類型接觸（forming-),沖頭和模具定義為目標面（主接觸面），工件定義為接觸面（從接觸面） 對于這些接觸類型中，模具網格無需貫通，因此減少接觸定義的復雜性 模具網格的方向必須一致 成形接觸選項基于自動接觸類型，因此功能十分強大,拉延筋（Drawbead）接觸,Drawbead contact通常用于板料成形，用于約束板料的運動 在板料成形過程中，通常會出現工件和模具失去接觸（如起皺） Drawbead co

27、ntact容許使用彎曲和摩擦阻力，用于確保工件在整個沖壓過程中與模具始終保持接觸。,CONTACT Card,*CONTACT_FORMING_ONE_WAY_SURFACE_TO_SURFACE $ CID CONTACT INTERFACE TITLE $ 2 BLANK/PUNCH $ SSID MSID SSTYP MSTYP SBOXID MBOXID SPR MPR 1 2 2 2 $ FS FD DC VC VDC PENCHK BT DT .125E+00 .000E+00 .000E+00 .000E+00 .200E+02 00.0000e+001.0000E+20 $ S

28、FS SFM SST MST SFST SFMT FSF VSF .000E+00 .000E+00 .000E+00 .000E+00 $ SOFT SOFSCL LCIDAB MAXPAR PENTOL DEPTH BSORT FRCFRQ 0 $ PENMAX THKOPT SHLTHK SNLOG 1 動靜態(tài)摩擦系數 接觸阻尼（VDC）：在接觸面垂直方向增加阻尼，實際頻率轉換為臨界頻率的百分比，用于消除振蕩 小的穿透 接觸厚度 接觸生死時間,參數選擇 接觸類型選擇（1）,很多種接觸類型滿足不同的應用需求 正確的接觸非常重要，不正確的接觸會導致: 大的穿透 沙漏模式 較長的CPU時間 計

29、算本身失敗 壞的仿真結果 對初學者來說，接觸的選擇是困難的,對鈑金成形分析，一般采用三種接觸 Type 1: CONTACT_FORMING_NODES_TO_SURFACE Type 2: CONTACT_FORMING_ONE_WAY_SURFACE_TO_SURFACE Type 3: CONTACT_FORMING_SURFACE_TO_SURFACE 第一二中算法一致，只是輸入文件不一樣 第三種同等對待板料和剛體，耗費更多的時間. 第三種接觸有時剛性過強，導致接觸不穩(wěn)定。 第二中在沖壓分析中最常使用,參數選擇 接觸類型選擇（2）,成形接觸特點 向前網格細分 剛體沒有厚度（如果相應參數

30、為+） 如果為負,工具容許接觸等距. 工具網格可以不連續(xù) 網格方向自動檢查 板料是從接觸面 工具是主接觸面 為每一個工具定義一個接觸界面,參數選擇 接觸類型選擇（3）,參數選擇 接觸類型選擇（5）,懲罰接觸理論 Placing normal interface springs between all penetrating nodes and the contact surface.,Contact force is proportional to the penetration,參數選擇 接觸類型選擇（6）,接觸不當造成的穿透,Large penetration cause bad resu

31、lt 大的穿透將引起不好的結果 得不到模具的細節(jié) Correlation with experimental measurement will not be good Over-predict sheet metal formability Large safety margin should be used 錯誤的應力分布 Cause bad springback prediction,參數選擇 接觸類型選擇（7）,參數選擇 接觸類型選擇（8）,Problem: Higher interface pressures results in unacceptable penetration So

32、lution: 減小時間步長 (DT2MS in *CONTROL_TIMESTEP) Scaling up the stiffness (SLSFAC in *CONTROL_CONTACT) 減小工具速度 較長的CPU時間 單元尺寸也許太大 網格細分頻率也許太低,智能網格細分的優(yōu)點,需要足夠小的來捕獲產品細節(jié) 太小的網格耗費更多的CPU和內存 智能網格細分可以選擇性的細分網格 初始網格較大可以采用較大的初始時間步.,智能網格細分的缺陷,缺點 每次細分都要I/O,降低計算速度. 太粗初始網格導致artificial bending and unbending,Too coarse initi

33、al mesh cause artificial deformation The following mesh refinement can not correct the error,*CONTROL_ADAPTIVE $ ADPFREQ ADPTOL ADPOPT MAXLVL TBIRTH TDEATH LCADP IOFLAG 4.0909E-04 .150E+02 2 40.0000e+001.0000E+20 1 $ ADPSIZE ADPASS IREFLG ADPENE ADPTH MEMORY ORIENT MAXEL 2.000 1 0 1.000 .500 0 0,網格細

34、分間隔太大,大單元進入成形半徑; ADPSIZE=0.4Rmin; LCADP自定義細分頻率;,Adaptive Meshing,顯式算法和隱式算法,顯式時間積分 沖擊, 穿透, 較大的動態(tài)效應 保證穩(wěn)定狀態(tài)需要小的時間步 運算的具體問題決定使用的最大時間步 隱式時間積分 靜態(tài),特征值，較低的動態(tài)效應 較大的時間步 線性求解 (剛度矩陣) 當K線性時，無條件穩(wěn)定，可以用較大的時間步 當非線性問題，用一席雷線性逼進來獲得求解（收斂需要較小的時間步，高度非線性問題，收斂困難）,隱式算法,隱式算法的優(yōu)點 類靜態(tài)和低頻反應問題 - 靜態(tài)加載的彈性問題計算成本小 使用大的時間步來使非線性問題求解順利 -

35、 比隱式需要較少的時間步 缺點 高度非線性問題求解非常困難甚至步收斂 - 材料非線性和失效 每個時間步消耗的CPU比較高 LS-Dyna Explicit/Implicit solver LS-DYNA version 950+ has an embedded implicit solver,時間步長,在求解過程中，LS-DYNA計算所需時間步時，檢查所有單元，用最小單元來決定時間步長。 其中N是單元數量。 為達到穩(wěn)定，采用一個比例因子，通產采用缺省值是.90或更小的值。 為了減少求解時間，可以使用較大的，保證穩(wěn)定性的時間步長。,Timestep,時間步長,臨界（或最?。r間步長: 其中C聲音

36、在3D連續(xù)介質中傳播的速度,Timestep,時間步長,聲波在不同材料中傳播速度例子 空氣 331 m/s 水 1478 m/s 鋼 5240 m/s 鋁 5328 m/s 鈦 5220 m/s 樹脂玻璃 2598 m/s,聲音在兩種最常見的金屬鐵，鋁中的傳播速度范圍是5000.m/s. 在兩種材料中,單元特征值是5.mm， 時間步長是1.s, 通常被設為行業(yè)標準。 合理的時間步要求最小的單元尺寸是5.mm, 車身的幾何形狀不能用有限單元方法完全表達出來。這是整車模型中限制最大的問題: 最小單元尺寸限制幾何細節(jié)或者充分的網格密度。,Forming Limit Diagram,BIAXIAL S

37、TRETCH (O-A Line): The sheet is stretched by equal amounts in both directions. The strains and stresses are maj = min , t = -2maj (in the thickness direction) maj = min 0, t = 0 PLANE STRAIN (O-C Line): The sheet is stretched in only one direction. The strain and stress states are (For isotropic material) maj = -t (in the thickness direction), min = 0 maj min 0 , t = 0 UNIXIAL TENSION (O-B Line): The sheet is stretched in one direction while it shrinks in the ot