混凝土徐變柔度函數逼近的混合智能優化方法

中圖分類號：TU528.1 文獻標志碼：A

本文引用格式：，陳夢成，謝力，等.混凝土徐變柔度函數逼近的混合智能優化方法[J].華東交通大學學報，2025，42（2）：46-53.

A Hybrid Intelligent Optimization Method for Approximating the Creep Compliance Function of Concrete

Ge Jing2， Chen Mengcheng12， Xie Li² ，Huang Hong12， Xu Kaicheng12，Fang Wei1，2 （1.State KeyLaboratoryof Performance Monitoring and Protecting ofRail Transit Infrastructure，East China Jiaotong University，Nanchang 33o013， China; 2.School of Civil Engineering amp; Architecture， East China Jiaotong University，Nanchang330013，China）

Abstract：The creep compliance function is usually expressed bya Prony series model when using finite element codes to calculate the creep effects of concrete structures.To address the physically meaningless phenomenon of negative，oscilatory and unstable parameters when using traditional single methods （such as the conjugate gradient method and the LM algorithm） to fit Prony series models to experimental creep function data，this paper proposesanovel hybrid intellgent optimization method for Prony series parameter identification.The method establishes an objective optimization function with penalty terms and combines a simulated annealing-genetic hybrid inteligent algorithm with nonlinear programming methods to constrain parameter identification within physically meaningful ranges.Subsequently，a simple and practical calculation formulation for the Prony series was proposed by parameter analysis. It was shown by numerical examples that the proposed method not only effectively eliminates the drawbacks ofconventional approaches，butalso gives the identification with relative erors below

3% .The Prony expression of creep function can be directly applied to the development of finite element software program in calculating the creep effect of concrete structures.

Key words： concrete material; viscoelasticity; creep function; prony series; hybrid intelligent optimization Citation format： GE J， CHEN M C， XIE L，etal. A hybrid inteligent optimization method for approximating the creep compliance function of concrete[J].Journal ofEast China Jiaotong University，2025，42（2）： 46-53.

徐變是混凝土材料的基本特性，其在現代大型混凝土工程結構中的時變效應表現尤為顯著[]，對結構的力學性能與安全性有重要影響。因此精準預測混凝土結構徐變是分析結構時變應力-應變響應的關鍵。

根據線性黏彈性理論中Boltzmann疊加原理，混凝土時變應力-應變關系可以表示為以下第一類典型Volterra積分方程[4]

ε（t）-ε^?（t）=∈t₀^tJ（t，τ₀）dσ（τ₀）

式中： χ_t 為混凝土齡期，d; τ₀ 為加載齡期，d; σ 為混凝土應力， MPa ε 為混凝土應變; J（t，τ₀） 為混凝土徐變函數或徐變柔量； ε⁰ 為與應力無關的非彈性收縮應變或溫度應變。對式（1）進行積分計算時，需要將徐變函數試驗值擬合成線性黏彈性數學模型。

求解混凝土結構徐變效應時，為提升數值積分計算效率，Bazant4使用了Prony級數表達形式，這一形式現也被廣泛應用于Abaqus和Comsol等商用有限元軟件[5-8]。求解之前，精準有效識別Prony級數模型參數是關鍵所在。目前國內外有很多確定Prony級數模型參數的研究報道，方法主要有預平順法，最小二乘法，共軛梯度法，連續延遲譜法和優化法等[46.9-14]。然而，Prony級數模型在徐變函數表征中存在若干理論局限性。首先，其指數型基函數因固有頻帶局限性，對試驗數據的離散特征表現出較高敏感性。而且，模型參數具有非唯一性特征，不同參數組合均可實現對徐變試驗值的擬合。盡管徐變彈性模量在整體上為正且單調遞增，但由于多解性的影響，其在某些時段的具體參數值可能出現負值或振蕩現象，這與其物理意義相悖。另外，傳統參數識別方法（包括共軛梯度法和LM法）存在對初始值差異較敏感的問題。

