前言：中文期刊网精心挑选了数值模拟范文供你参考和学习，希望我们的参考范文能激发你的文章创作灵感，欢迎阅读。
数值模拟范文1
2数值模型
2•1基本方程本文的模拟对象是铸轧过程中铝石墨半固态浆料的不可压缩稳态流动和传热的耦合问题.对于热流耦合这样的连续性问题,数学模型由连续性方程、动量守恒方程、能量守恒方程来描述.连续性方程
2•2凝固潜热的处理半固态金属铸轧成形过程中伴随凝固相变,因而不但要考虑传热,还要考虑凝固潜热的释放.本文采用等效比热容法,处理凝固潜热的释放问题.等效比热容法是假设凝固潜热在两相区之间释放平均,将凝固潜热处理为比热容的一部分.定义如下式中:L为凝固潜热,J/kg;ts,t1分别为凝固结束和开始温度,℃;C1,C2分别为固态和液态的真实比热容,J/(kg•℃).
2•3流变模型半固态金属的流变行为主要受剪切速率、固相体积分数和固相微粒形态的影响,而固相体积分数、固相微粒的形态又受到温度的影响,各种影响关系错综复杂.研究表明,半固态铝合金具有显著的剪切变稀特性,是一种非牛顿流体,其表观黏度与剪切速率的关系服从Power-law定律[5].采用了Carreau模型描述半固态金属的流变行为,其表达式如下
3几何模型和边界条件基于前述数学模型,结合实验过程中各种实际工艺参数及实验设备自身性能,本文采用有限元软件ANSYS对半固态铸轧复合过程进行了数值模拟.模拟对象的材料分别为08Al钢板和固相率为27%~46%的铝石墨半固态浆料.材料物性参数根据实际情况取值,其中铝石墨的凝固点为580℃.考虑到轧辊直径远大于复合板厚度,宽展忽略不计,仅以铸轧区域沿轧制方向纵切面建立二维模型作为研究对象.
3•1几何模型和模拟条件根据上述数值模型,以浇嘴内的扩散区和进入辊间后的铸轧区为计算域,对半固态铸轧复合工艺进行合理简化后建立的几何模型及划分的网格如图2所示.研究所用的模拟条件如下:轧辊直径320mm,辊缝宽度2•5mm,钢板厚度1•2mm,钢板预热温度510℃,浇注温度分别为610℃、620℃和630℃,铸轧速度0•4~1•0m/min.
3•2基本假设和边界条件铸轧复合过程中,铝石墨半固态浆料分别与轧辊和钢板发生换热并在轧辊出口处附近凝固,在铸轧区内,金属熔体同时存在固态与半固态两种形态,且涉及相变过程,传热流动现象极为复杂.为了建立描述铸轧过程中热流耦合问题的计算模型,适当简化了计算条件,采用了如下的基本假设:①铸轧过程是稳态进行的,经过初始过渡期后,工艺参数不随时间变化;②铝合金半固态浆料是不可压缩的,视为剪切变稀的非牛顿流体;③轧辊和轧件之间无相对滑动.边界条件的处理如下:1)入口边界.入口处半固态浆料流速为注流速度,温度为浇注温度,钢板温度为预热温度,浇口侧面与空气之间绝热,属于第一类边界条件:u,vx,y;uin口注流速度;Tin为半固态浆料浇注温度.2)铝石墨与浇嘴接触边界.因为浇嘴外包覆了绝热层,所以认为浇嘴对铝熔体有理想保温作用,此处设置为绝热边界浇嘴材料选用石墨,石墨浇嘴与铝熔体几乎不润湿,可有效减少流动阻力,因此在半固态浆料与浇嘴的接触面设置了壁面接触角为100°.3)铝石墨和轧辊的接触边界.根据基本假设,熔体与轧辊接触处的节点速度与轧辊表面线速度相同,方向为接触边界切线方向.此处为热流连续,但温度不连续的接触边界,属于第三类边界条件式中:x0,y0为辊心在x,y方向的坐标;κ为铝石墨导热系数;α为接触面对流换热系数;ω为轧辊角速度;Troll为辊套外表面温度.4)铝熔体与钢板接触面.此处属于第三类边界条件-κTn=α(T-Tsteel)(10)式中Tsteel为钢板上表面温度.5)钢板下表面冷却区.钢板下表面受到水雾冷却,此处设置为恒定的负热流密度,属于第二类边界条件-κTn=qw(11)式中qw为钢板与水雾接触面热流密度.6)铝熔体出口边界.此处为绝热边界式中ue为复合板出口速度.
4结果分析
4•1浇注温度对温度场的影响铸轧速度为0•7m/min,铝石墨半固态浆料浇注温度分别为610℃、620℃、630℃的条件下,计算得到的复合板铸轧区温度分布如图3所示.铝石墨半固态浆料进入铸轧区后,由于和水冷轧辊紧密接触,温度快速下降,靠近轧辊处首先形成凝固壳.从图3中可以看出,在凝固点以下,铝石墨温度梯度增大,温度的下降较未凝固前更为迅速.分析认为,这是由于熔体在到达凝固点之前需要释放大量凝固潜热,因而在相同的冷却条件下,凝固后的熔体相对于半固态熔体得以更快地释放自身热量.这一现象符合事实情况,也验证了等效比热法的准确性.在轧制出口处,钢板和铝石墨的温度基本一致.随浇注温度的升高,熔池内凝固界面的位置在逐渐向出口方向移动.图3(c)是浇注温度630℃时的温度分布云图,由于浇注温度过高,铝石墨在出轧辊时仍未凝固.如在此浇注温度下开展铸轧试验,将
4•2铸轧速度对温度场的影响保持半固态浆料浇注温度620℃,分别以铸轧速度0•4m/min、0•6m/min、0•8m/min、1•0m/min进行模拟,获得铸轧区温度分布云图如图4示.铸轧区总长度25mm,在铸轧速度为0•4m/min时,凝固前沿在铸轧区入口7mm处,此时凝固点过高,凝固壳较厚,容易导致铸轧力过大,可能出现轧卡现象,使得铸轧过程不能继续进行.在铸轧速度为1•0m/min时,由于铸轧速度过高,出辊时复合界面附近的铝石墨尚未凝固.根据模拟计算结果,铸轧度为0•6~0•8m/min时,凝固前沿在距离铸轧区入口13~20mm的范围内,既避免了凝固不足,也不会因凝固壳太厚导致轧卡,可以保证铸轧过程定进行,此计算结果与实验数据相吻合[6].图5对应于不同铸轧速度的复合界面上节点的温度分曲线图.铸轧速度越小,温度曲线斜率越大,即温下降速度越快.说明在铸轧辊冷却强度和铸轧区度不变的情况下,同样的浇注温度,铸轧速度越低,浆料与轧辊接触的时间越久,轧辊里的冷却水会走更多热量.同时,当复合界面温度下降到凝固点下,温度下降速度加快,单个曲线斜率增大,在铸速度为0•6m/min时较为明显.
数值模拟范文2
关键词：旧桥加宽；数值模拟；承载力验算
方案介绍
某项目沿线有一座右偏角为120度的一孔13m先张简支空心板，原路线全宽16.5m，需拓宽至25m，沿线桥梁均需加宽8.5m，13m先张简支空心板原规范标准板梁高度为55cm，新建部分按新规范标准板梁高为70cm，本文针对加宽后桥梁全桥结构采用梁格法建模，对全桥的承载力进行验算。推荐方案：原桥上、下部结构全部利用，只新建左幅8.5m宽上部及下部结构，新建下部结构与原桥下部结构齐平，原桥部分桥面通过调整调平层达到设计高程，新旧板之间通过植筋可靠连接；比较方案：原桥下部结构全部利用，上部结构右半幅利用，新建左半幅上部及下部结构，新建下部结构与原桥下部结构齐平，原桥部分桥面通过调整调平层达到设计高程，由于新、旧板铰缝处为中央分隔带，其间可留1cm的缝隙。两方案综合比较：经济合理的方案为将原桥完全利用，只在左侧新建需加宽部分桥梁，方案图如下：
主要材料指标
模型简介
该模型梁单元数量总计1501个，节点数量：872个。模型横断面图如下图：
模型横断面图
模型基本理论介绍
该桥右偏角为120度，与正桥相比，斜桥还具有如下特性：（1）沿宽度最大弯矩方向的变化，在边缘处与斜跨方向平行，在板的中央则接近垂直于桥台；（2）靠近钝角处出现上拱弯矩；（3）上部结构承受很大的扭转；（4）钝角角隅处出现较大的反力和剪力；（5）锐角角隅处出现较小的反力，还可能出现翘起。
分析结果
1.反力：在荷载作用下，钝角处反力较大，对于这种宽度大于上部结构跨度的结构，可以降低支座的柔性来调节各个支座反力的大小。
自重作用下各支座反力自重+二期恒载+公路1级车道荷载支座反力
2.挠度：以下两图挠度单位为mm，从下两图看出，旧板挠度与新板相比较大，在新旧桥相接处，板的挠度与其两侧板挠度相比略大，在旧桥加宽后，此桥按公路-Ⅰ级验算得新旧桥相接处跨中最大挠度为+0.9mm，新规范要求最大挠度不超过跨径的1/600，即21.7mm，可知旧桥位移满足现行规范的要求。
自重作用下挠度图 长期组合作用下挠度图
3.应力：旧桥按新规范荷载加载，按弹性阶段验算，梁单元跨中下缘仍为压应力，有一定的安全储备，应力验算也可满足现行规范的要求。
弹性阶段组合梁单元上缘应力验算弹性阶段组合梁单元下缘应力验算
结语
通过对同跨径不同梁高旧桥加宽后的桥进行验算后可知：按旧规范设计的13m跨径的空心板满足现行规范移动荷载为公路-Ⅰ级的要求，可完全利用，本结果对于其余跨径的简支空心板也有一定的参考价值。
参考文献：
[1]项海帆.高等桥梁结构理论[M]，北京：人民交通出版社，2001.
[2]汉勃利.桥梁上部结构性能
[3]李满红，肖清亮;浅谈改建公路中加宽桥梁的合理设计[J];东北公路;1994年04期
[4]叶生;旧桥整体加宽中若干问题的研究[D];合肥工业大学;2006年
作者：白延芳
学历：硕士研究生
数值模拟范文3
关键词：叶轮；浇注系统；数值模拟；有效应力
中图分类号：TG249.5 文献标识码：A 文章编号：1000-8365（2016）04-0805-04
熔模铸造工艺，被称为失蜡铸造或者精密铸造，是一个生产接近最终成型部件的完整体系。这种制造技术已经广泛地应用于铸造各种壳体、飞机配件等承受载荷的零件。对零件缩减重量和提高效率的研究使部件形状更加复杂，也更薄，也让熔模铸造变得更加复杂。复杂形状的零件熔模铸造的发展是通过尝试错误法或者是通过反复铸造试验来实现的，铸造参数在优良品质的铸件生产出来之前一直是变化的，但是这种方法非常昂贵也很费时[1-3]。数值模拟是另外一种确定铸造参数的方法，如果加工材料的准确性能数据已知，并且精确定义边界条件，铸造模拟将很可靠[4]。本文以某型铝合金叶轮为例，对铸件设计了侧注式和顶注式浇注系统，利用ProCAST有限元铸造模拟软件对其铸造过程进行数值模拟，有效地预测了铸造缺陷，从而为确定最佳工艺方案、优化铸造工艺参数，确保铸件质量，缩短产品试制周期，降低生产成本，提高生产效率提供了科学的依据。
1铝合金叶轮的结构特点
所研究的铝合金叶轮为薄壁零件，采用A356铝合金熔模精铸而成，其结构如图1。叶轮壳体内径108mm，外径112mm，厚度仅为2mm，8个叶片呈扭转形状均匀分布在壳体四周，叶片厚度不均匀，最薄处仅为1.2mm。因为叶片结构容易引起充型过程不平稳，产生紊流，出现缩松缩孔影响力学性能，应力过大导致裂纹缺陷，合格率较低，所以有必要对叶轮进行数值模拟研究，设计合理的浇注系统，为实际生产提供理论依据。
2原始设计方案的模拟与分析
2.1边界条件和热物性参数设计
启动软件VisualEnvironment8.0，在Vi-sualMesh模块，导入叶轮的igs格式文件，设置铸件和浇注系统网格长度均为2mm，划分好网格的浇注系统如图2，浇注合金为A356，耐火型壳材料采用FusedSilica，接触面类型为COINC，界面换热系数为1000W/m2•K，冷却类型Aircooling，浇注类型为GravityFilling。其热物性参数见表1。设置浇注温度为750℃，浇注速度为0.3m/s，型壳预热温度为400℃。
2.2数值模拟结果及分析
顶注式填充过程如图3。合金液首先到达型腔底部，逐渐向上充型的过程中同一水平位置的充型状态并不同，合金液面高低不平，这说明铝合金液面的平稳性较差。尤其在叶轮中间部位，紊流现象明显，容易造成合金液的氧化以及氧化物卷入合金液内部，影响铸件的质量。通过温度场缩松缩孔模拟结果可以看出，叶片心部出现了体积很大的缩孔缩松缺陷(如图4)，叶片根部与壳体连接部位1、2、3处也出现了221MPa的应力(如图5)，如果超过材料的屈服强度，有可能会出现微裂纹。如果叶片心部组织中存在较大缩松缩孔和微裂纹，会导致局部的应力集中，增大截面上所承受应力，因而在使用过程中，诱发早起的开裂，形成早期的裂纹源[9]。实际生产的带有缺陷的叶轮如图6。
3改进方案设计及模拟结果
从原始方案的设计和模拟结果可以看出，应当从三个方面进行优化：①金属液直接冲击叶片边缘，易引起金属液飞溅，产生紊流；②叶片中间部分有较大的缩孔缩松；③叶片根部和主壳体之间存在较大的应力，容易引起热裂缺陷。基于以上考虑，主要通过改变浇注系统来获得较好的充型。图7采用了侧注式浇注系统，在壳体上中下部位各设置一个浇口，模拟中的试验参数与原方案相同。改进后的浇注系统如图8所示，侧注式方案设计使得合金液在注入型腔时得到了一定程度的缓冲，对腔壁的冲击较小，无明显局部紊流现象发生。合金液从底部横浇道流入型腔后，开始依次向上充填进叶轮底部，与中部横浇道中合金液体汇合后再向上充型，最终充满整个型腔，实现了顺序充型。上中下横浇道同时注入较单一注入效率高，并且有利于气体的排出。在叶片边缘以未观察到明显的紊流，合金液面上升平稳，因此认为该方案的充型过程良好。横浇道和内浇口的设计有助于顺序凝固，让铸件得到有效的补缩，减小了叶片产生的缩孔缩松。由图9中可以看出，壳体的缩松缩孔缺陷较顶注式方案时减小，在叶片边缘与主壳体的结合部位有少量缩松缩孔出现，最大缩孔率约为0.69，铸件本体基本没有宏观缩孔缩松缺陷出现。从图10可以看出，叶片根部与壳体之间在凝固末期的有效应力降低为158MPa，有效减少了裂纹的产生。
4实验验证
基于仿真分析结果，采用侧注式的型壳和成品如图11和图12，铸件充型完整、轮廓清晰、表而质量好，无缩孔、缩松等明显的铸造缺陷。
5结论
（1）使用同样浇注参数的侧注式浇注生产的铝合金叶轮相比于顶注式浇注充型过程平稳，缩松缩孔体积较小，减小了叶片和壳体之间的有效应力。（2）实践表明，应用ProCast软件模拟铸件凝固过程，可以准确地预测缺陷类型、大小及位置、分析有效应力，为选择高效工艺方案提供了可靠的依据。
参考文献：
[1]PattnaikS,KarunakarDB,JhaPK.Developmentsininvestmentcastingprocess-Areview[J].JournalofMaterialsProcessingTech-nology,2012,212(11):2332-2348.
[2]ReedRC.Thesuperalloys:fundamentalsandapplications[M].Cambridgeuniversitypress,2006.
[3]MandziejST.Low-cyclethermal-mechanicalfatigueasanacceler-atedcreeptest[J].ProcediaEngineering,2010,2(1):349-358.
[4]HamiltonRW,SeeD,ButlerS,etal.Multiscalemodelingforthepredictionofcastingdefectsininvestmentcastaluminumalloys[J].MaterialsScience&EngineeringA,2003,343(s1-2):290-300.
[5]衣春雷.熔模铸造薄壁铸件工艺优化的实验及数值模拟研究[D].青岛:青岛理工大学,2012.
[6]张洁,王刚,郑福生,等.基于ProCAST的高压水泵壳熔模铸造工艺研究[J].铸造,2012，61(11):1324-1326.
[7]张培柳,曹韩学,郑会芬.基于数值模拟的箱盖件浇注系统的设计与优化[J].热加工工艺,2013，42(01):35-37.
[8]杨亮,张来启,林均品.高Nb-TiAl增压涡轮熔模铸造过程数值模拟[J].热加工工艺,2014,43(19):72-74.
[9]王芳,李宝宽,沈峰满.离心式中间包流场的大涡模拟研究[J].东北大学学报(自然科学版),2007(08):1143-1146.
数值模拟范文4
（青岛大学，山东 青岛 266071）
【摘 要】冷藏箱箱体内温度场与流场的分布均匀性与箱体的结构有着密切的关系，数值计算方法能够为冷藏箱结构优化提供一种省时省力的重要工具。近年来，随着计算机技术的快速发展，计算流体软件CFD越来越多的应用于冷藏箱的结构优化中。本文介绍了指出了冷藏箱研发以及使用当中遇到的一些问题，并介绍了国内外研究人员使用CFD优化箱体结构的经验。
关键词 冷藏箱；结构优化；数值计算方法；温度场；CFD
0 引言
随着社会的发展，人们生活水平的提高，对冷藏箱的需求也越来越大，同时对冷藏箱的要求也越来越高，不仅要求能够更好的保存物品的质量，而求要求能够节能和环保。如何设计能够更好地满足用户要求的冷藏箱是各个生产厂家面临的主要问题。目前冷藏箱仍然面临着一些问题，比如箱体内部温度场不均匀，某些部位温度超标，影响了物品的保存；压缩机的频繁开停机以及外界环境通过壁面和门封向箱体内部的漏热导致冷藏箱耗能的增加等问题，严重影响着冷藏箱的发展。因此如何更加有效快速的对冷藏箱箱体的结构进行优化，来改善冷藏箱的性能，是国内外学者以及公司设计人员的研究重点[1-2]。
许多企业在设计冷藏箱时还是采用传统的实验方法，不仅耗时耗力，而求很难达到预期的效果。随着计算机技术与数值传热学的结合，计算流体力学软件CFD得到快速的发展，也已运用到了冷藏箱的设计以及结构优化中去，不仅能够缩短设计周期，而且节省了大量的财力。目前，国内外的研究人员在设计冷藏箱时应用CFD有了很多的研究成果。
1 冷藏箱箱体内温度场与流场分析
冷藏箱箱体内温度场与流场的分布是否均匀，不仅影响着保存物品的质量，而且影响着冷藏箱的功耗。温度场与流场的优化如果仅仅靠做实验的方法来进行，不仅周期较长，而且需要大量的人力物力，不符合企业的利益；因为与温度场和流场的分布有关参数有多个，工作量太大，而且如果要测量流场需要布置大量的传感器，这样做又破坏了箱体内温度场与流场的分布[3]。数值模拟仿真能够较真实的模拟箱体内的温度场与流场的实际分布情况，已成为研究人员进行流场计算与优化的重要工具，为冷藏箱箱体的结构优化提供了依据。
为了能够使模型简单，考虑到箱体结构的对称性，很多学者对箱体的二维温度场与流场进行了研究。凌长明，陶文铨建立了冰箱非稳态自然对流换热的二维计算模型，分析三个互不相通空腔内的自然对流换热情况，计算了周期性非稳态工况下温度场和流场的分布[4]。
上海交通大学的丁国良等使用有限元软件FIDA7.6，模拟了冰箱二维稳态自然对流空气流场分布情况，并研究了内部热负荷、隔板以及蒸发器和门之间的间距、内部隔板的导热系数对箱体内部温度场和流场分布的影响[5]。
周湘江等针对冰箱周期性非稳态传热的特点，把导热、对流和辐射作为整体耦合求解，使用二维模型进行模拟仿真得出箱体内传热规律以及辐射的影响[6]。
吴小华采用FLUENT软件建立了冰箱三维几何模型，分别模拟网状和平板搁物架时箱体内温度场和流场的分布，得出了较为合理的结构，同时提出了优化的方向[7]。
于兵应用PHOENICS软件模拟了间冷式冰箱箱内温度场与流场的分布情况，在建立物理模型过程中，对其做出一些简化；通过分析得知，计算结果与实验结果基本吻合，并且可以忽略浮升力的影响[8]。
苏秀平等应用FLUENT软件计算了间冷式冰箱风扇区域流场，采用二阶k-ε湍流模型，风扇叶片区域使用多坐标模型，蒸发器区域使用多空介质模型，提出了风扇盖板的改进措施[9]。
俞炳丰等使用PIV技术测量了间冷式冰箱内部流场，并将采用k-ε紊流模型计算模拟得到的结果与测试结果对比分析，得出计算模型的正确，为箱体的结构优化提供了依据[10]。
2 冷藏箱研究中存在的问题
很多学者为了问题的简化，大都研究了冷藏箱的二维稳态情况，而很少考虑三维非稳态的情况，忽略了整体的影响，因此计算结果有片面性，计算结果的误差较大。数值计算中没有考虑箱体内的辐射换热，而因为箱体各壁面的温差较大，辐射换热对箱体内温度场与流场的分布有着重要的影响[11]。同时，门缝结构对箱体内温度场与流场分布的影响，无论是通过实验的方法还是数值计算方法，存在一定的难度，因此需要更有效的解法。
3 总结
从以上分析可以看出，CFD已在冷藏箱箱体结构设计以及优化发挥了重要的作用，但是还有很多值得我们去探索的地方：
（1）CFD在冷藏箱箱体温度场与流场的计算模拟中已经应用的非常成熟，但是很多都是为了简化问题模拟了箱体的二维稳态情况，对三维非稳态情况涉足较少，而下一步的工作需要考虑所有的箱体壁面辐射换热问题。
（2）模拟计算中通常将箱体作为一个空腔结构，忽略了搁架以及放置物品对冷藏箱内部换热和流场的影响，但实际中搁架及其上的物品对冷藏箱内部温度场和流场分布较大的影响。
（3）蒸发器的位置以及形式对箱体内换热的影响，通常考虑的是它对冷藏箱能耗以及本身换热性能的研究，而忽略了它的位置和形式对箱体的温度场与流场分布的影响。
参考文献
［1］乔洪涛,卢智利,丁国良,张春路.CFD在冰箱结构优化中的应用现状[J].低温工程,2003(4).
［2］Cezar O.R. Negrao, Christian J.L. Hermes. Energy and cost savings in household refrigerating appliances: A simulation-based design approach[J]. Applied Energy,2011,88:3051-3060.
［3］C.Conceicao Antonio, C.F. Afonso. Air temperature fields inside refrigeration cabins: A comparison of results from CFD and ANN mode lling[J]. Applied Thermal Engineering,2011,31:1244-1251.
［4］凌长明,陶文铨．冰箱内非稳态自然对流的二维数值模拟[J]．西安交通大学学报,1995,29(10):35-41.
［5］Guo Liang Ding, Hong Tao Qiao, Zhi Li Lu. Ways to improve thermal uniformity inside a refrigerator[J]. Applied Thermal Engineering,2004,24:1827-1840.
［6］周湘江,凌长明．家用冰箱内传热的非稳态数值研究[J]．中南工学院科技通讯,1999,15(1):21-24.
［7］吴小华,张璟,宋春节．冰箱冷藏室温度场和流场的仿真与优化[J]．北京石油化工学院学报,2006,1(3)．
［8］于兵,童灵,阙雄才.间冷冰箱三维气固耦合传热与流动数值研究[J].1998(7).
［9］苏秀平,陈江平,陈芝久,等．间冷式冰箱风扇区域流场的数值模拟和优化[J]．上海交通大学学报,2003,37(7):1133-1136.
［10］俞炳丰,费继友,孟祥兆,等．间冷式冰箱冷冻室内流场的PIV测试和计算模拟[J]．制冷学报,2003(2):32-36.
数值模拟范文5
关键词：基坑；土钉支护；FLAC3D；数值模拟
中图分类号：TU473
文献标识码：B
文章编号：1008-0422(2006)04-0118-03
收稿日期：2006-04-25
作者简介：梁琼(1967-)，男(汉族)，湖南涟源人，工程师，主要从事建筑结构设计。
1 前言
土钉支护是近几十年发展起来用于加固和增强边坡或开挖土体稳定的一种支挡技术，它是在原位土中自上而下设置细长、较为密集的金属杆件(土钉)、与土坡表面构筑的钢丝网喷射混凝土面层及被加固土体共同作用，形成一个自稳的和能支挡墙后土体的支挡结构。
土钉支护结构以其施工便捷、经济和支护效果好等诸多优点，在当前的基坑工程中应用得越来越广泛。然而目前对于土钉体、面层和被加固土体等的联合作用机理以及土钉支护对基坑边坡土移和应力的影。向还缺乏较深入的认识，由于深基坑开挖的现场难以对土钉支护参数进行系统的研究分析，因此采用数值模拟的方法对其进行研究分析，对基坑的支护设计具有较大的指导意义。
2 FLAC3D简介
FLAC全称是FastLagrangianAnalysisOfContlnua，是由美国Itasca咨询公司于1986年开发的显式有限差分法程序，90年代中期，Itasca公司在原有的二维分析软件的基础上开发出了FLAC3D。该程序采用拉格朗日差分公式来处理有限变形问题，计算过程中允许材料发生屈服及流变，适合于解决岩土工程中经常遇到的大变形，是一种理想的岩土工程计算软件。
3 工程概况
某工程为1栋11层写字楼，高度51m，框剪结构，设?层地下室，阀板基础，基础埋深7．Om，需开挖的基坑深度为7．Om。该基坑的土层情况及力学性质见表1：
4 FLAC数值模拟及计算分析
4．1土钉支护方案
本工程需开挖的基坑深度为7．Om，基坑围护结构采用土钉挂网喷砼支护体系，不设内支撑。共设6排土钉，每排土钉长8．4m，约为1．2H(H为基坑深度)。第一排土钉水平布置，其余土钉倾角取10’。土钉水平和垂直间距均为1．2m。成孔直径110mm。挡墙厚度为100mm。
基坑在开挖前进行了全面降水处理，土体沉降基本稳定，所以在进行数值模拟时不考虑地下水的作用。
整个基坑开挖和土钉支护的模拟分以下七步：
(1)模型的初始应力场的模拟；
(2)开挖到距离地面0．3m，打入土钉，保持土体稳定，注浆，喷射面层；
(3)开挖到距离地面1．5m，打入土钉，注浆，喷射面层；
(4)开挖到距离地面2．7m，打入土钉，注浆，喷射面层；
(5)开挖到距离地面3．9m，打入土钉，注浆，喷射面层；
(6)开挖到距离地面5．1m，打A-f-钉，注浆，喷射面层；
(7)开挖到距离地面6．3m，打入土钉，注浆，喷射面层；
(8)开挖到距离地面7．0m，喷射面层。
4．2计算模型及网格划分
土钉是在原位土体中插入截面较小的加筋体，具有三维作用。如果既考虑基坑的整体三维作用，又考虑土钉的三维局部作用，计算量和存贮量都太大，为此，本文对采用的三维模型作如下考虑：土钉一般是成排成列重复布置的，一列竖向土钉的加固范围为通过两侧水平间距中点的竖直截面间的土体，厚度等于水平间距，这块土体相对于通过土钉的竖直截面又是对称的，为减少计算量，取沿着基坑长度方向1．2m进行三维建模，如图1。
模型其它两个方向的尺寸根据地质条件和基坑开挖深度，经试算确定取基坑深度的2．5倍左右。模型四周与相邻的土体间有一定的相互约束，为简化分析，在模型的边界处施加适当的约束条件，底面为固定铰支，四个侧面分别为滚动支座，竖向方向没有约束，可以自由滑动。
4．3计算参数的选取
土层材料性质参数按表1―1选取。土钉支护结构中，土钉弹性模量值为El=2x10”，挡墙弹性模量EM为2．5x1010。泊松比rm=0．15，挡墙厚度为100mm。注浆体粘聚力都取2．2 x104Pa／m，土钉钢筋采用25螺纹钢，断面面积为4．9x10-4m2，屈服强度为200MPa。
4．4基坑分步开挖及土钉支护模拟
本例中基坑边坡选择土的应力应变模式为摩尔―库仑模式，这种模式概念明确、计算简便，对一般工程问题都有令人满意的精度；并且假设其为大变形。基坑的开挖和支护是逐步完成的，而支护体系的应力和变形往往又和施工过程紧密相关，所以为了较为真实可靠地分析支护体系的应力和变形，对支护工程施工的模拟是十分必要的。土钉支护用于基坑或边坡土体开挖时，从上到下分步修建，边开挖边支护，模拟这一施工过程的有限元网格划分如图2所示。
4．5结果分析
(1)基坑边坡位移场分析
随着开挖深度的不断增加，基坑水平位移与地面沉降逐渐增大；离基坑顶点越近水平位移越大，随着开挖的进行，水平位移等值线越来越陡，土体整体趋向于基坑内侧移动；在每步开挖完成后在坡顶处水平位移一直最大，所以在施工该部位时，应该特别注意，当开挖深度较大时，最好将其改为预应力土层锚杆。
基坑边坡的位移(变形)，特别是土钉边坡坡顶的水平位移值是控制基坑稳定性和满足环境要求的最重要和直接的指标。较合理的做法是选择基坑开挖终了时基坑边坡坡顶最大水平位移、坡顶最大沉降值、坡中点水平位移作为评价基坑稳定性和满足环境要求的指标。表2列出了模拟开挖各步时基坑顶点的水平位移值和竖向沉降值，并与实测值加以对比。
就本文实例来看，实际位移略小于数值模拟的结果。总体而言，实际发生的数值与FLAC数值模拟的结果基本吻合，说明采用FLAC可以对土钉支护的变形和土钉受力等做出较好的预测。
(2)基坑应力分析
基坑在开挖过程中，随着深度的不断增加，土体被扰动的程度不断加大，更大范围的土体受到影响，逐渐出现塑性屈服，受力状态也从初期的弹性应力状态转变为塑性应力状态。而随着支护结构的依次加入发挥作用，调整了土体的受力状况。土体在经过一段时间的应力重分布，达到稳定后，又恢复到弹性应力状态。
从图3中可以看出在土钉范围内土体的
应力水平较土钉范围外的土体应力水平迅速增加，而在坡脚处有应力集中现象，土钉支护后的土体可以看作是各向异性的复合材料，由于土钉的弹性模量远大于土体的弹性模量，在土体发生变形的情况下，土体的应变大于土钉的应变，土钉与土体的相互作用将在界面上产生摩阻力，这个摩阻力使土钉受拉，在土钉中产生拉力，使土体的侧向应力增大，在一定程度上弥补了由于土体开挖卸载引起的(侧向应力)的减小，即相当于在土体原有的应力基础上增加了一个Q3，使得土体的强度提高。
(3)土钉的轴力分析
假设在一定小的范围内(L)，土钉的应力分布是均匀的，其应力值大小相等。本例基于以上假设，在采用FLAC程序进行分析时，先将每根土钉设置为杆单元(ca-ble)，然后再以L＝L／20为单位划分成若干等距小单元，由此求得每根土钉的应力分布。土钉置入现场土体后，随着向下开挖或土体徐变即产生土体变形，于是通过土体与土钉之间的界面粘结力使土钉受拉。只要土体发生微小的变形，就可使土钉受力。土钉的拉力沿其长度变化，最大拉力部位随着基坑向下开挖，逐渐向靠近面层的端部转移，一般发生在土体的可能失稳破坏面上。当土钉长度较短时，土体破坏面可能移出上部土钉之外，这些钉中的最大拉力一般发生在钉长中部。不同深度位置上的土钉，其受到的最大拉力有很大差别，顶部和底部的土钉受力较小，靠近中间部位的土钉受力较大。而对于同一根土钉而言，其内力并不均匀，而是中部大、两头小，呈枣核状。如图4、图5所示。
图6给出的是每根土钉的最大应力位置。
由以上三图可见，土钉应力分布是很有特点的。首先它存在一个峰值，此峰值恰好是内部潜在的破裂面；其次是应力向峰值两侧递减，在土钉两端递减为零。由此可得，土钉长度只需满足设计要求即可，不应盲目增加土钉长度；在土钉端部，应力也递减很快，到达面层时已很小，所以土钉支护中面层虽然承受一定的内力，但受力不大，面层不是主要的受力构件。受拉应力值最大的土钉位于基坑的中下部，所以在开挖支护过程中应加倍注意中下部土钉的及时到位，密切注意基坑变形。
5 结论
(1)土钉支护基坑边坡可以有效地减小土体的侧向位移，提高承载力。土钉墙在受荷过程中由于土钉的作用推迟了塑性变形的发展阶段，呈现出渐进性变形与开裂破坏并存且逐步扩展的现象，不会发生像土边坡那样的突发性滑塌。
(2)应用FLAC程序对基坑土钉支护进行数值模拟是可行的，从分析得到的基坑边坡位移场、应力场和土钉的应力分布规律与实测结果基本吻合。
(3)由FLAC程序计算得到的单根土钉的应力分布趋势为中间大两端逐渐变小，最大拉力部位随着基坑向下开挖，逐渐向靠近面层的端部转移，最大拉力一般发生在钉长中部。不同深度位置上的土钉，其受到的最大拉力有很大差别，顶部和底部的土钉受力较小，靠近基坑边坡中间部位的土钉受力较大。
(4)为控制基坑变形，应该适当增加土钉长度，但不能过长。因为土钉长度过长的话，控制变形的效果并不明显，而且土钉拉力增加不多，强度得不到发挥，造成浪费。
参考文献：
[1]闫莫明，徐祯祥，苏自约．岩土锚固技术手册[M]．北京：人民交通出版社，2004．
[21乔俊宇，徐国元．深基坑土钉支护数值模拟[儿江西有色金属，2005(12)．
[31杨立超．土钉支护及其FLAC数值模拟[J].建筑技术开发，2003(9)．
数值模拟范文6
关键字：内域；数值计算；有限元
中图分类号：TU315.3
1.引言
虽然目前技术和计算设备发展十分迅速，计算能力不断提高，一些大型通用有限元软件已经具备十分强大的分析功能，它解决了地震反应的许多实际工程结构分析的问题。但对于一些大型复杂体系而言，空间离散的自由度数目非常庞大，数值稳定性的限制要求时间离散的步距也不能过大。这样，开展结构地震反应分析时所需要完成的时空四维数值计算的工作量将变的很大。在工程设计中，需要分析各种工况下结构的地震反应行为，对比不同的设计方案并做出优化决策，从而要求在设计期内多次完成结构的地震反应计算。在数值精度的基础上，保证系统的稳定和提高结构地震反应分析方法的核心计算效率。研究高效率的大型复杂体系地震反应数值分析方法，减少计算费用仍然是非常现实的考虑，它具有重要理论意义和实用价值。
地震反应分析的数学物理模型就是波动方程。波动数值模拟包含两个部分，一是对人工边界的数值模拟，二是对内域的数值模拟。这里仅讨论对内域的运动节点的数值模拟问题。
现有的内域波动的数值模拟方法，按能量等效式，可划分为三类。一类是空间有限元方法。所谓空间有限元方法是指的空间用有限元法进行离散，而使用一个直接的方式离散时间。第二类是时空有限元方法。这里是指对时间使用和对空间一样的有限元法方式离散。时空有限元的时空域是空间有限元的空间域在时间域上增加了一个时间维度，区别只是处理时间域的先后上，空间有限元是先处理空间问题，然后处理时间问题，而时空有限元是同时处理时间和空间的问题。第三类是微分求积方法，是直接的方式对空间和时间的离散。下面分别阐述这三类方法的应用发展过程――
2.空间有限元方法
在1956年，有限元的概念首次被Turner等人提出，最早应用于弹性力学平面应力问题上。1963年，Besseling、Melosh和Jones等人发现有限元法和基于变分原理的里兹法是等效的。有限元法在处理连续介质问题上比普通里兹法更有优势。随后几十年，在解决复杂工程问题上，有限元法得到广泛的应用。
波动方程是时空耦合的，基于广义Hamilton原理的波动有限元方法通常也是时空解耦的数值过程。传统有限元方法的离散过程通常包含两步，先进行空间离散，将微分方程转化为常微分方程，然后对时间进行离散，即在时域对常微分方程进行数值积分。由于时空耦合的数值过程包含过多的自由度，求解这类方程在实际工程中很难实现，建立时空解耦的波动数值分析方法是这方面重要的工作。最直接的做法是实现空间及时间域的解耦，通常是只建立空间的有限元离散方法，而时间采用直接的假设，最常用的是采用逐步积分的方式进行离散。
逐步积分法简单来讲就是把最终速度和位移由它们的初始值和一个积分表达式来表示。加速度历程的积分决定速度的变化，速度的积分决定位移的变化。换句话说，加速度控制了速度的变化，因而可以由这一步向前获得下一时间步。解答这类问题，第一步先考虑时间步内的加速度问题，假设加速度是如何变化的，依据加速度和位移的关系，得到关于时间步的递推公式。所谓的Euler-Gauss法就是假设在时间间隔内加速度为常数。而Newmark[1]法是加入系数从而可以改变初始和最终加速度的权重从而得到加速度的一种方法。Wilson- [2]方法是假定在时间步距内加速度为线性加速度的一种数值方法，用内插公式得到体系在下一刻的运动。α方法[3]是在Newmark方法的基础上，通过修改结构动力方程的时间离散形式得到的。Chung 和 Hulbert 发展了一种无条件稳定的隐式广义α法，它由三个参数控制数值损耗。Runge-Kutta[4， 5]方法是在一个时间步距中内插若干计算点，利用这些计算点上函数值线性组合来代替函数的泰勒展开中的高阶导数，从而提高精度阶。不同于两步信息预测，线性多步法[6]发展了多步信息来预测下一步，从而获得了更高的精度。
逐步积分是最主要的时域积分方法，而它最常规的做法是差分法。时域有限差分法（Finite Difference Method， 简称FDM），是地震波传播模拟最广泛被使用的一类方法[kelly―Marfurt，1990][7-10]。有限元差分是将微分方程中的微分项用相应的差商代替，从而将微分方程转化为代数形式的差分方程。
由微商和差商的定义可以知道，微分的有限形式是差分，而导数的有限形式是差商。而微分和导数是以极限形式表示。数值计算方法导数可以用差商的自变量趋近于零来代替，换句话说，位移对时间的求导可以用有限差分的方式得到，位移的一阶导数是速度，二阶导数加速度。当世界步距为等步长，得到中心差分。差商代替微分方程中的导数，就可以得到微分方程的有限差分形式。。较之传统有限元法，虽然在定义几何结构上不够灵活，但时域有限差分法具有显式计算的优势，计算效率高，计算精度高于显示有限元法。这些方法的缺点是精度不高，只有一阶或者二阶精度，难以模拟高频问题，这类无法避免算法阻尼，从而形成较大的误差。
为了克服低精度的问题，很多高精度的数值积分方法相继提出。不仅仅四阶，五阶精度，甚至更高精度的数值积分方法处在发展之中。Golley[11]为了得到三阶精度的格式，对时间域采用高斯点作为配置点。在哈密顿变分原理的基础上，Riff和Barch[12] 采用3次多项式构造插值函数，得到了有条件稳定四阶精度数值方法。Argyris[13]等在前人的基础上，采用Hermitian插值，得到无条件稳定四阶精度数值方法。Kujawski和Gallagher[14]从另外一个角度，利用广义最小二乘法，在无阻尼结构中，构造了一种四阶精度的无条件稳定积分格式。Tarnow和Simo[15]利用二阶精度算法的结果，在此基础上缩短时间间隔，构造了一种四阶精度算法。钟万勰[16-19]在1994年提出了精细时程积分法。在保守系统下，积分结果保持系统守恒量不变。1995年在之前工作之上，钟万勰提出了子域精细时程积分法，提高了计算效率，使工作量大大减小，存储量大幅减小，为精细法应用提供了基础。2000年，顾元宪[20]提出了增维精细积分法，改进了矩阵的运算，提高了精度。但局限条件较多。2004年，汪梦甫[21]利用精细积分方法基本原理，采用高斯积分方法，建立了更新精细积分方法，这种精细方法适应度高，为得到了广泛的应用提供了条件，并且提高了精度。理论上汪梦甫分析了精细方法得到任意精度的可能。2009年，富明慧[22-25]在汪梦甫研究的基础上提出了高效高精度广义精细积分法。
这类方法的困难在于不容易构造精度较高的时间离散模式，并且空间有限元在时域上每个时间步逐步推进，因而会产生误差累积[26-32]。
3.时空有限元方法
时空有限元最早由J.TOden，I.Fried和J.H.Argyris等人提出。利用哈密顿原理建立关于时间边界的变分原理。几十年来，在各个领域得到全面的发展。在波动问题，动力问题，流体问题等方向得到广泛的应用。
传统的数值方法假定时间和空间是相互独立的，这样的假定广泛应用在实践当中并且在数学上很好理解。但是，使用上述技术同时也导致了数学上的困难。因为有用的信息可能通过速度在结构传播，而没有分离的时空格式能规避这种类型的数值困难。这就要求对时间的离散和对空间离散一样使用有限元。例如，当结构是三维时，这种格式需要四维的维度来表示。从而需要对时间和对空间使用相同的离散方式。空间有限元对时间和空间分别离散，空间节点形成的网格只能处在同一时间下，形成的是规则网格。时空有限元同时对时间和空间进行离散，理论上可以把网格划成任意形状，不必考虑相同的时间值，可以灵活的划分。有限元方法区别于其它方法在于它利用能量等效原理将偏微分方程进行积分。得到方程的弱形式。恰当的变分形式是有限元是否能够利用能量等效的关键。而时空有限元能否成功取决于能否找到对应的变分原理。
R.Riff和M.Baruch等人建立了一种能同时求解所有变量的时间有限元，这种有限元格式借鉴了空间有限元，把整个时间域看作是空间的延续，从而能够求解不同时刻的变量。冯康是提出 罗恩在冯康的基础上，完善了哈密顿型变分原理，发展了非传统哈密顿变分原理。罗恩对比等价的正则方程和相空间变分原理，认为即使是等价的，但由于形式的不同，产生的算法并不一定会有相同的效果。其结果就是相空间变分原理计算效率更高，也更接近物理问题的本质特征。沿着这个思路，罗恩建立了非传统相空间哈密顿变分原理。刘世奎对哈密顿原理进行了推广，构建时间边界条件，建立了有两个参数的广义哈密顿变分原理，完成了从单一变量向多种变量的转变。卓家寿等人在哈密顿体系下，推导了几种变分原理的等价形式。罗恩在此基础上推导了类变量广义哈密顿变分原理，它包含了所有的条件。基于这种变分原理提出在时间子域上进行五次样条插值的时间子域法。 2007年钟万勰[38]发现时空有限元可以更高效的解决边界问题和多尺度刚度问题。2013年朱宝[39]在钟万勰的基础上进一步研究多尺度和边界问题，对其稳定性进行了研究。
获得高阶精度，只需要提高多项式基函数的次数，从理论上来说，时空有限元可以获得任意精度。空间域增加时间域之后，单元的几何性质简单，避免了空间有限元的复杂边界。让传统的有限元得到更广的应用。是一种有很大发展空间的数值方法。
4.微分求积方法
Bellman和Casti[40]在1971年提出微分求积法的基本原理。此后，微分求积法因为原理简单，广泛应用在工程问题中，微分求积法得到快速发展。
微分求积法即DQM方法，本质上来说，函数和它的导数在给点节点的值用全部节点的函数值乘以系数并求和来代替。从而让微分方程可以转化成关于节点函数值的一组代数方程组。由此可知，DQM是一种数值技术，它通常被用来解决初值和边界问题。从本质上来说，DQM是特殊的一种加权残值法，而且是高阶的有限差分法。DQM方法相对有限元方法而言，并不需要变分原理就可以求解微分方法。
从微分积分法的原理出发，可以发现影响数值精度主要由两个因素构成。一方面是权系数的值，另一个方面是选取合适的网格离散点。其中网格离散点的选取和假设试函数的模式可以确定权系数，从另一个角度来说，试函数的假设和网格点的选取是决定精度的关键因素。从而研究人员也沿着这个思路对微分求积方法进行了探索。利用多项式是有限元试函数选取的基本思路。Bert和Wang[41]为试函数来求权系数，此时构成线性方程组的系数矩阵是勒得蒙矩阵。但由于当离散点数目增多以后，勒得蒙矩阵会出现病态。所以出现很大的误差。后来。Quan[42， 43]用采用了Lagrange插值，得到了微分积分法一阶和二阶精度的公式。Bert和Striz[41]建立了HDQ方法，采用用不同于多项式的谐函数作为试函数，开阔了研究思路。由上可知，试函数的选取并不是单一的，可以从多个角度来选取。不仅是谐函数或者多项式，甚至是指数函数，对数函数等初始函数都可以作为试函数进行研究。根据需要选择混合初始函数来得到试函数是值得探索的方向。网格点的选择方面，研究发现一些问题对节点选取是很敏感的，等距网点因为使用方便而被先采用，但是结果发现得到的解不够理想。其实真实的状况让均匀网格模拟显得不够合理，发展非均匀网格更可能得到高精度方法。 Bellman[40]就用勒让得多项式零点进行了研究，发现用其作为网格点提高了精度。在这启发之下，Quan等研究了切比雪夫多项式零点，并且用之与其它正交多项式作了对比，发现切比雪夫多项式零点更有优势。此外，微分求积的研究的方向是更加具体的研究，Bellman[44]用微分求积法应用到初值非线性偏微分方法得到高效精确的解法。 在多维问题上，微分积分法也得到了应用，Civan[45]得到了多维积分微分方程。Bert[46]首先将这种方法用到结构力学问题的求解。2001年在DQM法则的基础上，Fung[47， 48]建立了一种不同于边值问题的动力微分方程解法，解决了初值问题的动力微分方程。并研究时间网点选择方式对数值精度和稳定性的影响，提出了一种高精度的时间网格离散方法。
微分求积法虽然发展的历史比较短，但是由于这种方法原理简单，精度高，计算效率高，处理数据方便等优点。在工程领域有广泛的应用。
5.结语
有限元方法一直在数值模拟中很占有重要地位，这种思想在理论研究和实际应用中发挥着很重要的作用。以有限元原理为基础，发展的新方法让数值计算展现出新的活力。但是数值模拟是一门深奥的学问，在理论上和实际应用中还有很多不完善的地方，需要克服的缺点还有很多，本文作者仅仅就自己所涉及的研究领域做了一些简单的论述。
参 考 文 献
[1]. Newmark， N.M. A method of computation for structural dynamics[C]. in Proc. ASCE. 1959.
[2]. Wilson E L. A computer program for the dynamic stress analysis of underground structures[R]. DTIC Document， 1968.
[3]. Trefethen L N. Finite difference and spectral methods for ordinary and partial differential equations[M]. Cornell University [Department of Computer Science and Center for Applied Mathematics]， 1996.
[4]. 易大义，陈道琦. 数值分析引论[M].杭州：浙江大学出版社，1998.
[5]. 胡健伟，汤怀民. 微分方程数值方法[M].北京：科学出版社，1999.
[6]. Chew W C. Waves and fields in inhomogeneous media[M]. IEEE press New York， 1995.
[7].Wilson E L， Bathe K. Numerical methods in finite element analysis[M]. Prentice-Hall， 1976.
[8]. Clough R W， Penzien J. Dynamics of structures[R]. 1975.
[9]. Wilson E L， Farhoomand I， Bathe K J. Nonlinear dynamic analysis of complex structures[J]. Earthquake Engineering & Structural Dynamics. 1972， 1（3）： 241-252.
[10]. 刘晶波，杜修力.结构动力学[M]. 北京：机械工业出版社，2005.
[11]. Golley B W. A TIMESTEPPING PROCEDURE FOR STRUCTURAL DYNAMICS USING GAUSS POINT COLLOCATION[J]. International Journal for Numerical Methods in Engineering. 1996， 39（23）： 3985-3998.
[12]. Riff R， Baruch M. Time finite element discretization of Hamilton's law of varying action[J]. AIAA journal. 1984， 22（9）： 1310-1318.
[13]. Argyris J H， Vaz L E， Willam K J. Higher order methods for transient diffusion analysis[J]. Computer Methods in Applied Mechanics and Engineering. 1977， 12（2）： 243-278.
[14]. Kujawski J， Gallagher R H. A generalized leastsquares family of algorithms for transient dynamic analysis[J]. Earthquake engineering & structural dynamics. 1989， 18（4）： 539-550.
[15]. Tarnow N， Simo J C. How to render second order accurate time-stepping algorithms fourth order accurate while retaining the stability and conservation properties[J]. Computer Methods in Applied Mechanics and Engineering. 1994， 115（3）： 233-252.
[16]. 钟万勰. 结构动力方程的精细时程积分法[J].大连理工大学学报.1994，34（2）：131-136.
[17]. 钟万勰. 计算结构力学与最优控制[M].大连：大连理工大学出版社，1993.
[18].钟万勰. 应用力学的辛数学方法[M].北京：高等教育出版社，2006.
[19]. Zhong W X， Williams F W. A precise time step integration method[J]. Proceedings of the Institution of Mechanical Engineers， Part C： Journal of Mechanical Engineering Science. 1994， 208（6）： 427-430.
[20]. 顾元宪，陈飚松，张洪武.结构动力方程的增维精细积分法[J].力学学报.2000，32（4）.
[21].汪梦甫，周锡元. 结构动力方程的更新精细积分方法[J].力学学报.2004，36（2）：191-195.
[22]. 富明慧，刘祚秋. 固体力学中的变分差分方法[J].中山大学学报：自然科学版.2001，40（2）：13-15.
[23]. 富明慧，林敬华，刘祚秋.结构动力分析的广义精细积分法[Z].第九届全国振动理论及应用学术会议论文摘要集.2007.
[24]. 富明慧，林敬华. 精细积分法在非线性动力学问题中的应用[J].中山大学学报：自然科学版.2008，47（3）： 1-5.
[25]. 富明慧，廖子菊，刘祚秋. 结构动力方程的样条精细积分法[J].计算力学学报.2009， 26（3）： 379-384.
[26]. 裘春航，吕和祥，钟万勰. 求解非线性动力学方程的分段直接积分法[J].力学学报.2002， 34（3）： 369-378.
[27]. 郑兆昌，沈松，苏志霄.非线性动力学常微分方程组高精度数值积分方法[J].力学学报.2003，35（3）： 284-295.
[28]. 胡海昌.弹性力学的变分原理及其应用[M].北京：科学出版社，1981.
[29]. 赵秋玲.非线性动力学方程的精细积分算法[J].力学与实践.1998，20（6）：24-26.
[30]. 王金东，高鹏. 波动方程的精细逐步积分法[J].力学季刊.2000，21（3）：316-321.
[31]. 张洪武.关于动力分析精细积分算法精度的讨 ）[J]. 力学学报. 2001，33.
[32]. 梁立孚胡海昌.一般力学中三类变量的广义变分原理[J].中国科学： A 辑. 2000，30（12）：1130-1135.
[33]. 罗恩.几何非线性弹性动力学中广义 Hamilton 型拟变分原理[J].中山大学学报（自然科学版）.1990， 29（2）：15-19.
[34]. 罗恩，黄伟江.相空间非传统 Hamilton 型变分原理与辛算法[J].中国科学：A辑.2002，32（12）： 1119-1126.
[35].罗恩，梁立孚，李纬华.分析力学的非传统 Hamilton 型变分原理[J].中国科学：G 辑.2007，36（6）： 633-643.
[36]. 梁立孚，罗恩，冯晓九.分析力学初值问题的一种变分原理形式[J].力学学报. 2007， 23（1）： 106-111.
[37].罗恩，朱慧坚，黄伟江.Hamilton 弹性动力学及其辛算法[J].中山大学学报（自然科学版）. 2003，42（5）： 131-132.
[38]. 钟万勰，高强. 时间-空间混和有限元[J].动力学与控制学报，2007，5（1）：1-7.
[39]. 朱宝，高强，钟万勰. 三维时间-空间混和有限元[J].计算力学学报，2013，30（3）：331-335.
[40]. Bellman R， Casti J. Differential quadrature and long-term integration[J]. Journal of Mathematical Analysis and Applications. 1971， 34（2）： 235-238.
[41]. Bert C W， Xinwei W， Striz A G. Differential quadrature for static and free vibration analyses of anisotropic plates[J]. International Journal of Solids and Structures. 1993， 30（13）： 1737-1744.
[42]. Quan J R， Chang C T. New insights in solving distributed system equations by the quadrature method―I. Analysis[J]. Computers & chemical engineering. 1989， 13（7）： 779-788.
[43]. Quan J R， Chang C. New insights in solving distributed system equations by the quadrature method―II. Numerical experiments[J]. Computers & chemical engineering. 1989， 13（9）： 1017-1024.
[44]. Bellman R， Kashef B G， Casti J. Differential quadrature： a technique for the rapid solution of nonlinear partial differential equations[J]. Journal of Computational Physics. 1972， 10（1）： 40-52.
[45]. Civan F， Sliepcevich C M. Differential quadrature for multi-dimensional problems[J]. Journal of Mathematical Analysis and Applications. 1984， 101（2）： 423-443.
[46]. Bert C， Jang S， Striz A. New methods for analyzing vibration of structural components[Z]. Structures， Structural Dynamics and Materials Conference， 28 th， Monterey， CA， Apr. 6-8， 1987 and AIAA Dynamics Specialists Conference， Monterey， CA. 1987，936-943.
[47]. Fung T C. Solving initial value problems by differential quadrature method―part 1： firstorder equations[J]. International Journal for Numerical Methods in Engineering. 2001， 50（6）： 1411-1427.
[48]. Fung T C. Solving initial value problems by differential quadrature method―part 2： secondand higherorder equations[J]. International Journal for Numerical Methods in Engineering. 2001， 50（6）： 1429-1454.
