1. 項目概述CRTSⅡ型軌道精細模型的定位這幾年做高鐵無砟軌道結構分析繞不開的一個問題就是怎么把列車移動荷載真實地加到有限元模型里而不是用靜態(tài)荷載糊弄過去。用Abaqus建立CRTSⅡ型軌道精細模型配合DLOAD子程序實現(xiàn)列車移動荷載的施加是目前工程和科研里比較主流的做法。這套方案既能拿到軌道板、底座板、CA砂漿層內部的應力應變細節(jié)又能反映移動荷載下結構動力響應的真實規(guī)律比單純在靜力分析里加一個固定力要靠譜得多。這個項目適合誰參考一類是做無砟軌道結構設計驗算和疲勞評估的工程師另一類是研究車輛-軌道耦合動力學、軌道結構損傷演化的研究生和科研人員還有剛接觸Abaqus二次開發(fā)、想搞清楚DLOAD子程序怎么用于移動荷載的同行。不管你屬于哪一類這套從建模到子程序再到后處理的完整流程都是可以直接落地復用的。我下面寫的內容全部來自一個實際項目的完整復現(xiàn)過程涉及幾何尺寸、材料參數(shù)、網格策略、子程序代碼和各類報錯排坑一步步來。2. 建模方案與結構簡化2.1 CRTSⅡ型軌道的結構組成與幾何參數(shù)CRTSⅡ型板式無砟軌道和CRTSⅠ型、CRTSⅢ型最大的區(qū)別在于軌道板是連續(xù)結構沿縱向有預應力筋貫穿中間沒有斷開的板縫。整體從下往上依次是底座板、CA砂漿調整層、軌道板、扣件系統(tǒng)和鋼軌。這種結構形式在高速鐵路線上很常見尤其是長軌鋪設的區(qū)段整體性好、剛度均勻對高速運行的列車來說平順性優(yōu)勢很明顯。建模前先把尺寸理清楚。鋼軌采用CHN60軌高176mm底部寬150mm扣件采用WJ-8B型小阻力扣件間距一般在630mm左右軌道板厚度一般取200mm寬度2550mmCA砂漿層厚度約30mm底座板厚度約300mm比軌道板寬一些具體寬度和配筋方式根據(jù)線下基礎類型略有差異。需要說明的是不同線路的局部尺寸會有調整建模前先拿到設計圖確定具體數(shù)值別直接抄本文的數(shù)值否則算出來跟實際結構對不上。如果做整段區(qū)間計算量會非常大。一般做法是取軌道結構的一段代表性長度來建模比如含4到5塊軌道板范圍的連續(xù)段總長10到20m然后施加電周期性邊界或者用黏彈性邊界模擬無限長的軌道基礎。我這里取了一個18m長的軌道段模型包含了完整的底座板、砂漿層、兩塊連續(xù)軌道板以及對應范圍的扣件和鋼軌。2.2 材料參數(shù)與單位體系的選擇Abaqus本身沒有固定的單位制需要自己保證量綱一致。軌道結構分析常用mm-N-s-tonne這套單位體系長度用mm力用N時間用s質量用tonne對應的應力單位就是MPaN/mm2密度單位是tonne/mm3。這套單位的優(yōu)勢在于和工程圖紙尺寸直接兼容不需要做數(shù)值換算。材料參數(shù)按常見工程取值來給。軌道板是C50混凝土彈性模量3.55×10? MPa泊松比0.2密度2.5×10?? tonne/mm3。底座板C30或C40混凝土彈性模量取3.0×10?到3.25×10? MPa。CA砂漿層是彈模相對低很多的材料一般在100到300 MPa之間泊松比0.2左右這層是整個軌道結構里比較薄弱的環(huán)節(jié)也是車致?lián)p傷分析的重點關注對象。鋼軌用鋼E2.06×10? MPa泊松比0.3密度7.85×10?? tonne/mm3?？奂到y(tǒng)用彈簧阻尼單元模擬豎向剛度一般在20到80 kN/mm之間阻尼20到40 kN·s/m。提示如果你的模型單位體系不同注意把速度單位也換算過來。比如mm-s單位下300 km/h要寫成83333.3 mm/s寫錯單位是最常見的低級錯誤。2.3 網格劃分與單元類型選擇網格策略是精細模型的關鍵取舍點直接決定計算精度和成本。鋼軌如果只關心荷載傳遞效果可以用梁單元但要想看到鋼軌截面內部接觸應力就得用實體單元。這里我建議用實體單元建模鋼軌軌頭接觸區(qū)網格加密到2到3mm軌底過渡區(qū)5到8mm這樣移動荷載下鋼軌的彎曲應力和接觸區(qū)應力都能算得比較準。軌道板和底座板都是大體積混凝土構件用C3D8R六面體縮減積分單元常規(guī)區(qū)域網格尺寸30到50mm扣件區(qū)域局部加密到10mm左右。CA砂漿層厚度只有30mm沿厚度方向至少要劃分2到3層單元否則彎曲響應出不來。網格數(shù)量控制在50萬到100萬之間配合隱式求解和子程序計算單次動力響應分析可在可接受的時間內完成。單元類型選C3D8R時要注意沙漏控制尤其是在沖擊荷載作用下縮減積分單元容易產生沙漏變形增加單元網格密度或使用增強沙漏控制選項可以緩解。軌道板受彎為主的受力模式下也可以考慮C3D20R二次單元提高彎曲精度但計算成本會明顯上升權衡下來C3D8R加合理網格密度是更劃算的方案。3. DLOAD子程序的編寫與核心實現(xiàn)3.1 DLOAD的工作機制與調用原理DLOAD是Abaqus標準求解器中用于施加與位置、時間相關分布荷載的接口。每個帶分布荷載的積分點都會在增量步開始時調用一次子程序Abaqus把當前計算時間的積分點坐標傳遞給子程序子程序根據(jù)坐標和時間判斷這個點是否處于荷載作用范圍內并返回對應的荷載值F。F的單位是壓強即力除以面積所以最終施加在鋼軌頂面上的荷載是通過一個壓力帶的形式來模擬移動輪載的。搞清楚這個機制后實現(xiàn)移動荷載的思路就很清晰了每個時刻計算一組輪對的空間位置然后判斷鋼軌頂面各積分點是否落在輪載作用區(qū)段內。落在區(qū)段內的點施加相應的壓力值區(qū)段外的點返回零。如果采用隱式求解這個荷載會跟隨時間步長更新從空間上看就是一組沿鋼軌縱向移動的壓力帶。一個需要提前說明的地方是DLOAD作用范圍是按當前荷載作用區(qū)域來識別的每個積分點只能屬于有荷載或者無荷載兩種狀態(tài)。如果直接寫一個if判斷荷載區(qū)段首尾會出現(xiàn)階躍變化也就是壓力從零瞬間跳到滿值這種突變在隱式求解里容易造成收斂困難。所以實際編寫時要給荷載區(qū)段的邊界設置過渡區(qū)域用線性過渡或平滑過渡的方式讓壓力值漸變上升和下降計算穩(wěn)定性會好很多。3.2 單個輪對移動荷載的數(shù)學表達先解決單輪對的實現(xiàn)問題。假設列車沿軌道縱向X方向行駛初始時刻輪對位于X0處運行速度為V那么在時間T時輪對的位置可以寫成Xwheel X0 V × T需要明確的是DLOAD子程序里TIME(1)指的是分析步的累計時間。多分析步情況下如果荷載只在第二個分析步開始施加要注意TIME(1)的基準點是從當前分析步還是整個分析開始計算可以先用寫入外部文件的方式確認一下時間基準避免出現(xiàn)荷載位置錯位的問題。在Abaqus代碼里判斷當前積分點是否在荷載作用區(qū)內的邏輯如下如果ABS(COORDS(1) - Xwheel)小于等于荷載分布半長L/2那么這個積分點位于荷載帶內F取為輪載壓強否則F取0。輪載壓強怎么定假設軸重為14噸單個輪載為70kN即7×10? N。荷載分布區(qū)長200mm鋼軌頂面荷載作用寬度按50mm估算壓力帶的承載面積就是200×5010000mm2對應壓強F70000/100007 MPa。這個壓強值在上述單位制下恰好為7 N/mm2。完整的單輪對DLOAD子程序如下SUBROUTINE DLOAD(F,KSTEP,KINC,TIME,NODE,NOEL,NPT,LAYER, 1 KSPT,COORDS,JLTYP,SNAME) C INCLUDE ABA_PARAM.INC C DIMENSION COORDS(3), TIME(2) CHARACTER*80 SNAME C REAL*8 V, X0, XWHEEL, LZONE, PRESS, DLOADWIDTH C C 參數(shù)定義 V 83333.3D0 ! 列車速度 mm/s對應300km/h X0 1000.0D0 ! 初始輪對位置 mm LZONE 200.0D0 ! 荷載分布區(qū)長度 mm DLOADWIDTH 50.0D0 ! 荷載分布區(qū)寬度 mm PRESS 70000.0D0 / (LZONE * DLOADWIDTH) ! 輪載壓強 N/mm2 C C 當前時刻輪對位置 XWHEEL X0 V * TIME(1) C C 判斷積分點是否在荷載作用區(qū)內 IF (DABS(COORDS(1) - XWHEEL) .LE. LZONE / 2.0D0) THEN F PRESS ELSE F 0.0D0 END IF C RETURN END這里把荷載定義成一個200mm長的均布帶實際車輪鋼軌接觸斑沿縱向也就十幾毫米200mm是一個等效分布長度目的是在網格尺寸不小的情況下也能保證荷載帶覆蓋至少一個單元的積分點避免壓力只落在個別單元上導致局部應力失真。網格越密這個分布長度可以越接近真實接觸斑尺寸。3.3 多個輪對與整列車荷載的擴展實現(xiàn)高速列車一個轉向架帶兩輪對兩輪對之間軸距通常為2.5m一節(jié)車廂兩端的轉向架中心距約17.5m實際編組車里每個輪對的絕對位置都會隨時間變化。擴展寫法是把所有輪對初始位置存成數(shù)組每個輪對按同樣的速度移動然后任意積分點只要落在任何一個輪對的作用區(qū)段內就施加對應的輪載壓強。具體實現(xiàn)里要注意一軸兩端輪對分布在兩根鋼軌上每根鋼軌只承擔左側或右側的輪載所以子程序施加時按鋼軌位置區(qū)分。如果模型的鋼軌編號和坐標固定可以直接在子程序里判斷COORDS(2)或者COORDS(3)來區(qū)分是哪根鋼軌。多輪對擴展的推薦寫法是用循環(huán)例如REAL*8 DIST(4) DATA DIST /0.0D0, 2500.0D0, 17500.0D0, 20000.0D0/ C xBase X0 V * TIME(1) F 0.0D0 C DO I 1, 4 XWHEEL X0 V*TIME(1) - DIST(I) IF (DABS(COORDS(1) - XWHEEL) .LE. LZONE/2.0D0) THEN F PRESS GOTO 100 END IF END DO C 100 CONTINUE RETURN END注意DIST數(shù)組存的是相對首輪對的偏移量xBase相當于首輪對的當前位置后面每個輪對的位置在此基礎上減去偏移量。一個常見的錯誤是直接把所有輪對的絕對位置寫死這樣車一動起來輪對之間的間距就不對了。3.4 子程序的編譯驗證與調試技巧在Abaqus中使用子程序前先確認Fortran編譯環(huán)境和Abaqus版本匹配。過一遍這個流程安裝Intel Fortran Compiler和Microsoft Visual Studio配置好環(huán)境變量后在命令行執(zhí)行abaqus verify -user_std如果顯示successful則說明編譯鏈路是通的。這個驗證步驟不要跳過否則經常在提交任務時報一堆找不到編譯器的錯誤浪費時間又查不到根因。調試DLOAD子程序最直接的辦法是在子程序里把關鍵變量寫入外部文件比如每調用一次就記錄當前節(jié)點坐標、時間、計算出的XWHEEL和F值。我在實際調試中是把這些信息寫入一個文本文件然后導入Excel里檢查荷載帶的位置時序是否與理論值一致。這個方法雖然笨但往往幾分鐘就能定位到問題。另外一個調試技巧是先做一個靜態(tài)驗證把速度設為0讓輪對固定在一個位置提交一個靜力分析步看看鋼軌變形和應力分布是否對稱合理。對稱性檢查能快速發(fā)現(xiàn)模型坐標系錯誤、荷載作用位置偏移等問題比直接上動態(tài)分析好查得多。4. 邊界條件、接觸設置與求解控制4.1 層間接觸與約束策略CRTSⅡ型軌道層間連接是建模中影響結果很大的環(huán)節(jié)。鋼軌和扣件之間、扣件和軌道板之間采用彈簧阻尼單元連接一般用Spring2/Dashpot2單元單獨建立扣件系統(tǒng)替代實際的扣件部件。這樣做的好處是可以通過調整彈簧剛度和阻尼參數(shù)來模擬不同扣件類型比如WJ-8B和WJ-7型差異就直接改參數(shù)不需要重新建模。軌道板與CA砂漿層、CA砂漿層與底座板之間可以采用綁定約束Tie來簡化處理。但如果研究目標是軌道板與砂漿層的離縫損傷就必須用帶損傷本構的界面單元或者面面接觸這樣才能模擬層間拉應力超過粘結強度后的脫開行為。這個選擇取決于你的研究目標不要盲目追求精細。用Tie約束時要注意主面和從面的網格密度協(xié)調從面網格應比主面細一些或至少相當否則約束面上會出現(xiàn)應力集中和偽振蕩。CA砂漿層本身是薄弱層在Tie處理后雖然不會脫開但應力結果相對均勻適合做整體響應分析。4.2 邊界條件的合理截斷軌道結構的縱向尺度遠遠大于建模范圍如果直接把有限長度的模型兩端約束死會產生嚴重的邊界效應移動荷載接近端部時結果失真。正確處理辦法是采用半無限域近似或者黏彈性邊界。最簡單的方案是把底座板底面固結在長度方向兩端外側再加一段過渡底座板并在端面施加彈性地基彈簧來模擬周圍土體與相鄰結構的約束作用。在動力學計算中還可以在端部加黏性邊界即通過阻尼單元模擬能量的逸散避免反射波在模型里來回彈跳導致結果振蕩。對18m長的模型我給底座板底面全部固結縱向兩端設置彈性彈簧彈簧剛度根據(jù)地基系數(shù)和等效面積估算。如果做的是具體線路評估最好按實際線下基礎條件來標定這部分參數(shù)。4.3 分析步設置與求解器參數(shù)隱式分析里移動荷載是強非線性輸入分析步參數(shù)設置直接影響收斂性和計算效率。建議把分析步設置成固定增量步長一般取荷載帶走過一個單元長度所需時間的1/5到1/10。比如網格尺寸20mm速度83333mm/s走過一個單元需要0.00024s增量步取2×10??到5×10??s比較合適。時間增量步太大荷載跳變劇烈容易不收斂步長太小計算時間成倍增加。我實際試算下來的經驗是先用一個較粗的網格和較大的增量步跑通全流程確認結果合理后再加密網格并細化步長不要一上來就追求極限精度。阻尼方面需要特別注意。軌道結構的實際阻尼遠小于一般建筑結構瑞利阻尼的Alpha和Beta參數(shù)要根據(jù)結構自振頻率來標定不要隨便取默認值?？梢韵茸瞿B(tài)分析獲取軌道結構的一階豎向彎曲頻率再用頻率值反推阻尼系數(shù)。質量阻尼Alpha對低頻響應影響大剛度阻尼Beta對高頻振蕩影響大給得太高會把高頻輪軌動力響應抹平給得太低又會出現(xiàn)數(shù)值振蕩需要反復對比。5. 常見報錯與排坑實錄5.1 CPU數(shù)量超過許可限制的報錯Abaqus 提交并行任務時報錯“the number of cpus (20) exceeds the number of cpus available”是很多新手容易卡殼的地方搜索量也一直很高。這個錯誤核心原因是兩種一是求解器分配的CPU數(shù)量超過了當前許可證允許的核數(shù)二是軟件讀取的系統(tǒng)邏輯核數(shù)與實際可用核數(shù)不符。后者在Windows系統(tǒng)下比較常見比如虛擬機環(huán)境只分配了部分邏輯核Abaqus卻識別到了更多。處理方法在Job模塊點擊Edit把Parallelization里的CPU數(shù)量改小一般先設2或者4跑通流程確認沒問題再逐步增加。同時可以命令行執(zhí)行abaqus informationlicenses查看許可證授權的核心數(shù)。如果是虛擬機或遠程桌面環(huán)境檢查系統(tǒng)CPU親和性設置是否限制了Abaqus實際可用的核心數(shù)。5.2 安裝后無桌面啟動文件和許可證不能啟動“Abaqus安裝后桌面沒有啟動程序文件”和“Abaqus許可證不能啟動”屬于安裝配置階段的高頻問題。桌面沒有快捷方式一般不是安裝失敗而是安裝程序沒默認創(chuàng)建圖標。解決辦法是找到安裝目錄下的啟動腳本比如CAE的bat文件或者Exec文件夾里的abaqus.bat直接雙擊運行或者手動創(chuàng)建快捷方式指向該腳本。還有一種情況是CAE啟動時依賴的Python環(huán)境路徑配置錯誤檢查環(huán)境變量PYTHONHOME是否被其他軟件改寫。許可證不能啟動的原因比較多集中在幾個方向許可證服務器服務沒起來、環(huán)境變量LM_LICENSE_FILE和服務器的端口設置不對、防火墻阻斷了Abaqus License Server的通信。先確認許可證服務已在服務管理器里啟動再用命令行執(zhí)行abaqus licensing檢查當前許可證狀態(tài)。如果服務器是遠程的確認客戶端環(huán)境變量里填寫的端口和主機名與服務器設置一致注意端口號必須和服務器配置的端口完全匹配。5.3 DLOAD子程序編譯與運行期故障DLOAD子程序最常見的編譯錯誤是找不到Fortran編譯器。Abaqus版本和Intel編譯器版本之間兼容性要求很強不匹配就會出現(xiàn)“cannot find ifort”或類似的報錯。先執(zhí)行abaqus verify -user_std驗證整套編譯鏈路這是最快定位問題的方法。如果驗證失敗對照Abaqus官方兼容性表格重新安裝匹配的編譯器版本。運行期還有一個很隱蔽的問題DLOAD子程序里的局部變量沒有初始化。Fortran中未初始化的局部變量在不同編譯環(huán)境下可能是隨機值導致F輸出異常。強烈建議子程序入口處把F默認為0所有局部變量顯式賦值。這類問題排查起來特別耗時因為模型網格、材料參數(shù)都沒問題但荷載就是不對。5.4 荷載帶階躍導致的計算不收斂這個問題前面提到過但值得專門拿出來說。DLOAD子程序用if判斷實現(xiàn)的荷載帶邊界上是從0直接跳到滿值在隱式求解器中容易造成應力波傳播異常和收斂迭代次數(shù)激增。我實際遇到的案例是同樣的模型和材料參數(shù)加了漸變過渡的荷載帶后計算時間縮短了一半還多而且結果更平滑。推薦做法是把判斷條件改成按相對位置計算過渡系數(shù)比如離荷載帶中心越遠荷載值按線性或余弦曲線遞減到0。在子程序里用一個過渡半寬定義比如過渡區(qū)取20mm那么F PRESS × max(0, 1 - |x - xwheel - LZONE/2| / TRANSWIDTH)。注意單獨處理荷載帶前后兩個邊界不要寫死對稱邏輯否則頭尾過渡不對稱。6. 結果解讀與模型驗證經驗6.1 軌道板與鋼軌動力響應判讀移動荷載算完之后第一步是看鋼軌的豎向位移時程。單輪荷載下鋼軌最大動位移一般在1到2mm量級如果速度提高后位移明顯增大且伴隨高頻振蕩說明輪軌動力作用增強結果在物理上說得通。軌道板的彎曲應力重點關注板底受拉區(qū)因為無砟軌道損傷最容易從板底開裂開始。后處理時沿軌道板縱向取幾條路徑輸出彎矩應力的分布對比不同時刻的應力峰值位置可以識別出列車輪載作用下軌道板的受荷循環(huán)特征。這里注意區(qū)分靜載作用和動載沖擊作用產生的應力增量動載沖擊導致的應力增幅一般在10%到30%之間如果遠超這個范圍要檢查阻尼參數(shù)是否給得過大。6.2 模型驗證的幾條關鍵指標模型做出來對不對不能只看云圖顏色好看得有對照依據(jù)。幾條可用的驗證路徑第一與理論解析解對比鋼軌在集中力作用下的彈性彎曲位移可以用Winkler地基梁公式估算對比有限元結果和理論值的偏差如果超過10%優(yōu)先檢查扣件剛度和網格密度第二與文獻中類似參數(shù)的CRTSⅡ型軌道實測數(shù)據(jù)對比重點關注軌道板加速度峰值區(qū)間和鋼軌動位移范圍第三做收斂性驗證用兩套不同粗細的網格計算同一工況如果關鍵響應偏差在5%以內說明網格密度足夠。6.3 模型擴展的方向與建議這套模型框架可以非常方便地擴展。想研究鋼軌波磨與輪軌力的關系可以修改子程序里的輪載表達式把車輪扁疤或軌道不平順的影響加進去。想分析CA砂漿層離縫擴展可以把砂漿層單元換成內聚力模型給界面一個損傷起始強度和斷裂能配合DLOAD移動荷載反復掃掠就能模擬疲勞累積損傷的演化過程。想考慮橋上無砟軌道則需要在底座板下方增加橋梁梁段和支座的建模計算量會再次提升但方法完全一致。我在實際做這個項目時最大的體會是不要把精力全放在追求模型“多精細”上而是先明確研究問題需要的精度等級。純粹算整體動力響應CA砂漿層簡化成Tie就能得到很好的結果要研究層間損傷就必須上內聚力模型和精細網格。工具就擺在那里關鍵是舍得花時間在子程序調試和模型驗證上這兩個環(huán)節(jié)做扎實了后面出結果和分析都是水到渠成的事。