本文提出了一種新的Prony級數模型參數識別的智能優化方法。其基本原理是以最小二乘準則為數學基礎，建立目標優化函數，使徐變函數試驗數據值與Prony級數理論模型計算值之差平方和極小化。建立了模擬退火-遺傳混合優化策略（SA-GA）與非線性約束規劃協同的Prony級數模型參數識別機制。相較于傳統共軛梯度法和LM法，本文方法能夠有效克服延遲強度識別結果出現負數和不穩定的問題，而且識別收斂快、精度高，可用于現有有限元商用軟件二次開發，進而計算分析混凝土結構徐變效應。

1基本概念和數學公式

1.1線性黏彈性本構關系與Prony級數

當應變為小應變且不出現逆應變和循環應變時，混凝土材料可以近似為線性黏彈性材料，其力學響應遵循疊加原理。基于此理論框架，單軸時變混凝土體系的本構描述可通過式（1）的積分方程實現。這個積分方程具有Stieltjes積分含義，允許用有限個階段應力描述分段連續應力過程。當混凝土應力 σ（τ） 全過程已知，應變可以直接從方程（1）計算得到。然而，方程（1）不大適合一般數值計算，因為時間通常離散為有限個增量 Δt_n ，最終形成一個有關 Δσ_n 和 Δε_n 的準靜態結構分析本構方程。每計算一步應力 Δσ_n ，對前 n-1 步 Δσ 進行求和、存儲，需要大量的計算機存儲空間和CPU運行時間。

為了解決上面問題，Bazant針對老化混凝土，提出了用Prony級數表達式擬合模擬方程（1）中的徐變函數 ，具體表達式如下

式中： ， E_μ^-1（τ₀） 和 τ_μ 分別為玻璃態柔量、延遲強度和延遲時間；玻璃態柔量為瞬態柔量，即： 。表達式（2）與如圖1所示的廣義Kelvin鏈模型的力學模型相關聯

Prony級數的重要優勢是能夠描述黏彈性材料寬帶現象，且因其基函數為指數函數可有效提升計算效率。將式（2）代入方程（1），得到基于廣義Kel-vin鏈模型的應力-應變本構方程，即

其中

式（4）可以看作為式（3）的內變量，描繪Kelvin元件當前時刻之前的應變變化過程。

1.2 模型參數識別優化目標函數

假定徐變函數試驗測試點值或預平滑試驗曲線上點值用 J^Test（t_r，τ₀） ， r=1，2，…，N_t 表示，那么根據最小二乘準則，在加載齡期 τ_0m，m=1，2，…，N_τ0 選定的情況下，與試驗值等效的理論曲線模型Prony級數表達式（2）中的系數 E_μ^-1（τ₀） （注意： E₀^-1（τ₀） 沒有計入）， μ=1，2，…，N_μ ，可以使用以下最小化標量函數表達式進行識別優化

式中第2項＼～第4項代表罰函數，起到光滑數據的作用，目的是使 E₁^-1（τ_0m） ， E₂^-1（τ_0m） ，…， 之間差平方最小，同時降低 E_μ^-1 對擬合數據離散性的敏感度。 ω_?1 ， ω₂ 和 ω₃ 取值不應該過大，應該盡量小，一般小于1，可以預先給定（取決于計算經驗，本文取ω=@=ω=0.0001），也可作為參數優化識別。

2模型參數智能優化識別基本思想

2.1基礎工作準備

2.1.1 延遲時間均勻化離散點選擇

在固定加載齡期 τ_0m 條件下，基于Prony級數J^Prony（t，τ₀） 識別模型參數（材料參數） E_μ^-1（τ₀） ，或者擬合給定的試驗數據 時，存在數學層面上的困難。其根本原因是Prony級數系數 E_μ^-1（τ₀） 和 τ_μ 之間的依存關系不穩定，完全依賴于給定的擬合數據。若直接從試驗數據確定延遲時間 τ_μ ，將因 τ_μ 的多解性（不同的 τ_μ 解均可等效逼近試驗數據曲線）導致Prony級數模型參數求解方程呈現高度病態。理論分析表明，混凝土徐變函數延遲時間譜完全可以用離散化方法進行表征，而這些離散值對應著不同的時間細分。因此，需預先構建松弛時間譜的優選策略。實際上，在時間對數尺度上均勻選擇τ₁，…，τ_Nμ 值最佳。根據Bazant4計算經驗，在延遲時間對數軸上按以下方式均勻選定 τ_μ 值，足以保證參數識別的計算精度

