1、ANSYS 和 MIDAS 混凝土徐变模拟比较简述：本文主要对比 ANSYS 和 MIDAS 这两种有限元软件在模拟混凝土收缩徐变上的差异，包括计算精度、计算方式、计算时间等方面。计算模型为 10m 长的 C50 方形柱顶施加1kN 的集中力，柱截面为 1m 1m。1. 混凝土徐变混凝土徐变是混凝土结构在长期荷载作用下随着时间的增长混凝土中产生的应变变化目前尚未对混凝土徐变有比较统一的说法，在此不去讨论具体有何说法，关键在于理解混凝土徐变与应力是有关系的。而通常我们计算结构时大部分是按照线性徐变处理的。2. 混凝土徐变本构关系2.1 老化理论本构关系根据迪辛格尔法可知徐变函数可定义为在 t0
2、时刻作用于混凝土的单位应力(即 t0=1)至时刻所产生的总应变。如采用徐变系数 的第一种定义，则可表示为：(,0)如采用第二种定义，则可表示为：3. ANSYS 立柱计算模型由于 ANSYS 并没有专门板块来混凝土徐变模拟，故而需要借助金属蠕变的计算机理来等效模拟混凝土徐变效应。ANSYS 提供两种方法计算徐变：显式计算和隐式计算。显式计算需要细分较多的时间步长，计算时间长；隐式计算计算精度高，计算时间短。但是在实践中也发现，涉及到单元生死情况时，隐式计算可能出现异常现象。下面将会对这两种方法进行详细的比较。3.1 ANSYS 显式计算显式计算对时间步长是有要求的，尤其是在徐变系数曲线变化剧烈
3、的时间段需要细分子步以减小误差和帮助收敛。因而，时间步长的划分方式、时间点的数目对计算结果都会有较大的影响。（1） 等间距时间步长和对数时间步长假设混凝土的龄期是 7 天，徐变变化速率为 0.005，考虑收缩徐变 10 年（3650 天） ，若 3650 天时刻的徐变系数为 1，那么按照等间距时间步长划分，则时间步长间距， （3650-7）/500=7.286。按照对数时间步长划分，若采用 30 个数据点，具体数据如下所示。表 1 对数时间步长数据表编号 时间 编号 时间 编号 时间 编号 时间 编号 时间1 7 11 24.46454 21 85.50194 31 298.8236 41 1
4、044.3692 7.933082 12 27.7256 22 96.89913 32 338.656 42 1183.583 8.990541 13 31.42135 23 109.8155 33 383.798 43 1341.3494 10.18896 14 35.60973 24 124.4537 34 434.9573 44 1520.1475 11.54712 15 40.35642 25 141.043 35 492.936 45 1722.7796 13.08632 16 45.73582 26 159.8437 36 558.6431 46 1952.4217 14.8306
5、9 17 51.83229 27 181.1504 37 633.1087 47 2212.6738 16.80758 18 58.7414 28 205.2973 38 717.5005 48 2507.6179 19.04799 19 66.57147 29 232.6629 39 813.1415 49 2841.87610 21.58704 20 75.44528 30 263.6763 40 921.5311 50 3650（2） 徐变变化速率徐变变化速率影响着徐变变化曲线的陡缓，将会对不同徐变变化速率值进行比较：0.001、0.002、0.003、0.004、0.005、0.006
6、、0.007、0.008、0.009、0.010。（3） 计算结果对比两种时间步长划分方式和不同徐变变化速率柱顶徐变 10 年位移计算结果如下表所示：表 2: 柱顶徐变 10 年计算结果对比（单位:m）等间距时间步长 对数时间步长徐变变化速率 理论值计算值 相对误差 计算值 相对误差0.001 5.721e-7 5.71e-7 0.19% 5.66e-7 1.07%0.002 5.795e-7 5.78e-7 0.26% 5.74e-7 0.95%0.003 5.797e-7 5.77e-7 0.47% 5.75e-7 0.81%0.004 5.797e-7 5.77e-7 0.47% 5.7
7、5e-7 0.81%0.005 5.797e-7 5.76e-7 0.64% 5.74e-7 0.98%0.006 5.797e-7 5.75e-7 0.81% 5.74e-7 0.98%0.007 5.797e-7 5.74e-7 0.98% 5.74e-7 0.98%0.008 5.797e-7 5.74e-7 0.98% 5.74e-7 0.98%0.009 5.797e-7 5.73e-7 1.16% 5.74e-7 0.98%0.010 5.797e-7 5.72e-7 1.33% 5.74e-7 0.98%从上表可以看出，不同的时间步划分方式对结果产生较大的差异，等距时间步长随着徐
8、变变化速率的增大精度不断降低，对数时间步长则随着徐变变化速率的增大精度先上升后趋于平稳。且对数时间步长的计算时间要短，精度也能满足工程要求，且在较大徐变变化速率区间，采用对数时间步长更容易获得较好结果。故而，建议采用对数时间步长进行显式计算。3.2 ANSYS 隐式计算隐式计算也需要区分等间距时间步长和对数时间步长两种划分方式。不过这里不对隐式计算进行详细的探讨，光是一个简单的立柱还不能够说清楚 ANSYS 模拟徐变问题。下面将会讨论考虑施工过程的两端固结梁的徐变问题。4. ANSYS 两端固结梁计算模型只考虑显式计算。仍然考虑等间距划分和对数划分两种方式，比较这两种划分方式的计算结果、收敛情
9、况、计算耗时等。此处附上显式计算命令流，不作详细说明，因为显式计算碰到收敛问题无法解决！! 计算悬臂梁转化为固定梁的徐变效应finish/clear/prep7Ec=3.45e10! time array*dim,tt,array,50,1*vread,tt(1),tt1,txt(1F8.3)! creep coefficient=0.005nn=50*dim,fi,array,nn,1*do,i,1,nnfi(i)=1-2.7182818*(-0.005*(tt(i)-7)*enddo*dim,C1,array,nn,1f1=0k1=3dt=0*do,i,1,nnC1(i)=(fi(i)-f
10、1)/(1+fi(i)/(tt(i)-dt)f1=fi(i)dt=tt(i)mp,ex,k1,Ecmp,prxy,k1,0.2tb,creep,k1tbdata,1,C1(i),0,1,0k1=k1+2*enddo*dim,C2,array,nn,1f1=0k2=4dt=0*do,i,1,nnC2(i)=2*(fi(i)-f1)/(1+2*fi(i)/(tt(i)-dt)f1=fi(i)dt=tt(i)mp,ex,k2,Ecmp,prxy,k2,0.2tb,creep,k2tbdata,1,C2(i),0,1,0k2=k2+2*enddoet,1,23keyopt,1,6,0r,1,1,1/1
11、2,1! 左边材料属性mp,ex,1,Ecmp,prxy,1,0.2tb,creep,1tbdata,1,0,0,0,0! 右边材料属性mp,ex,2,Ecmp,prxy,2,0.2tb,creep,2tbdata,1,0,0,0,0! 建立有限元模型*do,i,1,19n,i,(i-1)/2*enddo*do,i,1,19n,i+19,(i-1)/2+9*enddotype,1mat,1*do,i,1,18e,i,i+1*enddotype,1mat,2*do,i,18,35e,i+2,i+2+1*enddod,1,alld,38,allesel,s,1,18sfbeam,all,pres,
12、1000esel,s,19,36sfbeam,all,pres,1000/soluoutres,all,allallsel,allcrplim,0.25bfunif,temp,100time,1e-6solvek1=3k2=4*do,i,1,nnesel,s,1,18mpchg,k1,allk1=k1+2allsel,allesel,s,19,36mpchg,k2,allk2=k2+2allsel,alltime,tt(i)nsubst,10solve*enddo隐式计算命令流如下：! 计算悬臂梁转化为固定梁的徐变效应finish/clear/prep7Ec=3.45e10! time arr
13、ay*dim,tt,array,50,1*vread,tt(1),tt1,txt(1F8.3)! creep coefficient=0.005nn=50*dim,fi,array,nn,1*do,i,1,nnfi(i)=1-2.7182818*(-0.005*(tt(i)-7)*enddo*dim,C1,array,nn,1f1=0k1=3dt=0*do,i,1,nnC1(i)=(fi(i)-f1)/Ec/(tt(i)-dt)f1=fi(i)dt=tt(i)mp,ex,k1,Ecmp,prxy,k1,0.2tb,creep,k1,11tbdata,1,C1(i),1,0,0,0k1=k1+2
14、*enddo*dim,C2,array,nn,1f1=0k2=4dt=0*do,i,1,nnC2(i)=2*(fi(i)-f1)/Ec/(tt(i)-dt)f1=fi(i)dt=tt(i)mp,ex,k2,Ecmp,prxy,k2,0.2tb,creep,k2,11tbdata,1,C1(i),1,0,0,0k2=k2+2*enddoet,1,188sectype,1,beam,rect,a1secdata,1,1keyopt,1,1,0keyopt,1,3,3! 左边材料属性mp,ex,1,Ecmp,prxy,1,0.2tb,creep,1,11tbdata,1,0,0,0,0,0,0! 右
15、边材料属性mp,ex,2,Ecmp,prxy,2,0.2tb,creep,2,11tbdata,1,0,0,0,0,0,0! 建立有限元模型*do,i,1,19n,i,(i-1)/2*enddo*do,i,1,19n,i+19,(i-1)/2+9*enddotype,1mat,1*do,i,1,18e,i,i+1*enddotype,1mat,2*do,i,18,35e,i+2,i+2+1*enddod,1,alld,38,allesel,s,1,18sfbeam,all,pres,1000esel,s,19,36sfbeam,all,pres,1000/soluoutres,all,allrate,offallsel,allcrplim,0.25bfunif,temp,100time,1e-6solvecpintf,allrate,onk1=3k2=4*do,i,1,nnesel,s,1,18mpchg,k1,allk1=k1+2allsel,allesel,s,19,36mpchg,k2,allk2=k2+2allsel,alltime,tt(i)nsubst,40solve*enddo5. MIDAS 计算模型
