摘要：通过对混凝土粘弹性相的任意徐变度函数C(t-t′)进行Laplace变换，导出了混凝土材料的连续阻尼谱函数，可以使Bazant，Z.P.固化徐变模型适合于具有任意徐变规律的混凝土材料，避免了Bazant算法中根据经验选择阻尼时间的做法。
关键词：固化徐变理论 离散阻尼时间 连续阻尼谱 指数算法
1996年以前，工程上对温度徐变应力的计算一般采用“松弛系数法”[1]，这对于均质或满足比例变形条件的非均质是合适的。1996年以后，新颁《混凝土设计规范》[2]建议考虑大坝的非均质性以及材料参数的时间依赖性，按混凝土徐变度方程计算温度应力，即根据不同的加载过程，以适当的时间步长，利用线性徐变理论的叠加原理，逐步计算坝体温度应力。这类方法以文献[3]提出 的“混凝土徐变应力分析的隐式解法”最为成熟与实用，它既能有效地节省存储空间，也能考虑时间步长的变化，为大体积混凝土计算、的仿真计算奠定了基础。
混凝土非线性徐变理论的研究在我国六十年代初就开始了[4～6]，但还没有形成统一的、有影响力的理论。所以本文的研究放在已在国际上产生影响的Bazant Z. P.非线性徐变理论-混凝土固化徐变理论[7]的基础上。
1 Bazant固化徐变理论
1.1 混凝土粘弹性相徐变及其求解 Bazant固化徐变理论是将弹性理论、粘弹性理论与流变理论结合起来，模拟由于水泥不断水化、固相物不断增多、混凝土宏观物理力学性质随时间不断变化的新理论。这一理论最大的特点是将混凝土宏观材料参数对时间的依赖性，归结为混凝土材料的粘性相与粘弹性相体积不断增多(粘性相与粘弹性相的物理性质不变)、非承力相体积(如孔隙、胶体、水等)不断固化的结果(弹性相体积不变)，因此也称为混凝土固化徐变理论。该理论与用某一类函数模拟宏观上混凝土徐变度的做法不同，是从微观物理概念出发，直接推导出宏观上混凝土徐变度的表达式，导出了徐变应力控制方程。
—在Bazant固化徐变应力控制方程中，在任意时刻，混凝土的总应变ε应满足：
ε=σ/E0+εc＋ε0，εc＝εv＋εf (1)
公式中，εc为混凝土的徐变应变，εv为混凝土粘弹性相徐变，εf为混凝土粘性相流动徐变。ε0为各种附加应变，包括混凝土的自生体积变形、混凝土的温度变化、混凝土微裂缝的扩展等引起的应变。σ/E0为混凝土弹性相应变。式(1)中，除了εv比较复杂外，其它应变都比较简单，不是本文研究的对象。
混凝土粘弹性相徐变εv没有龄期效应，只与持荷时间有关，可以用一系列串联的Kelvin固体单元来模拟[8-9]。根据Kelvin固体的串联模型，第μ个Kelvin单元的平衡条件为：
(2)
式中：Eμ、ημ分别为第μ个Kelvin单元的弹性模量和粘滞系数，rμ为第μ个Kelvin单元的应变，r为粘弹性相的总应变，σ为混凝土宏观应力。将式(2)分别求解，然后再求和，得到在不变应力作用下，混凝土粘弹性相任意时刻的应变为：
(3)
这时，如果混凝土粘弹性相徐变度函数服从对数幂函数分布［10］，就可以用快速收敛Dirichlet级数来逼近它，即令：
(4)
式(3)与式(4)表达的物理意义相同，形式相当，两者只在常数项有区别，其转化关系为Eμ=1/q2Aμ。常数Aμ需要根据试验资料按最小二乘法确定；阻尼时间常量τμ如果也由试验资料确定时，将导致一个病态方程组的求解[11]，最好根据计算经验取值。根据Bazant的计算经验，第1个Kelvin单元的阻尼时间τ1及Kelvin单元的个数N要根据我们感兴趣时间范围来选择，尤其是τ1的选择，要经过试算，第μ个Kelvin单元的阻尼时间τμ则可取为对数时间坐标，即τμ=τ110μ-1(μ=1，2，…，N)。当τμ、Eμ一定，第μ个Kelvin单元的粘滞系数ημ也就完全确定了，即
ημ＝Eμτμ (5)
1.2 混凝土粘弹性相徐变度函数服从对数幂函数分布时的有关系数 现在要针对具体材料徐变度函数分布，确定算法中的有关系数。首先给出对数函数logξ及指数函数ξn的Dirichlet级数展开式。
(6)
(7)
式中：bμ(n)为查表算得的常数[8]。
对于对数幂函数ln(1+ξn)，当ξ《1时，ln(1+ξn)≈nlnξ；当ξ1时，ln(1+ξn)≈ξn。为了得到ln(1+ξn)的Dirichlet级数展开式且符合Kelvin固体的一般规律，Bazant教授将式(4)改写为：
式(8)、(9)、(10)中的有关系数，如c、z、bμ等均按试验参数，利用Levenberg Marquardt算法[12]优化而得。其中：z、b1与n的关系见表1.当0.05≤n≤0.25时，在我们感兴趣的时间范围内，如0.25τ2≤ξ≤0.5τN，Dirichlet级数逼近原函数的误差在1%以内[7]。
表1 函数ln(1+ξn)的Dirichlet级数展开式中的两个系数
2 混凝土徐变度的连续阻尼谱函数
如前所述，在不变的单位应力作用下，混凝土的柔度函数为：
J(t，t′)=q1+C(t，t′) (11)
式中：q1为瞬时弹性应变，C(t，t′)为混凝土的徐变度。对于没有老化性质的Kelvin固体，徐变度仅为持荷时间的函数，即
（12）
q2＝1，为了避免阻尼时间τμ选择上的任意性，令[13]：
（13）
其中，L*(τ)=L(τ)/τ，L(τ)为徐变度函数的连续阻尼谱函数。将它代入上式，有：
（14）
再令τ=1/ζ，则d(lnτ)=-d(lnζ)，式(13)即转换为：
如果记
（16）
则式(15)就变成：
C(ξ)＝-f(ξ)+f(0) (17)
显然，f(ξ)为核函数ζ-1L(ζ-1)的Laplace变换式，ξ为转换变量。对式(17)进行Widder[14]变换，即得:
（18）
且
（19）
利用式(17)，f(0)为常数，便有
（20）
L(τ)即为待求的混凝土徐变度的连续阻尼谱函数。C(k)(kτ)为徐变度函数的k阶导数。将式(20)代入式(14)就可以用连续阻尼谱表示混凝土的徐变度了。这与式(4)表示徐变度的离散方法是完全不同的，它适合各种徐变度函数，但其基本前提是混凝土徐变度拟合函数的k阶导数存在。
3 对数幂函数表达的徐变度函数的连续阻尼谱的离散方法
式(20)对非老化材料是普遍适用的。由于在求解大型的徐变应力问题时，常采用有限差分方法，所以，在时间域上还需对式(14)进行离散，方可将式(14)应用于数值分析。下面以对数幂函数表达的徐变度函数为例，说明离散的方法。
如果公式C(ξ)=q2ln[1+(ξ/λ0)n］中的λ0＝1d[7]，那么，徐变度函数变为：
C(ξ)=q2ln(1+ξn) (21)
对于k=3，按式(20)求得的对数幂函数表达的徐变度函数的连续阻尼谱为：
（22）
由于n为很小的正常数，式(22)可以简化为：
（23）
将式(14)的时间lnτ离散，Δ(lnτμ)=ln10Δ(logτμ)，并将积分以求和近似代替，则
（24）
Aμ=L(τμ)ln10Δ(logτμ) (25)
在式(24)、(25)中，ξ为混凝土的持荷时间，是由数值计算的时间步长确定的。τμ是混凝土的阻尼时间，是反映材料徐变特性的一种参数。数值计算中，可根据公式(24)拟合连续函数的光滑程度取值。经验表明，当Δ(logτμ)=1时，正好取对数时间坐标，曲线也足够光滑，见图1。
图1 两种阻尼谱公式对徐变度函数的拟合效果
4 连续阻尼谱函数的合理性检验
为了比较Bazant离散阻尼谱函数和本文提出的连续阻尼谱函数对混凝土粘弹性相徐变度的拟合效果，就文[15]求出的沙牌碾压混凝土的拟合系数q2，在0.1～10000d时间范围内，计算了公式(8)、(21)、(24)的具体值，见表2.其图形表示见图1.总的感觉是由Bazant公式(8)逼近公式(21)(24)拟合公式(21)有稍高的精度。
在利用公式(8)时，取τ2=1，N=7，τ1=0.1τ2；b1、z按n=0.1在表1中线性插值，在利用公式(24)时，只涉及到n和q2，与τμ的取值无关。但为了和式(21)、式(8)比较，在图1中，取对数阻尼时间步长Δ(logτμ)=log10=1。
从表2和图1均可看出，在1≤ξ≤10000(d)范围内，将连续函数展开成Dirichlet级数的误差在4%以内；在1≤ξ≤1000(d)范围内，误差小于1%，这与Bazant的研究结果吻合；当ζ&<1d或ξ&>10000d时，连续阻尼谱函数公式的误差较大，这一问题有待进一步研究。
同时，为了说明离散阻尼谱函数的非唯一性，我们就τ2＝1和τ2＝0.1这两种情况，在τ1≈(10-5～10-1)τ2范围内讨论了τμ的取值对Dirichlet级数精度的影响，见图2和图3。
观察图2和图3就可发现，如果要用离散的阻尼谱函数拟合混凝土粘弹性相的徐变规律，需要采用“试算法”，当τ1取值合适，效果将很好，否则达不到1%的精度。从这一意义上讲，连续阻尼谱合适的优越性是非常明显的。
说明1：式(8)，τ1=0.1τ2；2:式(8)，τ1=0.01τ2；3:式(21)；4:式(8)τ1=0.001τ2；5:式(8):τ1＝0.0001τ2
图2 τ2=0.1时公式(8)对公式(21)的拟合效果比较
说明:1：式(18)，τ1=0.1τ2；2：式(8)τ1＝0.01τ2；3式(21)；4:式(8)，τ1=0.001τ2；5：式(8)，τ1=0.0001τ2
图3 τ2=1时公式(8)对公式(21)的拟合效果比较
5 小结
为了避免Bazant算法中根据经验选择Kalvin固体的阻尼时间的做法，通过对粘弹性相的任意徐变度函数C(t-t′)进行Laplace变换，导出了混凝土材料的连续阻尼谱函数，并根据沙牌碾压混凝土徐变度函数的拟合资料，检验了这种方法的合理性。这样，就可以使Bazant，Z.P.固化徐变模型适合于具有任意徐变规律的混凝土材料。