τ_μ=τ₁?10^μ-1，μ=1，2，…，N_μ

其中， τ₁ 取值應該和時間對數 log（t-τ₀） 軸上混凝土徐變函數曲線開始出現上升的時間點或者研究所需的延遲時間譜下限保持一致，而 τ_Nμ 取值則應與時間對數 log（t-τ₀） 軸上徐變函數曲線開始出現下降的時間點或者研究所需的延遲時間譜上限保持一致。

2.1.2計算時間均勻化離散點確定

同樣，在計算時間對數 log（t-τ₀） 軸上對時間進行均勻離散，為 t_r，r=1，2，…，N_t ，每10年中有4個值為最佳選擇，即 ，為提高計算結果精度，按以下形式取值

需要說明的是，試驗測試值通常不可能按式（7）均勻分布，必須內插使其在 （t-τ₀） 內均勻分布。

2.2 基本方程

將式（6）和式（7）代人式（5），Prony級數模型參數 E₁^-1（τ_0m） ， E₂^-1（τ_0m） ，…， E_?Nμ^-1（τ_0m） 的識別問題轉化為求解下面標準優化問題

式中： x 為設計變量； 為任意選定值，代表設計變量 E_μ^-1（τ_0m） 的上界。

2.3 MATLAB程序實現

為解決Prony級數模型參數智能優化識別問題，基于式（8），給出了智能優化過程Matlab軟件程序設計計算結構框架流程圖（圖2）。首先依據

2.1.2中給出的混凝土持荷時間均勻化理論，對徐變函數試驗數據在時間離散節點進行采樣；其次根據模擬退火遺傳算法SAGA程序輸入參數要求輸人采樣點值，以徐變函數采樣點數據值與徐變函數Prony級數理論模型計算值之差的平方的和最小為優化目標函數，計算得到局部最優解；然后，在同一目標下，以該局部最優解作為初始值，輸人有約束非線性規劃方法設計程序進行計算，得到全局最優解；最后，輸出求解結果。在參數識別優化過程中，使用了SAGA程序自編Matlab程序，結合Matlab軟件中有約束非線性規劃算法的函數“fmincon\"命令。

3應用實例與計算分析

3.1徐變函數試驗數據

