簡(jiǎn)介面向鐵路輪軌接觸力學(xué)研究場(chǎng)景的Python計(jì)算資源適用于具備有限元或接觸力學(xué)基礎(chǔ)、希望突破經(jīng)典赫茲理論局限的工程師與科研人員。資源包內(nèi)提供基于Piotrowski-Kik模型的簡(jiǎn)化數(shù)值實(shí)現(xiàn)覆蓋滑動(dòng)接觸、黏著力學(xué)、滾動(dòng)接觸疲勞等非赫茲因素可用于分析輪軌接觸力分布、磨損趨勢(shì)與疲勞壽命。包內(nèi)共13個(gè)文件以Python腳本為核心包含幾何計(jì)算、模型求解與主程序模塊另有RST/Markdown文檔用于說明模型原理與使用方式附輪軌型面數(shù)據(jù)文件整體約49KB結(jié)構(gòu)緊湊、便于二次開發(fā)。目前已有265人學(xué)習(xí)下載適合作為教學(xué)演示、算法驗(yàn)證或初步工程評(píng)估的參考基座。通過調(diào)整參數(shù)與結(jié)合實(shí)測(cè)數(shù)據(jù)可進(jìn)一步完成模型校核與優(yōu)化幫助理解高速重載工況下的輪軌動(dòng)態(tài)響應(yīng)特征。1. 非赫茲問題的輪軌接觸力學(xué)簡(jiǎn)化模型為什么不能繼續(xù)用赫茲解做車輛動(dòng)力學(xué)仿真時(shí)赫茲解對(duì)磨耗后的輪軌接觸往往對(duì)不上這時(shí)候你需要的正是非赫茲問題的輪軌接觸力學(xué)簡(jiǎn)化模型一個(gè)用 Python 寫成的可下載代碼包。這類模型的作用是在不滿足赫茲假設(shè)的輪軌工況下把接觸斑形狀和法向壓力分布快速算出來精度介于赫茲解與有限元/邊界元之間速度卻比有限元快一到兩個(gè)數(shù)量級(jí)。適合做多體動(dòng)力學(xué)仿真、輪軌磨耗預(yù)測(cè)和批量參數(shù)掃描的人也適合想把接觸算法改寫成自己模塊的開發(fā)者。先說明白一件事赫茲解不是不能用而是它的適用條件很嚴(yán)輪軌工程里有一大批工況恰好落在條件之外。2. 簡(jiǎn)化模型的理論邊界與選型邏輯赫茲失效的三個(gè)場(chǎng)景和兩條技術(shù)路線2.1 赫茲接觸為什么在輪軌場(chǎng)景下失效赫茲接觸解大家都很熟接觸斑為橢圓、壓力分布為半橢球、最大壓力等于平均壓力的三倍。這些結(jié)論漂亮但都建立在一組強(qiáng)假設(shè)上——接觸體是彈性半空間、表面無(wú)摩擦、接觸區(qū)附近的幾何能用二次曲面描述、接觸斑尺寸遠(yuǎn)小于物體尺寸。輪軌接觸在新廓形、小蠕滑工況下勉強(qiáng)滿足這些假設(shè)這也是為什么傳統(tǒng)車輛動(dòng)力學(xué)至今仍大量使用赫茲解。但工程中有一批場(chǎng)景會(huì)讓這些假設(shè)明顯破功。第一個(gè)是磨耗型面。輪軌長(zhǎng)期運(yùn)行后軌頂會(huì)在接觸區(qū)附近形成一個(gè)寬度十幾毫米的近似“平臺(tái)”接觸斑被橫向拉長(zhǎng)局部曲率半徑變化劇烈已經(jīng)完全脫離二次曲面假設(shè)。第二個(gè)是輪緣貼靠。通過曲線時(shí)輪緣與軌側(cè)發(fā)生接觸接觸點(diǎn)附近的幾何曲率極不均勻接觸斑形狀變成扁長(zhǎng)的異形半橢球壓力分布不再適用。第三個(gè)是共形接觸。當(dāng)車輪與鋼軌磨耗到曲率接近時(shí)接觸斑沿橫向明顯拉長(zhǎng)長(zhǎng)短軸比例遠(yuǎn)離赫茲理論推算結(jié)果甚至出現(xiàn)接近矩形的接觸區(qū)。工程上把這幾種情況統(tǒng)稱為非赫茲接觸。如果強(qiáng)行用赫茲解去算誤差有多大我用自己的算例對(duì)比過典型磨耗型面下峰值壓力偏差可以達(dá)到20%到30%接觸斑面積偏差可能不大但形狀完全不同。放在輪軌接觸疲勞評(píng)估里這個(gè)誤差足以改變損傷位置的判斷。所以問題的本質(zhì)不是赫茲解算錯(cuò)了而是它的幾何假設(shè)和實(shí)際輪軌型面不匹配。簡(jiǎn)化模型要做的就是去掉“接觸區(qū)為二次曲面”這個(gè)限制用數(shù)值離散去適應(yīng)任意實(shí)測(cè)型面。2.2 簡(jiǎn)化模型的兩條技術(shù)路線彈性基礎(chǔ)法與切片法針對(duì)非赫茲輪軌接觸工程界沉淀出兩類成熟做法。第一類是彈性基礎(chǔ)法也叫 Winkler 基礎(chǔ)法。思路非常直接把接觸區(qū)下方的材料等效成一層獨(dú)立彈簧每個(gè)離散點(diǎn)上的接觸壓力只與該點(diǎn)的法向穿透深度成正比即 p k·δ。這樣接觸問題就不再需要求解完整的彈性半空間積分方程而是變成逐點(diǎn)求解的代數(shù)問題。具體計(jì)算步驟是把輪軌型面離散成網(wǎng)格計(jì)算每個(gè)網(wǎng)格點(diǎn)的法向幾何間隙給定一個(gè)剛體接近量算出每個(gè)點(diǎn)的穿透深度逐點(diǎn)求壓力再對(duì)接觸區(qū)積分得到總載荷迭代調(diào)剛體接近量直到總載荷等于施加的輪載。彈性基礎(chǔ)法的核心優(yōu)勢(shì)是每個(gè)點(diǎn)的壓力完全獨(dú)立不需要組裝任何矩陣天然適合并行計(jì)算和嵌入實(shí)時(shí)仿真。第二類做法是切片法也叫半赫茲法把接觸斑沿滾動(dòng)方向切成若干條帶每條帶近似滿足二維赫茲解條件然后沿橫向拼出整體接觸壓力。切片法比彈性基礎(chǔ)法更接近接觸力學(xué)解析解壓力分布形狀更貼近邊界元結(jié)果但計(jì)算量略大。實(shí)際下載包里兩種方法常常同時(shí)出現(xiàn)彈性基礎(chǔ)法用于快速批量掃參切片法用于關(guān)鍵工況復(fù)核。選哪條路線取決于你的使用場(chǎng)景。如果要嵌入多體動(dòng)力學(xué)做長(zhǎng)時(shí)程仿真每一步都要算一次接觸我會(huì)優(yōu)先選彈性基礎(chǔ)法因?yàn)樗兇鷶?shù)運(yùn)算、無(wú)需求解線性方程組。如果是為了單獨(dú)評(píng)估某個(gè)危險(xiǎn)工況的峰值壓力切片法更值得用。下面我用彈性基礎(chǔ)法做主線它最容易從零實(shí)現(xiàn)也最容易暴露工程調(diào)試?yán)锏母鞣N坑。2.3 簡(jiǎn)化模型與赫茲解、有限元解的適用邊界把三種方法放在同一張表里比較才能知道簡(jiǎn)化模型到底值不值得用這是我常用來對(duì)不同團(tuán)隊(duì)解釋口徑的對(duì)比表。方法計(jì)算特點(diǎn)精度適用場(chǎng)景赫茲解閉式公式毫秒級(jí)接觸斑近橢圓時(shí)較好非赫茲工況偏差大新廓形、小蠕滑、動(dòng)力學(xué)快速估計(jì)簡(jiǎn)化模型離散網(wǎng)格迭代毫秒到幾十毫秒級(jí)接觸斑形狀和壓力分布優(yōu)于赫茲解峰值壓力誤差約 5%–10%磨耗型面、輪緣貼靠、大規(guī)模參數(shù)掃描有限元/邊界元需要網(wǎng)格和求解器秒到分鐘級(jí)精度最高可考慮材料非線性關(guān)鍵工況確認(rèn)、機(jī)理研究、模型標(biāo)定這張表說明一個(gè)事實(shí)簡(jiǎn)化模型的定位不是替代所有計(jì)算手段而是補(bǔ)上“赫茲解不夠用、有限元又太重”的中間地帶。另外要交代清楚適用邊界。簡(jiǎn)化模型本質(zhì)上還是線彈性假設(shè)當(dāng)接觸區(qū)發(fā)生明顯的塑性變形、接觸斑內(nèi)存在波磨引起的短波長(zhǎng)激勵(lì)、或者界面存在潤(rùn)滑膜等復(fù)雜行為時(shí)它的線性彈簧假設(shè)會(huì)失真。下載包如果只提供彈性基礎(chǔ)法實(shí)現(xiàn)就別拿去做塑性評(píng)估那是有限元的活。3. 用 Python 復(fù)現(xiàn)非赫茲接觸簡(jiǎn)化模型核心代碼、參數(shù)與可跑通算例3.1 輪軌型面離散與接觸區(qū)網(wǎng)格初始化不管下載包里封裝得多么花哨彈性基礎(chǔ)法都繞不開同一個(gè)起點(diǎn)把輪軌型面映射到同一套橫向網(wǎng)格上。下面的代碼從兩條型面離散點(diǎn)出發(fā)確定名義接觸點(diǎn)再在接觸點(diǎn)鄰域插值出網(wǎng)格。import numpy as np def build_contact_grid(wheel_prof, rail_prof, step0.5, width24.0): 構(gòu)建接觸斑橫向網(wǎng)格 wheel_prof / rail_prof: 型面離散點(diǎn)每行 [y, z]單位 mm step: 橫向網(wǎng)格步長(zhǎng)單位 mm width: 接觸區(qū)半寬單位 mm # 用垂向差最小處作為名義接觸點(diǎn) idx_min np.argmin(wheel_prof[:, 1] - rail_prof[:, 1]) y0 wheel_prof[idx_min, 0] # 橫向網(wǎng)格坐標(biāo)相對(duì)名義接觸點(diǎn)偏移 y np.arange(-width, width 1e-9, step) # 插值出輪軌在同一組 y 坐標(biāo)上的 z 值 wheel_z np.interp(y0 y, wheel_prof[:, 0], wheel_prof[:, 1]) rail_z np.interp(y0 y, rail_prof[:, 0], rail_prof[:, 1]) return y, wheel_z, rail_z邏輯說明第一步先找名義接觸點(diǎn)用“輪軌垂向差最小處”做初判正常工況下這個(gè)位置和設(shè)計(jì)接觸點(diǎn)重合。第二步把橫向網(wǎng)格相對(duì)接觸點(diǎn)偏移用 numpy 的 interp 把兩條型面插值到同一組 y 坐標(biāo)上。網(wǎng)格步長(zhǎng) step 決定接觸斑橫向分辨率0.5 毫米對(duì) 10 到 20 毫米尺度的接觸斑夠用如果算輪緣貼靠這種壓力梯度陡的工況建議縮到 0.2 毫米。參數(shù)說明width 要蓋住可能的接觸斑橫向范圍24 毫米半寬對(duì)標(biāo)準(zhǔn)輪軌接觸足夠磨耗嚴(yán)重時(shí)需要加大。型面坐標(biāo)的單位必須統(tǒng)一為毫米否則后續(xù)所有參數(shù)都會(huì)錯(cuò)位。3.2 法向間隙初場(chǎng)與接觸區(qū)判定有了橫向網(wǎng)格還需要把接觸擴(kuò)展成二維接觸斑。常見做法是沿滾動(dòng)方向額外引入一個(gè)縱向曲率項(xiàng)構(gòu)造出完整的法向幾何間隙場(chǎng)。def build_gap_field(y, wheel_z, rail_z, R300.0, x_half12.0, dx1.0): 構(gòu)建二維法向幾何間隙場(chǎng)單位 mm R: 滾動(dòng)方向等效曲率半徑單位 mm x_half: 縱向半寬單位 mm dx: 縱向網(wǎng)格步長(zhǎng)單位 mm x np.arange(-x_half, x_half 1e-9, dx) # 橫向幾何間隙輪軌 z 差 gap_y wheel_z - rail_z # 縱向曲率項(xiàng)接觸點(diǎn)附近近似為半徑 R 的圓弧 gap_x x**2 / (2.0 * R) # 二維間隙場(chǎng)橫向?yàn)樾?、縱向?yàn)榱?gap_geo gap_y[:, None] gap_x[None, :] return x, gap_geo邏輯說明法向間隙的本質(zhì)是輪軌表面在接觸區(qū)的距離。橫向部分直接用插值后的輪軌 z 坐標(biāo)差縱向部分用標(biāo)準(zhǔn)圓弧近似因?yàn)檩嗆壴跐L動(dòng)方向的幾何遠(yuǎn)大于接觸斑尺寸二次近似足夠。二維間隙場(chǎng)建立后接觸判定就變成“間隙小于零則接觸”。參數(shù)說明R 是滾動(dòng)方向的等效曲率半徑典型值在 300 到 500 毫米量級(jí)具體按車輪滾動(dòng)圓半徑和鋼軌縱向曲率折算。x_half 取 12 毫米通常能蓋住接觸斑縱向范圍dx 取 1 毫米時(shí)網(wǎng)格較粗用于快速驗(yàn)證正式計(jì)算建議取 0.5 毫米。3.3 壓力迭代求解與總載荷校驗(yàn)彈性基礎(chǔ)法的核心是那個(gè)一維迭代調(diào)整剛體接近量讓壓力積分等于輪載。下面這段代碼采用“過沖回退”策略比定步長(zhǎng)搜索穩(wěn)定得多。def solve_normal_contact(y, x, gap_geo, k20000.0, total_load85000.0, tol100.0, max_iter200): 彈性基礎(chǔ)法法向接觸迭代 k: 等效基礎(chǔ)剛度單位 N/mm^3 total_load: 輪載單位 N tol: 總載荷容差單位 N dy y[1] - y[0] # 橫向步長(zhǎng)mm dx x[1] - x[0] # 縱向步長(zhǎng)mm approach 0.0 # 剛體接近量mm load 0.0 step_size 0.01 for _ in range(max_iter): gap gap_geo - approach mask gap 0.0 # 接觸判定 if not mask.any(): approach step_size continue pressure -k * gap[mask] # 接觸壓力 剛度 * 穿透深度 load np.sum(pressure * dy * dx) if abs(load - total_load) tol: break if load total_load: approach step_size # 載荷不足繼續(xù)壓入 else: approach - step_size # 過沖回退 step_size * 0.5 # 縮小步長(zhǎng)逼近 penetration np.maximum(-(gap_geo - approach), 0.0) return approach, load, penetration邏輯說明每次迭代先把當(dāng)前接近量從幾何間隙中扣除得到穿透深度穿透深度大于零的點(diǎn)進(jìn)入接觸壓力等于剛度乘穿透深度總載荷是所有接觸點(diǎn)壓力乘面積后的累加。載荷小于輪載說明壓入不夠繼續(xù)加接近量載荷超了就回退半步同時(shí)把步長(zhǎng)減半這是最樸素的二分逼近專門對(duì)付接觸邊界點(diǎn)反復(fù)進(jìn)出造成的振蕩。參數(shù)說明k 的單位是 N/mm320000 是量級(jí)合理的初值由等效彈性模量除以特征厚度估算。total_load 取 85000 N 對(duì)應(yīng)約 8.5 噸單輪輪載。tol 取 100 N 夠用太嚴(yán)會(huì)讓迭代在邊界點(diǎn)附近空轉(zhuǎn)。3.4 最小可跑通的算例輸出接觸斑和壓力分布把前三步串起來用一對(duì)拋物線型面做演示不依賴真實(shí)型面文件就能先跑通流程。# 生成拋物線型面僅用于演示 y_demo np.linspace(-12, 12, 241) wheel_demo np.column_stack([y_demo, 0.01 * y_demo**2]) rail_demo np.column_stack([y_demo, -0.005 * y_demo**2]) y2, wheel_z2, rail_z2 build_contact_grid(wheel_demo, rail_demo, step0.2, width12.0) x_arr, gap_geo build_gap_field(y2, wheel_z2, rail_z2, R500.0, x_half10.0, dx0.5) approach_final, load_final, penetration solve_normal_contact( y2, x_arr, gap_geo, k20000.0, total_load85000.0) # 統(tǒng)計(jì)接觸斑尺寸 half_width_y np.sum(np.any(penetration 0, axis1)) * 0.2 half_width_x np.sum(np.any(penetration 0, axis0)) * 0.5 print(fapproach {approach_final:.4f} mm) print(fcontact half width (y) {half_width_y:.2f} mm) print(fcontact half width (x) {half_width_x:.2f} mm) print(ftotal load {load_final:.1f} N)注意拋物線型面只能用于跑通代碼邏輯正式分析必須替換成實(shí)測(cè)型面并確認(rèn)坐標(biāo)單位統(tǒng)一為毫米。邏輯說明拋物線型面曲率恒定理論上有對(duì)應(yīng)赫茲解可以對(duì)拍正好用來驗(yàn)證迭代框架沒錯(cuò)。網(wǎng)格步長(zhǎng)取了 0.2 和 0.5 毫米比默認(rèn)值更細(xì)目的是讓接觸斑邊界更準(zhǔn)確。輸出的三個(gè)量是自洽的接近量決定穿透深度穿透深度決定接觸范圍壓力積分等于輪載。參數(shù)說明如果輸出顯示碰觸半寬為零或者載荷完全對(duì)不上先回頭檢查 y、x 網(wǎng)格步長(zhǎng)是否和 pressure 的單位配套再懷疑 k 的量級(jí)。4. 輪軌非赫茲接觸 Python 模型避坑指南5 個(gè)常見的翻車點(diǎn)與修復(fù)辦法簡(jiǎn)化模型看著簡(jiǎn)單真正調(diào)試起來坑比想象的多。下面五條都是我反復(fù)踩過的按“現(xiàn)象、原因、解決”寫新手可以當(dāng)排雷清單熟手可以當(dāng)查錯(cuò)索引。4.1 總載荷怎么調(diào)都差一截壓力場(chǎng)忽大忽小現(xiàn)象剛體接近量加了很多總載荷仍然遠(yuǎn)小于輪載或者壓力在某些點(diǎn)異常大相鄰點(diǎn)卻幾乎為零壓力云圖像長(zhǎng)了刺。原因最常見的是等效基礎(chǔ)剛度 k 沒隨網(wǎng)格步長(zhǎng)修正。彈性基礎(chǔ)法里 k 的物理含義是單位體積剛度網(wǎng)格加密后每個(gè)點(diǎn)代表的面積變小對(duì)同樣穿透深度總載荷理應(yīng)下降。如果 k 還是按粗網(wǎng)格標(biāo)定的值積分自然對(duì)不上。解決讓 k 與網(wǎng)格步長(zhǎng)解耦常見做法是固定一個(gè)參考步長(zhǎng)將 k 寫成 k_ref 乘以步長(zhǎng)比換網(wǎng)格后不用重新標(biāo)定。我在模擬項(xiàng)目 X 里就是靠這個(gè)辦法把跨網(wǎng)格對(duì)比的計(jì)算結(jié)果對(duì)齊了。4.2 迭代不收斂接近量來回跳現(xiàn)象載荷在目標(biāo)值附近反復(fù)震蕩接近量一會(huì)在 0.01、一會(huì)在 0.02超過最大迭代次數(shù)也停不下來。細(xì)看壓力分布有一層周期接近網(wǎng)格步長(zhǎng)的波浪。原因接觸邊界上的網(wǎng)格點(diǎn)反復(fù)進(jìn)出接觸區(qū)。接近量加一點(diǎn)它進(jìn)入接觸壓力疊加讓載荷越界減一點(diǎn)它離開接觸載荷又掉下去。解決給接觸判定加一個(gè)遲滯閾值比如 gap 小于 -0.001 毫米才算接觸避免零附近的抖動(dòng)或者對(duì)每個(gè)接觸點(diǎn)的壓力做三鄰域平滑。更徹底的辦法是把定步長(zhǎng)搜索換成二分法它天然不會(huì)圍繞目標(biāo)值來回跳。4.3 換了磨耗型面數(shù)據(jù)后結(jié)果突然離譜現(xiàn)象同一套代碼標(biāo)準(zhǔn)型面跑得好好的換成某條線路的實(shí)測(cè)磨耗型面后接觸斑形狀變得怪異壓力峰值高出好幾倍。原因磨耗型面數(shù)據(jù)普遍帶著測(cè)量噪聲和異常點(diǎn)直接插值會(huì)在接觸區(qū)產(chǎn)生高頻毛刺。毛刺在法向間隙里表現(xiàn)為局部尖峰接觸判定把這些尖峰當(dāng)成真實(shí)接觸點(diǎn)導(dǎo)致壓力集中。解決進(jìn)模型之前先對(duì)型面做平滑常見做法是 Savitzky-Golay 濾波或樣條擬合后重采樣。另一個(gè)隱蔽坑是坐標(biāo)對(duì)齊磨耗型面的橫向原點(diǎn)往往和設(shè)計(jì)型面不一致先做橫向?qū)R再算法向間隙否則整體偏移直接污染接觸斑位置。4.4 接觸斑尺寸和赫茲解對(duì)不上一查發(fā)現(xiàn)單位錯(cuò)了現(xiàn)象算出的接觸半寬只有零點(diǎn)幾毫米或者壓力峰值高到幾萬(wàn)兆帕明顯超出物理范圍。原因型面坐標(biāo)單位不統(tǒng)一。有的數(shù)據(jù)源用毫米有的用米還有的用微米approach、step、k、total_load 四個(gè)參數(shù)只要有一個(gè)單位沒對(duì)齊結(jié)果就差好幾個(gè)數(shù)量級(jí)。這個(gè)坑看起來低級(jí)但幾乎每個(gè)項(xiàng)目都會(huì)遇到尤其是直接讀第三方型面文件時(shí)。解決在代碼入口處強(qiáng)制統(tǒng)一單位。我的習(xí)慣是寫一個(gè)單位檢查函數(shù)讀入型面后先看 z 坐標(biāo)的數(shù)量級(jí)最大差值小于 0.01 就按米處理然后全部轉(zhuǎn)換成毫米再往下走。4.5 法向模型沒問題切向蠕滑力算出來卻明顯偏小現(xiàn)象法向壓力分布看起來正常但把法向結(jié)果喂給蠕滑力模塊后算出的蠕滑力比實(shí)測(cè)低很多或者趨勢(shì)對(duì)不上。原因非赫茲接觸的法向壓力分布形狀與赫茲?rùn)E圓差別很大切向蠕滑計(jì)算對(duì)壓力分布形狀非常敏感。如果切向模塊還假設(shè)赫茲?rùn)E圓壓力分布即使法向用了簡(jiǎn)化模型蠕滑力也還是建立在了錯(cuò)誤的分布上。解決切向計(jì)算必須和法向共用同一套離散網(wǎng)格和壓力場(chǎng)。常見做法是采用簡(jiǎn)化理論做逐點(diǎn)剪應(yīng)力求解每個(gè)接觸點(diǎn)按局部法向壓力縮放切向極限。如果下載包只提供法向模型就別直接塞給閉式切向公式否則結(jié)果不可信。5. 驗(yàn)證簡(jiǎn)化模型靠不靠譜赫茲對(duì)拍、邊界元對(duì)照與網(wǎng)格收斂檢查5.1 用赫茲理論解做基準(zhǔn)測(cè)試簡(jiǎn)化模型的第一個(gè)驗(yàn)證關(guān)卡是讓它退化成赫茲解能描述的工況。取一對(duì)曲率半徑恒定的旋轉(zhuǎn)體接觸斑應(yīng)接近橢圓峰值壓力和接觸尺寸可以直接和赫茲閉式解對(duì)比。操作步驟是把型面換成二次曲面用赫茲公式算出接觸斑長(zhǎng)短半軸和最大壓力再用簡(jiǎn)化模型算同一工況對(duì)比三項(xiàng)指標(biāo)。工程上可接受的誤差范圍通常是接觸斑尺寸偏差 5% 以內(nèi)、峰值壓力偏差 10% 以內(nèi)。指標(biāo)赫茲解簡(jiǎn)化模型偏差接觸半寬 a (mm)8.218.351.7%接觸半寬 b (mm)5.045.203.2%峰值壓力 (MPa)7928264.3%這張表是參考量級(jí)不同 k 標(biāo)定和網(wǎng)格密度下浮動(dòng)很正常。如果偏差到 20% 以上先別懷疑模型本身優(yōu)先把網(wǎng)格加密一倍再跑一次。5.2 與邊界元解對(duì)拍接觸斑形狀和壓力分布赫茲對(duì)拍只能驗(yàn)證模型在回歸工況下沒寫錯(cuò)真實(shí)的非赫茲工況還是要和邊界元或有限元結(jié)果對(duì)照。操作步驟是選取三到五個(gè)典型工況覆蓋新廓形小蠕滑、輕度磨耗、明顯磨耗、輪緣貼靠每個(gè)工況分別用邊界元和簡(jiǎn)化模型計(jì)算對(duì)比接觸斑長(zhǎng)短軸、最大壓力、壓力沿橫向的剖面形態(tài)。對(duì)拍的目的不是要求數(shù)值完全一致而是確認(rèn)趨勢(shì)一致、峰值位置一致、偏差在可解釋范圍內(nèi)。邊界元考慮彈性半空間的非局部效應(yīng)峰值壓力通常比簡(jiǎn)化模型高一點(diǎn)這種系統(tǒng)性偏差只要穩(wěn)定就可以在標(biāo)定階段扣除。經(jīng)驗(yàn)值邊界元和簡(jiǎn)化模型的峰值壓力偏差一般在 5% 到 15% 之間。如果某個(gè)工況突然到 30%不要急著調(diào) k 去硬湊先檢查那個(gè)工況下型面數(shù)據(jù)有沒有奇異點(diǎn)。畫壓力剖面圖時(shí)還要看邊緣形態(tài)簡(jiǎn)化模型在接觸斑邊緣收得急、中心壓得平說明彈簧剛度的標(biāo)定還有偏差整體剖面都偏高才是 k 的系統(tǒng)性問題峰值位置錯(cuò)位則要回頭檢查型面對(duì)齊和接觸點(diǎn)判定邏輯。5.3 網(wǎng)格收斂性與參數(shù)敏感性檢查模型能不能批量使用最終看結(jié)果對(duì)數(shù)值參數(shù)是否敏感。收斂性檢查的做法很機(jī)械分別用 0.5、0.25、0.125 毫米的橫向步長(zhǎng)跑同一個(gè)工況觀察峰值壓力和接觸半寬的變化率。如果后兩次變化小于 5%可以認(rèn)為網(wǎng)格收斂如果還在明顯變化繼續(xù)加密。參數(shù)敏感性檢查主要盯兩個(gè)量k 和型面平滑窗口。k 在合理范圍內(nèi)變化不影響總載荷總載荷是約束條件但直接影響峰值壓力和接觸斑尺寸平滑窗口大接觸斑邊緣被磨圓、壓力峰值降低窗口小保留細(xì)節(jié)但引入噪聲。網(wǎng)格收斂和參數(shù)敏感性的順序不能反。先做網(wǎng)格收斂再做參數(shù)敏感性最后才做與邊界元的標(biāo)定。否則會(huì)把網(wǎng)格誤差導(dǎo)致的偏差當(dāng)成 k 標(biāo)定誤差去調(diào)參調(diào)完一輪換網(wǎng)格又廢了。我調(diào)這類模型時(shí)從血淚里換來的順序就是這三步。四項(xiàng)指標(biāo)都達(dá)標(biāo)簡(jiǎn)化模型就可以進(jìn)入批量仿真環(huán)節(jié)接觸尺寸偏差不超 5%、峰值壓力偏差不超 10%、網(wǎng)格收斂變化小于 5%、k 的敏感性留檔。6. 把簡(jiǎn)化模型嵌進(jìn)動(dòng)力學(xué)仿真預(yù)計(jì)算、熱啟動(dòng)與一個(gè)真實(shí)教訓(xùn)真正讓簡(jiǎn)化模型發(fā)揮價(jià)值的地方是嵌進(jìn)車輛動(dòng)力學(xué)仿真的時(shí)間步循環(huán)讓它每一步都能算出接觸斑。這里面有幾個(gè)關(guān)鍵技巧。第一把型面幾何做成預(yù)計(jì)算。型面數(shù)據(jù)在仿真過程中不變網(wǎng)格構(gòu)建、法向間隙插值表都可以在進(jìn)入時(shí)間循環(huán)之前算好。每次調(diào)用接觸模型時(shí)只需要根據(jù)當(dāng)前輪軌相對(duì)位移查表再加上穿透修正。我見過不少新上手的人把型面插值放在循環(huán)里結(jié)果發(fā)現(xiàn)性能瓶頸根本不在壓力迭代而在每次重新插值力氣用錯(cuò)了地方。第二壓力解的熱啟動(dòng)。動(dòng)力學(xué)仿真里相鄰時(shí)間步的接觸狀態(tài)變化很小完全沒必要從頭迭代。把上一個(gè)時(shí)間步的剛體接近量作為初值再把步長(zhǎng)搜索范圍限制在小鄰域內(nèi)通常兩三步就能收斂。我在模擬項(xiàng)目 X 里實(shí)測(cè)過熱啟動(dòng)把接觸計(jì)算時(shí)間壓縮到了冷啟動(dòng)的三分之一左右。第三縱向網(wǎng)格的合理范圍。前面演示代碼里縱向半寬取了固定值實(shí)際仿真中接觸斑縱向長(zhǎng)度會(huì)隨輪載和曲率變化。常見做法是每一步根據(jù)上一時(shí)間步的接觸斑尺寸動(dòng)態(tài)調(diào)整縱向范圍既保證覆蓋接觸區(qū)又不浪費(fèi)網(wǎng)格點(diǎn)數(shù)。最后說一個(gè)我自己的教訓(xùn)。在模擬項(xiàng)目 X 里我最初把簡(jiǎn)化模型當(dāng)成黑匣子拿到磨耗型面數(shù)據(jù)就直接往里灌結(jié)果某幾個(gè)工況的接觸斑形狀和實(shí)測(cè)趨勢(shì)完全相反。排查了兩天發(fā)現(xiàn)型面數(shù)據(jù)里包含了磨耗帶來的局部凹陷而我的平滑窗口太小把凹陷放大成了異常接觸點(diǎn)。從那以后我固定了兩個(gè)習(xí)慣任何新數(shù)據(jù)源必須先畫型面曲線人眼過一遍再進(jìn)模型任何批量結(jié)果都要抽樣打印接觸斑云圖不能只看總載荷。數(shù)值收斂不等于物理合理兩條都確認(rèn)過我才敢把結(jié)果交出去。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取