為闡述本文方法，徐變函數試驗數據采用Bazant等提出的混凝土徐變收縮預測B3模型[15進行計算。具體表達式為

J^B3（t，τ₀）=q₁+C₀（t，τ₀）+C_d（t，τ₀，t_sh）

式中： t 為混凝土的齡期，d; τ₀ 為加載齡期，d， t_sh 為混凝土干燥開始的齡期，d; q₁ 為單位應力下產生的瞬時應變； 為基本徐變度； 為干燥徐變度，相關項具體計算式及徐變函數試驗數據計算中采用的混凝土參數參考文獻[15]。 q₁ ， 的表達式如下

式（11）基本徐變度的Prony級數模型，可以通過參數智能優化識別寫成如下形式

其中，

3.2 智能優化識別過程輸入參數

智能優化識別中常見的輸入參數見表1。在數值識別過程中，式（12）中Prony級數的模型項數 N_μ 選為7項，各項對應的延遲時間 τ_μ 按式（6）方式選取，見表2；計算時間選定方式按式（7）進行。

3.3Prony級數模型參數識別結果分析討論

3.3.1不同識別方法結果比較

分別使用共軛梯度法、LM法、無罰項和有罰項智能識別優化法對Prony級數模型參數延遲強度進行擬合，識別對象為加載齡期為28d時B3模型基本徐變度中延遲強度第3部分 E_c，μ^-1（28） （見式（13）），具體模擬結果見表3。從相關系數來看，4種方法擬合精度都不錯，但延遲時間較小 （0lt;τ_μ?1） 時，前兩種方法識別結果出現負數，失去物理意義，說明智能識別優化法優于這兩種方法。這一結論可從表4得到進一步證實，表4是擬合徐變函數試驗實測數據得到的識別結果。在實際數值擬合過程中，當識別參數給定初始值相同時，共軛梯度法和LM法識別結果穩定，但當給定初始值不同時，兩種方法識別結果變得不穩定。智能識別優化法識別結果基本穩定，有罰項智能識別優化穩定性更佳。

3.3.2Prony級數模型參數加載齡期譜

從3.3.1分析可知，有罰項智能識別優化方法識別結果最優。因此，基于該方法計算7個模型參數延遲強度 E₁^-1（τ₀），…，E₇^-1（τ₀） 隨加載齡期的變化規律，并使用以下3種不同模擬方法對這些規律進行模擬。具體模擬結果討論如下

1）冪函數

圖3給出了7個模型參數延遲強度隨加載齡期變化規律的冪函數擬合結果與優化識別計算結果比較。從圖中擬合效果來看，前5個模型參數幾乎完全吻合，后2個模型參數則出現些許誤差。

2）指數函數

用指數函數模擬的加載齡期對7個模型參數延期強度的影響規律與用有罰項智能識別優化法計算的結果對比結果，如圖4所示。從圖4可以看出，指數函數法的擬合效果正好與冪函數法相反，即用指數法的計算結果，前5個模型參數與優化識別計算結果之間存在些許差別，后2個模型參數則幾乎完全一致。

3）冪函數與指數函數混合函數法

為發揚冪函數法和指數函數法的優點，克服它們的缺點，前5個模型參數用冪函數、后2個模型參數用指數函數的混合擬合方法分析7個模型參數延遲強度隨加載齡期變化的規律。具體計算結果如圖5所示。從圖5比較可以看出，冪函數和指數函數混合擬合法效果最好，兩者結果吻合得非常好。

用冪函數和指數函數混合法擬合計算得到的模型參數延遲強度的表達式如下

佳擬合效果，計算得到的徐變柔量函數值與B3模型預測值具有良好一致性，相對誤差在 1% 以內。

4結論

將式（16）代入式（12），即可得到基本徐變度函數。這個基本徐變度函數很容易通過開發標準接口模塊實現與主流商業有限元軟件的無縫對接，從而提高混凝土結構徐變效應計算分析效率。

綜上所述，用指數函數法模擬徐變函數Prony級數模型中延遲強度參數效果最差，冪函數法稍好，冪函數和指數函數混合法最佳。

3.3.3 B3模型驗證

基于3種模型參數延遲強度擬合方法，圖6給出了本模型在28d加載齡期條件下14個不同持荷時間點的徐變函數計算值與B3模型的理論計算結果對比。對比分析表明：用指數函數法擬合延遲強度加載齡期譜，計算得到的徐變柔量函數值呈現較大偏差，與B3模型預測的徐變函數值相比低估 7% 左右，這可能是指數函數模擬延遲強度加載齡期譜效果差所致；用冪函數擬合延遲強度加載齡期譜，計算得到的徐變柔量函數值，仍有小幅低估，相對誤差在 3% 以內；冪函數與指數函數混合函數法表現出最

本文提出了一種用Prony級數理論模型逼近徐變函數試驗數據的有罰項混合智能優化方法（SA-GA⁺NLP ），并編制了相應的Matlab程序。通過研究分析和討論，得出以下結論。

1）通過與B3模型理論計算結果及徐變函數試驗數據對比驗證表明，相較于共軛梯度法、LM法，本文提出的含罰項混合智能優化方法具有顯著優勢：僅需簡單參數輸入即可實現計算，無需復雜非線性運算及初始值設定，有效解決了延遲強度識別結果出現負數和不穩定現象的問題，而且收斂速度快、準確、有效、可行。

2）本文提出的延遲強度加載齡期譜徐變函數Prony級數模型，基于冪-指數函數混合表達，相較于單一冪函數或指數函數的表達，對混凝土徐變函數的預測精度誤差在 1% 以內。其數學表征形式在混凝土結構徐變效應計算分析中，良好適用Abaqus，Ansys，Comsol等主流有限元分析平臺的二次開發需求。能夠有效節省計算機存儲空間和提高計算效率。

參考文獻：

[1]CHENZM，GAO F，HUJY，et al. Creep and shrinkage monitoringand modellingof CFST columnsina super high-rise under-construction building[J].Journal of Building Engineering，2023，76：107282.

[2]許開成，黃凌娟，李鵬清，等.基于縮尺模型的高鐵橋梁 預應力損失試驗研究[J].華東交通大學學報，2024，41 （2）：48-55. XUKC，HUANGLJ，LIPQ，etal.Experimental study onprestress lossofhight-speed railwaybridgesbased on scaled model[J]. Journal ofEast China JiaotongUniversity，2024，41（2）：48-55.

[3]ZHANGH，SUC，CHENXH，etal.Calculation of mass concrete temperature and creep stressunder the influence oflocal air heat transfer[J].CMES-Computer Modeling in Engineeringamp; Sciences，2024，140（3）：2977-3000.

[4]BAZANT ZP，WU S T.Dirichlet series creep function foraging concrete[J]. Journal of the EngineeringMechanics Division，1973，99（2）：367-387. tion[J].Atomic Energy Science and Technology，2024，58 （8）：1689-1696.

[14]LYUHJ，YEWL，TANYQ，etal.Inter-conversionof thegeneralizedkelvinand generalizedmaxwell modelparameters via a continuous spectrum method[J].ConstructionandBuildingMaterials，2022，351：128963.

[15]BAZANT ZP，BAWEJA S.Justification and refinements of model B3 for concrete creep and shrinkage 1.statistics and sensitivity[J].Materials and Structures，1995，28（7）： 415-430.

[16]劉成才.大跨徑鋼管砼拱橋的時間、幾何、溫度非線性空 間分析[D].重慶：重慶交通學院，2004. LIUC C.The analysis about time、geometry、temperature to long-span concrete filled steel tubular arch bridges[D]. Chongqing： Chongqing Jiaotong University.

[5]HRISTOV JY.Linear viscoelastic responses： the prony decomposition naturally leads into the Caputo-Fabrizio fractional operator[J].Frontiers in Physics，2018，6：135.

[6] LOPEZ-CAMPOS JA， SEGADE A， CASAREJOS E， et al.Behaviorcharacterization of viscoelasticmaterialsfor the finite element method calculation applying Prony series[J]. Computational and Mathematical Methods， 2019， 1（1）： e1014.

[7] XU Q W， ENGQUIST B， SOLAIMANIAN M， et al. A newnonlinear viscoelastic model and mathematical solution of solids for improving prediction accuracy[J].Scientific Reports，2020，10（1）： 2202.

[8]DIAS F，PAULLO MUNOZ LF，ROEHL D.A numerical model for basic creep of concrete with aging and damage onbeams[J].Applied Mathematical Modelling，2023， 121： 185-203.

[9]WANG Y M， SHANG L， ZHANG P P， et al. Measurement of viscoelastic properties for polymers by nanoindentation[J].Polymer Testing，2020，83：106353.

[10]黃永輝，劉愛榮，傅繼陽.基于PCG法的B3模式基本徐 變度的Dirichlet級數擬合[J].廣州大學學報（自然科學 版），2013，12（3）：30-33. HUANG YH，LIU AR，FU JY. Dirichlet progression fitting of B3 model’s compliance function for basic creep based on PCG method[J]. Journal of Guangzhou University（Natural Science Edition），2013，12（3）： 30-33.

[11] JIRASEK M， HAVLASEK P. Accurate approximations of concrete creep compliance functions based on continuous retardation spectra[J]. Computers amp; Structures，2014，135： 155-168.

[12] BARRIENTOS E， PELAYO F， NORIEGA A， et al. Optimal discrete-time Prony series fitting method for viscoelasticmaterials[J].MechanicsofTime-DependentMaterials，2019，23（2）：193-206.

[13]向華偉，榮華，范興朗，等.混凝土徐變柔度函數的高效 逼近方法[J].原子能科學技術，2024，58（8）：1689-1696. XIANG H W，RONG H，FAN XL，et al. Efficient Approximation method for concrete creep compliance func

第一作者：葛競（1999一），男，碩士研究生，研究方向為鋼-混組合結構耐久性數值分析。E-mail：gejing1204@foxmail.como

通信作者：陳夢成（1962一），男，教授，博士，博士生導師，研究方向為鋼-混組合結構設計理論及耐久性和安全性評估。E-mail：mcchen@ecjtu.edu.cn。

（責任編輯：李根）