时间：2023-06-01 08:51:34
开篇：写作不仅是一种记录，更是一种创造，它让我们能够捕捉那些稍纵即逝的灵感，将它们永久地定格在纸上。下面是小编精心整理的12篇数值方法，希望这些内容能成为您创作过程中的良师益友，陪伴您不断探索和进步。
关键词数值模拟；数值流形；边坡；物理网格
中图分类号：C35文献标识码： A
1 引言
数值流形方法（Numerical Manifold Method）是石根华先生基于有限覆盖技术通过研究非连续变形分析法的基础上提出的新型数值方法。
数值流形方法应用现代数学“流形”的覆盖技术，将连续体的有限单元法、不连续体的非连续变形分析法和解析方法统一起来，是一种更高层次的计算方法[1][2]。其基本思想是将许多只占局部区域的、相互重叠的有限覆盖连接成一个覆盖系统，用这个覆盖系统去覆盖求解的区域，在各个有限覆盖上构造覆盖函数形成总体函数去逼近求解域的真实场函数。
对于连续的材料区域，数值流形方法采用分片可微分的在材料的局部区域定义的有限个覆盖位移函数来描述材料的总移场，这些覆盖之间是部分重叠的，覆盖位移函数在重叠区域通过权函数进行叠加，能够反映连续材料位移场的连续性。
对于岩土工程重的不连续材料区域，各区域之间覆盖函数没有相互的重叠，使其所构成的总移函数在不连续的材料之间也是不连续的。由此，数值流形方法可以达到统一解决连续和非连续问题的目的，采用相应的覆盖函数能够较精确地计算材料体的位移、变形和应力，不连续面两侧材料的相对位移及材料的大变形、大位移等[3][4]。
2 数值流形方法的基本原理
数值流形法是一种基于有限覆盖技术的数值方法。有限覆盖包括数学覆盖和物理覆盖两层含义，数学覆盖是数值流形法中的基本网格，这些数学覆盖被物理边界切割而形成物理覆盖，物理覆盖的重叠区域即形成流形单元。将物理覆盖上的位移函数结合起来形成计算域上的全域位移近似函数，似函数就是形成此单元的若干个相互重叠的物理覆盖上的近似函数的加权平均，然后即可利用最小势能原理形成整体平衡方程。
其基本流程可做如下描述（如图2.1）：
2.1覆盖函数
数值流形方法以物理覆盖为覆盖函数的构造区域，而以流形单元作为覆盖函数的插值区域，由此来形成求解域上的整体函数。
对于平面问题，规定物理覆盖上的覆盖函数为：
其中，为完全多项式的基本级数函数；对于完全0、1、2阶覆盖函数，m分别为1、3、6。
2.2 总移函数
各物理覆盖的覆盖函数在它们的交集――流形单元上，通过权函数对其进行加权平均构成流形单元的总移函数。设分析域内有个物理覆盖，每个流形单元有个物理覆盖(个流形单元)，每个物理覆盖有个未知系数(广义自由度)，则物理覆盖函数为：
2.3 总体平衡方程
数值流形方法较全面地考虑了参与体系平衡的各种平衡力项。对于不同的问
题，参与平衡方程的各势能源不同，总的来说主要有：①应变能势能；②初应力势能；③点荷载势能；④体荷载势能；⑤惯性力势能；⑥用于不连续变形分析的接触弹簧力势能和摩擦力势能等。
在得到流形单元上的总移函数后，就可建立弹性力学边值问题中的能量泛函表达式。系统的总势能为
其中: 为单元应变能； 为初应力势能； 为点荷载势能； 为体荷载势能； 为惯性力势能； 为约束形成的势能。因为每个物理覆盖有个未知系数，由式给出的系数矩阵的子矩阵是一个的矩阵。和是的子矩阵，表示物理覆盖的个未知系数(节点的自由度)。
3 边坡工程的数值流形计算分析
设某边坡的一剖面如图3.1所示， ，数学覆盖采用有限元三角形网格。利用数值流形方法模拟边坡的滑动，对边坡进行稳定性分析。具体参数如下：
表3.1边坡算例计算条件
计算条件 计算参数 计算条件 计算参数
边坡高度 10m 边坡底边长度 10m
弹性模量 28GPa 泊松比 0.27
摩擦角 14o 接触面粘结力 0.01MPa
容重 26KN/m3
图3.1某边坡裂隙剖面图
3.1覆盖分析
利用有限元网格对物理网格覆盖，生成231个流形单元，如图3.2，流形方法中的单元通常具有不规则的形状，基于三角形有限元网格所定义流形单元均具有三个物理覆盖，即三个物理覆盖的交定义了一个流形单元，这三个覆盖可看作流形单元的3个节点.在无不连续面域内，相邻流形单元沿公共边有相同节点；有不连续面或材料界面分离的两个流形单元有不同的节点。
图3.2边坡有限元网格覆盖
3.2变形分析
考虑到可能影响边坡稳定的因素，如降雨使土体重度增加、受河流冲刷、地下水活动、地震及人工切坡等因素影响，在算例中施加水平方向20KN数值方向1300KN的体力作用。
取步长为0.00006s，在上述因素的影响下分别考虑300步、600步、900步、1200步及1500步时滑动变形过程，如图3.3所示。
图3.3边坡滑动变形过程
根据变形模拟可见，1200时步是边坡已经破坏，而且每个块体本身也有不同程度的变形。数值流形方法对模拟连续―非连续问题展现出了明显的优越性。
针对类似边坡的连续―非连续变形问题，分析表明，块体界面特性模拟与实际岩体还有一定差别，由于所遇到的大量情况是从连续性逐渐过渡到非连续性，而流形方法在岩体结构处采用不同的物理覆盖来描述岩移的不连续性，从而使得能更有效的模拟连续与不连续性的耦合问题。
4 总结
工程岩体十分复杂，它既不是连续体，又不是完全的离散体，既受工程、环境等影响引起小变形、大变形，更具有沿各类不连续面(结构面)相对变形、滑动、运动的大位移(含变形位移和刚移)，因此能够解决实际问题的数值流形方法无疑是一种较为适宜的方法，它综合考虑了工程岩体的复杂特性。
参考文献
[1]石根华著, 裴觉民译. 数值流形方法与非连续变形分析[M]. 北京：清华大学出版社, 1997
[2]赵光明编著.无单元法理论与应用[M].合肥: 中国科学技术大学出版社,2009 .
Lighthill利用保角变换的方法首先提出了二维翼型的反设计方法，Hicks，Murman和Henne等人将此方法发展为可应用机设计的工程设计方法。后Campbell等提出过一种带约束的直接迭代的表面曲率（CDISC）方法，Yu将其与N-S解算器耦合形成了一种翼型和机翼的设计方法。波音公司则将此方法发展成工程应用的设计方法，并广泛地应用于波音的B757，B777和B737NG等型号的设计过程，取得了很好的效果。例如在B777研制中由于使用了反设计方法，仅经过三轮机翼的设计便取得了满意的结果，使风洞实验的机翼模型大大少于过去B757和B767设计时的数目，充分表明了该设计工具的作用。可以说，反设计方法曾对民机设计起过革新性的推动作用；但反设计方法也有其固有的弱点（参见文献[13]的附录D）：首先，对于高度三维的流动要找到“好”的压强分布很困难；其次，不能保证所得结果为最优，即既具有高速巡航低阻的特性又在非设计条件下具有可接受的性能；最后，其他学科的约束会导致反复迭代。
低可信度CFD模型的数值优化方法
随着计算能力和数值优化方法的快速发展，应用基于CFD的数值优化方法于民机设计得到了很大的发展。这一方法的应用也从低可信度CFD模型开始，逐渐发展到采用先进的N-S方程解算器。波音公司发展了一种耦合TRANAIR[16]（一种全速势方程的有限元方法，可参见文献[13]附录B）和梯度优化方法的数值优化气动力设计方法，并在1992年形成了TRANAIR优化器的雏形[17]。经过近十年的改进，得到了一个适用于位势流/边界层耦合飞行条件的气动力优化设计工具[18-20]，具有多点优化设计能力，可处理高达600个几何自由度和45000个非线性不等式的约束条件（图1表示了TRANAIR优化过程示意图）。作为一个例子，图2给出了采用该软件对机翼/发动机短舱设计计算前后压强分布的对比，图a和图b分别表示了设计前后等马赫数线的分布。可以看出图a中挂架处出现激波；图b中短舱附近的机翼表面上消除了由于短舱干扰形成的激波。算例结果表明该设计软件可以处理很复杂的飞机/发动机综合设计问题。
高可信度CFD模型的数值优化方法现代优化算法可以分为依赖和不依赖梯度的方法两大类。
1.依赖梯度的优化算法
目前可用的大多数依赖梯度的数值优化方法都是从控制理论出发的，Jameson是此类方法的先驱者之一。尽管最初是由Pironneau提出利用控制理论进行椭圆方程系主控的外形优化的[21-22]，但Jameson首先提出了通过控制理论自动进行外形优化的伴随方程方法[23]并应用于跨声速流动。后来，Jameson和他的合作者，还有其他研究者，大力发展此方法，从全位势方程到Euler/N-S方程，从无粘设计到有粘设计，甚至从气动设计到气动/结构的耦合设计，形成了大量文献[24-36]。此方法不同于一般梯度优化方法之处在于它将外形作为一个自由表面，促使流动解和最终优化的外形同时趋于收敛，因而使优化方法具有很高的效率（其基本思想可参见文献[13]附录D）。
2.不依赖梯度的优化算法
最早无需梯度的优化算法有Powell（共轭方向法）[37]和Nolder-Mead的单纯形法[38]。最近Sturdza还应用后者于空气动力的设计[39]。近二十多年来人们更多地使用诸如模拟退火法[40]和遗传算法(GeneticAlgorithm-GA)等的搜索方法，特别后者更为人们所关注。Holland利用进化理论创造了遗传算法[41]（可参阅文献[13]附录D），即模仿生物的自然选择进行搜索以寻求最优解。与传统的搜索和优化方法相比，遗传算法具有下述4个特点[42-45]：1）不是直接作用于参变量集本身，而是对参变量集的某种编码运算。2）不是对单个点而是对多个点构成的群体进行搜索。3）直接计算适应值（函数），无需导数和其他辅助信息。4）利用概率转移原则，而非传统优化方法中的确定性原则。已有愈来愈多的研究和民机研制机构表现出了对这种随机寻优方法的浓厚兴趣，也已出现了不少利用遗传算法进行翼型或机翼优化计算的文献[46-56]。
3.对高可信度CFD模型数值优化方法的要求
分析最近十余年中出现的大量基于Euler/N-S方程的数值优化方法和文献，可以看出多数仍表现为学院式的探讨，提供可直接用于工程设计的方法和工具显得尚很有限，尽管已开始向这方面努力。这可能是因为：1）只是近几年来随DPW研讨会等的进行，数值模拟才可以比过去更正确地估算阻力值。2）工程界的空气动力外形优化需要在高维搜索空间中进行并存在大量的非线性约束，使优化问题十分复杂且计算开销巨大；3）巨大的计算量要求很丰富的计算资源和很长的计算时间，这与工程问题要求的迅速反馈相悖。
因此要使基于CFD的空气动力优化方法和软件成为日常的工程设计手段和工具需解决如下技术关键：1）具有建立准确计算诸如升力、阻力、力矩等敏感气动特性的正确流动模型的能力。比较现有的气动力优化方法可知，大多数方法还在使用不完善的流动模型，如基于Euler方程，甚至全位势方程等。虽然它们在一定条件下，如巡航小迎角飞行状态，可以提供合理的结果，但工程应用常要求准确地估算出阻力、俯仰力矩等敏感的气动特性，要求可计算整个飞行包线的飞行状态以及不同的复杂的几何外形等，这只能通过求解N-S方程来实现。顺便指出，有些文献（如文献[28]）虽以N-S方程为主控方程，但优化时的伴随运算子却是在没有考虑粘性流动的假设下得出的（参见文献[28]第6节）。为了提高计算准确度,最好在离散N-S方程时使用高阶的差分算子[53-54]。2）具有寻求全局最优的能力。通常基于梯度的算法容易陷入局部最优，而遗传算法等随机搜索的方法则具有取得总体最优的优点。3）能有效地处理大量几何和气动力的非线性约束。优化问题的最优解常常是位于不同维超曲面（hyper-surface）的交汇处，遗传算法不同于基于梯度的方法，不限于目标函数的光滑扩展，可应用于多重约束的情况[53-54]。4）可应用于不同的几何外形和设计条件。5）扫描高维搜索空间的计算有效性高，以满足设计周期和研制成本的要求。遗憾的是这正是遗传算法的主要缺点，即估算适应函数的高代价。可以采用多处理器上的有效并行计算来大大减少计算时间[57]，或在估算适应函数值时采用近似模型，如降阶模型[54,58]或响应面模型[50]等。
数值优化方法的发展现状和验证研究#p#分页标题#e#
1.空气动力优化设计计算的系列研讨会
近年来CFD学术界和航空业界都十分关注计算阻力的精度问题，这也是CFD应用于工程设计时所面临的第一个具有挑战性的计算。AIAA的应用空气动力学专业委员会在各方支持下，自2001年开始举行了DPW（DragPredictionWorkshop）系列会议[59]，参与者都用N-S方程求解相同的几何外形（翼/身组合体，翼/身/短舱/挂架的复杂组合体等），得到了一个巨大的计算结果数据集，可与已有的已经过修正的风洞试验值比较。由于参与的计算者所采用的数值方法、湍流模型、计算网格形式及数目等各不相同，此数据集可用作分析和讨论各种因素对CFD计算结果的影响。该系列会议至今已举行了5届，对推动和提高CFD计算阻力的精度很有意义。文献[13]的附录C中给出了前3届结果的分析和讨论。鉴于DPW系列会议的成功，AIAA应用空气动力学专业委员会针对CFD面临的第二个挑战---计算三维高升力外形的最大升力CLmax，于2009年发起并组织了类似的高升力计算研讨会，其第一次会议(HiLiftPW-I)已于2010年6月在美国举行,文献[60]是该次会议的总结。在上述工作的基础上，2013年1月AIAA又进一步在其ASM会议过程中形成了以加拿大McHill大学Nadarajah教授为首的空气动力优化设计讨论组，作为空气动力优化设计计算系列研讨会实际的组委会。讨论组讨论了：1）建立可供在一个有约束的设计空间中测试气动优化方法的一组标准算例。2）举行研讨会的时间。与会者一致认为，由于工业界对基于CFD的气动外形数值优化方法有强烈的需求，优化方法和工具的研制也已有了相当的发展，可以以类似于DPW的研讨会形式，通过对一系列复杂气动外形的优化，来评估现有的寻求最小阻力外形的各种优化方法的能力，并将结果向工业界/研究机构公布。与会者还认为第一次会议从二维和三维机翼外形开始是合适的，并请加拿大的与会者准备标准算例。第一次会议拟于2013或2014年的AIAA应用空气动力会议期间举行。
2.先导性的研究
事实上为准备此研讨会，波音的Vassberg，斯坦福的Jameson，以色列的Epstein及Peigin等三方从2007年起即开始了先导性的研究（pilotproject）,以积累经验和发现问题。三方用各自己开发的优化软件（MDOPT，SYN107，OPTIMAS）对第三届DPW会议的测试机翼DPW-W1独立地作优化计算[61,62]。波音研制的MDOPT[63]（也可参见文献[13]的1.7.3节）可使用响应面模型（InterpolatedRe-sponseSurface—IRS）的数值优化格式[64]，也可直接从流场解计算设计变量的灵敏度代替IRS模型完成优化。其流场解软件为TLNS3D[65]，计算网格点为3582225。Jameson开发的SYN107采用基于梯度的“连续”伴随方程方法[23,31]，其流场解软件为FLO107，计算网格点为818,547。
以色列航空公司开发的OPTIMAS采用降阶模型的GA算法，流场解软件为NES[66-68]，计算网格点为250,000。对三方独立优化后所得的外形再用不参与优化的流场解软件OVERFLOW[69]作评估计算，计算网格点数为4,000,000，以便能准确地计算阻力。结果表明，4个分析软件计算得到的阻力增量值的分散度在Ma=0.76时为5counts(1count=0.0001)，Ma=0.78时为10counts，因此很难确定哪个优化后外形最优。但从Ma=0.76，C=0.5单设计点的阻力改进结果（表1）[61]看，OPTIMAS优化后的04外形明显优于MDOPT优化后的M5和SYN107优化后的S4。文献[61]还讨论了从比较中可吸取的经验和教训。
一种基于高可信度CFD模型的数值
优化方法的构造本节将以OPTIMAS为例对如何满足可应用于工程实践的高可信度CFD模型数值优化方法的要求做一说明。
1.优化方法的构造及其特点
OPTIMAS是将遗传优化算法和求解全N-S方程的分析算法相结合的一种有效并鲁棒的三维机翼优化方法。1)其全N-S方程的流场并行解算器NES[66-67]基于高阶低耗散的ENO概念（适用于在多区点对点对接网格中的多重网格计算）[66,71]和通量插值技术相结合的数值格式，采用SA湍流模型，可快速准确地完成气动力计算，因此具有计算大量不同流动和几何条件的鲁棒性。作为例子图3和4给出了ARA翼身组合体Ma=0.80，Re=13.110时的升阻极线和CL=0.40时的阻力发散曲线[68]，使用的网格点数分别为，细网格（3lev）：900,000，中等网格（2lev）：115,000。可见升阻极线直到大升力状态的计算与实验都很一致。对比图中还给出的TLNS3D在细网格（2,000,000）中的计算值可见，无论升阻极线或阻力发散曲线NES的都更优。作为数值优化软件的特点之一是其在流场解算器中首次使用了高精度格式。2)优化计算的遗传算法中采用了十进制编码、联赛选择算子[42]、算术交叉算子、非均匀实数编码变异算子[72]和最佳保留机制。为解决搜索时总体寻优耗时大和求解N-S方程估算适应函数代价高的问题，在寻优过程中估算适应函数时采用当地数据库中的降阶模型[54,58]获取流场解（当地数据库是在搜索空间中离散的基本点处求解全N-S方程建立的），并以多区预测-修正方法来弥补这种近似带来的误差。多区预测-修正方法即在搜索空间的多个区域并行搜索得到各区的优化点，再通过求解全N-S方程的验证取得最优点。为保证优化的收敛，寻优过程采用了迭代方法。3)在整个空间构筑寻优路径（图5），扩大了搜索空间和估算适应函数的区间[54]。4）为提高计算效率，OPTIMAS包含了五重并行计算：Level1并行地求解N-S方程Level2并行地扫描多个几何区域，提供多个外形的适应函数的计算（level1隐于level2中）。Level3并行的GA优化过程（level3隐于level4中）。Level4并行地GA搜索多个空间。Level5并行地生成网格。5）采用单参数或双参数的BezierSpline函数对搜索空间参数化;并基于优化外形与原始外形的拓扑相似自动地实现空间网格的快速变换。
2.优化设计的典型结果
文献[53]~文献[58]给出的大量算例充分表明了OPTIMAS优化软件的优异性能。本文5.2中给出了其优化三维机翼的性能，这里再补充两例。1）翼身组合体整流（fairing）外形的优化文献[73]讨论了某公务机翼身组合体机翼外形优化的单点和多点设计两者性能的比较。结果表明，多点优化设计能同时保证设计的巡航状态时，和高Ma数飞行，起飞等非设计状态时的良好性能。文献[74]进一步讨论了翼身组合体整流外形的优化设计。流动的复杂性（三维粘流/无粘流强相互作用的流动区域）和几何的复杂性（三维非线性表面）使整流外形的设计经历了传统的试凑法，基于Euler解的试凑法等，最后才发展为现代完全N-S解的数值优化方法。文献[74]采用了这种方法，先作机翼外形优化，再作整流外形优化，然后再作机翼优化，整流外形优化，……依次迭代，直至收敛。优化中用双参数的BezierSpline函数将整流外形参数化，所得搜索空间的维数ND=（2N-2）*（M-1）决定的参数化整流外形与实际外形的差别在M=10，N=4，ND=54时可准确到0.3mm（满足工程需求）。计算网格数为90万。表3给出了设计条件和约束，表4给出了设计点的阻力值比较。由表4可知，GBJ2的减阻为16.7，50%DC，GBJFR1的减阻为10.7，32.1%DC，GBJFR2的减阻为5.9，两次优化机翼的减阻总计为22.6，67.9%DC，优化机翼和优化整流外形减阻作用分别约占2/3和1/3，可见整流外形的优化也是十分重要的。约束则使减阻损失4.6（如GBJFR3-GBJFR1）。图6至图9分别为原始外形，GBJ2，GBJFR2和GBJFR4的整流处等压线分布，可见整流外形的优化消除了原始外形和GBJ2中存在的激波。图10和图11分别给出了Ma=0.8时升阻极线和CL=0.4时阻力发散曲线的比较，可见优化设计不仅对设计点，对非设计状态也都有好处。2）翼身融合体飞机气动外形的优化[75]优化设计以英国克朗菲尔德大学设计的BWB外形[76]为出发外形，该外形的主要设计点为，。在数值优化计算中还考虑了，的第二个设计点和，（起飞状态）的第三个设计点。几何约束有剖面相对厚度，前缘半径，后缘角，每个剖面的樑处还附加两个厚度约束，其中上标b表示出发外形，*表示优化外形，下标i表示第i个剖面。附加的空气动力约束为对俯仰力矩的规定。采用Bezier样条描述几何外形，总设计变量为93个。表5给出了设计计算各状态的条件和约束，其中是权系数。表6给出了优化计算结果。#p#分页标题#e#
单点优化的BWB-1结果与文献[77]的结果相比较可见，文献[77]采用Euler方程的无黏优化使阻力降低了26counts；而这里的BWB-1全N-S方程优化使阻力降低了52counts，显示了此黏性优化方法的优点。比较有、无俯仰力矩约束时优化得到的BWB-2和BWB-1表明，尽管BWB-1阻力降低的效果突出，但其值过大，出于稳定性考虑而不能接受；BWB-2的阻力虽比BWB-1大了1.9counts，却满足了力矩的要求。表6中的双点优化设计（BWB-4），使第三设计点（低速状态）的达到1.671（消除了BWB-2达不到设计要求1.63的缺点），且基本保持了主设计点的阻力收益，为196.6。然而BWB-4在时的阻力达216.6，高于BWB-2的213.4，表明需要三个设计点的优化设计（BWB-3）。BWB-3在时，为202.5（比两点设计值减小了14.1），同时满足了其它两个设计点的性能要求。图12至图15给出了所有设计状态和时的极曲线，时的阻力发散曲线和时的随迎角α变化的曲线。由图可见，时所有优化设计的极曲线都非常接近，相比于原始外形的极曲线，性能有了很大改进；时也一样，特别是三点优化设计的BWB-3，优点更明显。阻力发散曲线也都有了很大改进，在前所有的总阻力基本保持常值，单点与两点优化的阻力发散点接近，而三点优化的可达附近。由图15可知，没有考虑低速目标的BWB-1和BWB-2具有较低的，将低速目标计入设计状态的BWB-3和BWB-4所得的皆优于原始外形的。上述结果表明三点优化设计具有最佳的优化效果和总体最好的气动性能。Fig.15LiftcoefficientCLvsangleofattackatMa=0.2最后，上述各优化结果在（主设计点）时的阻力值基本相同，但几何外形却差别不小，由此可见，外形阻力优化问题没有唯一解[75]。上述计算是在具有456GBRAM，114MB二级高速缓存的机群环境下通过“过夜”方式完成单点优化设计，在1.5-2天的计算时间内完成三点优化设计的，计算时间可满足应用于工程设计的要求[75]。
结束语
随着电子技术和通信技术的发展，无线通信以及遥测遥控系统被广泛应用于工业、农业、航空、航海等各个领域中。出海口及内陆河道作为航海航运重要的一部分，其管理维护方法及管理质量对我国航运业的影响至关重要。发展至今，电子通信产品的可靠性越来越高，成本越来越低，这使得航道管理维护自动化、数字化的实现成为可能。GPS(全球定位系统)是美国国防部于1973年开始研制的卫星全球导航定位系统，主要为其海陆空三军服务。近几年来已逐步应用于民用设施及测绘技术中，同时美国军方逐步放松对民用GPS设备的限制，使得民用GPS达到了比较高的定位精度。利用GPS对航道航标等设备进行位置遥测与监控是一种比较理想的方法。本文以航标监控的具体要求为标准，把整个航道管理区域内需监控的目标物组成一个GPS遥测网，并利用各种滤波方法消除相应的误差，提高了遥测数据的准确性。
1 GPS OEM板与航道GPS遥测网
1.1 GPS OEM板
GPS OEM板是GPS接收机中一个重要的组成部分，它具有成本低、体积小、重量轻、产品种类多、性价比高等很多优点，因此被广泛应用于定位及导航领域中。它的定位精度已经能达到几十米，甚至可以达到10米以内的精度。本课题所用到的Thales集团导航定位公司的GPS OEM B12就是一款性价比很高的产品。
1.2 航道监测
航道是交通网络中一个重要组成部分，其安全质量直接影响着整个交通系统。以前航道部门专门在航道的堤岸、桥头、故障物旁边安装各种航标灯作为警戒导航装置，各种船只可以根据航标灯光及其闪动频率来确定自己的航向。至于航标的维护，则是航道部门每隔一定时间派巡航船只对各航标灯进行目测和实测。因为航道中航标灯比较多，这就使得这种巡航航道的维护方式操作繁琐，运作维护成本高，安全质量低。
1.3 航道GPS遥测网
航道中航标遥测网主要是对水标(抛锚在水中的航标)进行遥测以便对其位置进行实时监控(其系统原理图如图1所示)；而岸标(固定在堤岸上的航标)由于其位置不变所以无需GPS遥测。GPS在航标遥测网中的实际任务就是实时测量航标灯所在位置，并与预先划定的位置范围进行比较，如果漂离出所标定的范围，即通过GSM网发送警报信息给监控中心，以便于监控中心采取相应措施。这将就可以排除航标灯因船只碰撞、水流冲击等原因而漂离引起事故。而每个航道管理区域内有成百个水标，因此在提高安全质量的同时也需考虑成本投入。根据航道的具体要求，其精度并不需要精确到米级以下，因此不需要价格昂贵的高精度GPS接收机及测量仪。同时将GPS OEM板与水标进行捆绑，可以以相对较低的成本取得高质量的管理效果。本系统使用的是法国Thales公司生产的B12 GPS OEM板模块，它具有并行的12个接收通道(即同时可以接收12颗定位卫星传送的星历信息)。
2 误差分析、数值处理及控制流程
2.1 误差分析
GPS测量的误差主要包括卫星部分、信号传播、信号接收等各个方面带来的误差，但从性质上来讲可以归纳为系统误差和随机误差两部分。其中系统误差主要包括卫星的星历误差、卫星钟差、接收机钟差以及大气折射的误差等。随机误差主要包括信号的多路径效应等。虽然系统误差比随机误差要大些，其消除主要靠接收机本身[1]，但是它总是有一定的规律可循的，所以采取一定的措施进行处理对整个系统的可靠性都是非常重要的。由于水面多路径效应比较严重，所以使用精密相位中心、具厄流圈的测量天线是消除由于水面环境所引起误差的一个重要方法。
2.2 数值处理
针对各种误差，测量技术中已应用了各种滤波方法来消除或减弱各种误差的影响，例如中值滤波法、算术平均滤波法、进退递推滤波法等。通过大量的测量试验与观察分析发现，随着时间的不同、卫星分布状态的改变以及天气的变化，GPS所读数据都有不同曲线方向的飘移，但是其分布状态接近于正态分布，所以采用一些滤波方法对数据进行处理对整个测量系统精度的提高至关重要。以下是系统中所用到的几种滤波方法。
中值滤波法：即对所测三个数据进行排序，去掉最大和最小的一个，取中间值作为测量值。基于这种思想，本文在终端控制器上电初始化的时候连续测量n(可调)次经纬度数据并将它们从小到大进行排队，去掉最大的m次数据和最小的m次数据，以中间的n-2m次数据作为基准，并存于一个存储单元。由于航道遥测系统对实时性要求并不高，所以把n尽量取得大些。设n次所读数据和为Xn，经排序后最小m次数据和为XmMIN，最大m次数据和为XmMAX，则：
Xsum=Xn-XmMIN-XmMAX
把Xsum存于存储单元作为后续处理方法的和基准。 算术平均滤波法：即采样一定量的数据，然后对其求平均值作为测量估计值，这样可以使得偏离真值的正负误差相消，从而使测量值更接近真实值。本课题将前面所取得的n-2m次测量数据作算术平均，且存于固定的算术平均值存储单元，并根据以后所读数据进行实时修正。这样有：
X=(Xswn)/(n-2m);Xi=(Xsumi)/(n-2m).
其中，X是初始化时所求平均值，作为一个平均基准存于存储单元。Xi是每读一次数据所求平均值，作为位置评估值应用于位置飘移判断控制中。
进退递推滤波法：前面两者都是读取一定数据以后再作后处理，而测量过程中必须对所测数据进行实时处理。所以，所测量经纬度的变化趋势必须反应出来，以便航标因为意外而漂出所给定范围时能实时向监控中心发送警报信息，从而进行修正。本文根据实验与观察的结果，采取进一新数退一平均数的进退递推滤波方法，即：
Xswni=Xsum_i-1+Xi-1+xi
限幅滤波法：在测量过程中，常常会碰到偏离中值较远的粗大误差。这对经过前面几种滤波法处理后的数据基准会产生较大的冲击，限幅滤波法就是针对这一思想的。设定一个阈值，当所测数据与基准数据比较后，差值超过阈值就认为是粗大误差并舍掉。但是本课题中如果航标灯因意外而漂出很远，就必须能识别出来，而不能当粗大误差全部舍掉。所以在控制程序中专门设计了一计数器对舍掉比率进行计数，如果舍掉比率大于某一值则重新初始化，即重新读取n-2m次的和基准及其算术平均基准。
图2、图3、图4分别是对利用Visual Basic6.0开发的数据采集与处理程序采集的10小时GPS数据进行几种数据处理后的坐标示意图(其中，横坐标、纵坐标分别表示经、纬度)。从这三个图中可以看出，从图2到图4，数据收敛性依次增强，可见综合几种滤波法于数据处理中，将大大减少误差，提高系统精度。
2.3 控制流程
摘要：数值仿真是一类基于计算机求解，针对实际工程问题采用仿真模拟，从而进行数值实验的方法。在工科类课程中引入数值仿真方法，不但有助于学生深入理解复杂工程现象，提升实践能力，还有助于学生基于数值仿真的现代工程分析方法，提升学习兴趣与创新思维能力、理论学习水平与工程实践能力，培养早期科研思维。
关键词：数值仿真；工科课程；应用研究
数值仿真是一类基于计算机求解，针对实际工程问题采用仿真模拟，从而进行数值实验的方法，已成为与实验研究、理论研究相并列的研究科学问题的必备方法。基于数值仿真，可通过图像显示所研究问题发生的物理、化学过程，如天气预报中，即可通过计算机仿真的方式预测某地区的温度变化过程，具有直观性强、成本低、可研究问题范围广的优点。如今，数值仿真方法已被应用于力学、机械工程、能源动力工程等工科类专业课程教学中，为复杂或不易开展的实物实验提供数值实验途径，加深学生对课程中复杂的物理、化学原理的理解。
一、数值仿真应用于工科类课程理论教学的具体方法与特点
数值仿真是利用数值计算原理与计算机求解技术开发出的现代科研方法。基于所求解问题的物理、化学原理，采用有限元等数值分析方法，将连续的物理、化学问题离散求解，获取所求解对象的物理、化学量分布特性。目前，数值仿真方法广泛应用于科学研究，如高温燃烧热场的温度分布、高超音速流动的流场压力分布、天气温度预测等。在工科类课程教学中，由于理论、概念多，物理、化学过程复杂，常制约着学生的理解，而数值仿真方法可通过图像显示的方式展示所研究问题发生的物理、化学过程，具有可操作性强、结果显示直观的特点，因此，可将数值仿真方法引入工科类课程教学，提升学生对课程知识的理解。一般的数值仿真过程分为预处理、求解、后处理三部分，学生在掌握基本数理知识的基础上，即可理解数值仿真的基本流程。在预处理中，需建立所分析问题的几何模型，并给定求解条件。例如，分析钢结构塔在承受一定重量时是否具有足够的强度，需首先基于分析软件建立钢结构塔的外形几何模型。在求解环节，基于数值计算方法，利用分析软件中求解仿真模型以获取需掌握的物理、化学量。在后处理环节，完成分析结果的可视化，以图像形式直观展示分析结果，如天气预报中展示的温度云图，即是后处理结果。在工科类课程中，涉及大量有关物理、化学现象的理解，而某些环节无法通过直观实验复现，也无法通过理论分析获得结果，对学生的理解带来困难。例如，在锅炉设计中通过改进折焰角的尺寸，即可优化锅炉炉膛内部的燃烧热场，而炉膛内部燃烧的温度无法通过理论计算获取，也无法直接进行实验测试，但可基于数值仿真获取其燃烧温度分布图，直观展示折焰角尺寸调整对热场分布的影响，有利于学生通过实践操作，获取分析结果，加深理解。
二、数值仿真应用于理论教学的具体实例———材料力学弯曲应力求解教学
材料力学是机械、土木、能源等工科类专业的必修专业基础课程，是结构设计的必要基础，基于材料力学基本原理，可分析部件的强度、刚度及稳定性，为工程结构设计提供指导。材料力学中的杆件弯曲问题是重要的教学知识点，当矩形截面杆件受弯时，杆件横截面上的正应力呈线性分布，在横截面上、下边缘正应力最大，在中部中性轴处最小。在理论求解中，需首先根据给定约束条件与外部受力特点，分析杆件受到的弯矩，然后利用理论公式求解正应力。弯曲应力求解为材料力学中的重点与难点内容，区别于杆件的轴向拉压与扭转问题，弯曲问题的应力分布更复杂，学生在学习时不易获取直观理解，故可借助数值仿真方法。利用受力求解的数值仿真软件，在预处理环节，首先建立杆件的几何模型，根据求解要求设置约束，其后在给定位置施加受力，最后，指定杆件所使用的材料并由数值仿真软件自动完成求解。该问题为静力学求解问题，求解精度较高。在完成求解后，进入后处理环节，在该环节可通过软件直接显示梁在弯曲作用下的变形，也可以通过图像显示梁的应力分布。藉此，直观清晰地向学生描述了相关力学概念，展示了复杂的受力情况，使学生对课程理论学习产生更直观的认识，既提升了学习兴趣，又进一步加深了对课程内容的理解。
三、数值仿真应用于理论教学的结论
利用数值仿真方法，可以直观、清晰、准确地显示复杂物理、化学现象的发生过程，而工科类课程中涉及大量物理、化学概念的讲解，引入数值仿真教学环节，不但有助于学生深入理解相关概念，提升对课程的掌握程度，还有助于学生基于数值仿真的现代工程分析方法，提升学习兴趣与创新思维能力、理论学习水平与工程实践能力，培养早期科研思维，助力学生今后发展。
作者：李一飞
关键词:汽车覆盖件;数值模拟;有限元;冲压成形
中图分类号:TG386文献标识码:A
文章编号:1009-2374 (2010)21-0042-02
汽车覆盖件一般由钢板冲压而成,冲压成形是一种非常复杂的力学过程,用传统方法很难求解。近年来,随着计算机软硬件技术、图形学技术、人工智能技术、板料塑性变形理论和数值计算方法等的发展,以及与传统的工艺/模具设计技术的交叉集成开创了利用CAD/CAM/CAPP技术和CAE数值模拟分析技术进行覆盖件成型工艺设计的新领域。板料冲压过程的计算机分析与仿真技术已能在工程实际中帮助解决传统方法难以解决的模具设计和冲压工艺设计难题,如计算金属的流动、应力应变、板厚、模具受力、残余应力等,预测可能的缺陷及失效形式,如起皱、破裂、回弹等。在汽车覆盖件的设计中采用数值模拟技术能从设计阶段准确地预测各种工艺参数对成形过程的影响,进而优化工艺参数和模具结构,缩短模具的设计制造周期,降低产品生产成本,提高模具和冲压件产品质量。
1冲压成形数值模拟理论
板料成形过程及特点决定了其成形是涉及几何非线性、材料非线性和边界条件非线性的弹塑性大变形力学问题,如果用传统的理论分析方法来分析成形过程是不可能的,甚至根本无法实现。长期以来,国内外学者对板料成形性能、成形过程中应力、应变分布的研究基本建立在实验或经验公式的基础上。随着有限元数值模拟理论技术的发展,人们开始把眼光转移到其在汽车覆盖件冲压成形的应用上来。经过多年的研究,板料成形有限元技术在材料本构关系、单元技术、接触算法、求解格式等方面得到了发展。
1.1本构关系
目前在汽车覆盖件冲压过程进行分析中,凸模、凹模及压边圈在冲压过程中的变形小,通常采用刚体材料模型,而对于板料大多采用弹塑性本构关系,对于不同的金属有不同的弹塑性模型可以选择。
建立弹塑性本构关系模型首先要解决复杂受力情况下屈服状态以及屈服后的塑性流动,解决复杂受力情况下屈服状态就要建立屈服准则。冲压成形领域中经常采用的屈服准则有:von Mises屈服准则、Hill屈服准则以及3参数Barlat屈服准则。在早期的冲压分析中,板料被假设为各向同性材料,因此经常采用von Mises屈服准则,后来随着有限元的发展,研究人员证明板料是各向异性的。Hill提出了用二次函数来描述正交各向异性材料的塑,即Hill屈服准则。但Hill屈服准则却无法正确分析分析多晶体塑性材料,因此人们进一步研究建立了许多屈服函数和屈服准则来描述多晶体塑性材料。例如Barlat等人提出了一种形式化的方法来描述多晶体材料的屈服准则。
经实验分析表明:当厚向异性系数r较小时,使用Hill屈服准则建立的材料模型,计算结果误差很大,甚至大于使用von Mises屈服准则的材料模型。而采用3参数Barlat屈服准则进行分析时,则能够得到满意的结果。当厚向异性系数r较大时,则Hill准则和3参数Barlat屈服准则都能获得正确的结果。3参数Barlat屈服准则的结果要优于Hill准则,vonMises屈服准则结果最差。因此在汽车覆盖件冲压成形分析中3参数Barlat屈服准则是最常用的材料模型。
1.2单元技术
用于冲压成形有限元分析的单元有三种:基于薄膜理论的薄膜单元、基于板壳理论的壳单元和基于连续介质理论的实体单元。薄膜单元格式简单,但忽略了弯曲效应,因而只适用于分析胀形等弯曲效应不明显的成形过程。在薄板壳的成形分析中,又因为薄膜理论是二维理论,因此薄膜单元只适合二维成形问题分析。实体单元虽然考虑了弯曲效应和剪切效应,但由于计算时间太长,除非板料厚度非常大的情况下,一般在汽车覆盖件成形分析中不采用实体单元。基于板壳理论的壳单元不仅考虑了弯曲效应和剪切效应,而且板壳单元是处理薄板三维变形的工具。因此,在汽车覆盖件成形分析中常采用壳单元。
对于薄壳单元,人们提出采用Kirchhoff理论和Mindlin理论其应力或应变状态进行简化。Kirchhoff理论需要构造C1连续性插值函数,在三维分析中构造C1连续性插值函数是非常困难的,构造的壳单元效率也很低,因此在冲压成形分析中不采用基于Kirchhoff理论的C1型壳单元。Mindlin理论采用位移和转动独立插值的方法,从而使问题简化。近年来人们开发了很多种基于Mindlin理论的壳单元,例如BT壳单元,由于其计算结果准确、计算效率高,因此常用来建立汽车覆盖件成形分析中板料的有限元模型。
1.3接触算法
板料变形时,接触发生的时间和位置随着接触体的变形而改变,用有限元处理接触问题时必须建立正确的接触问题模型。接触界面的处理实际上找出所有接触对及状态,然后计算每个接触对的作用力。前者需要解决的是接触点、接触区域的搜索及接触状态,后者需要解决的是接触区域间法向接触力和切向摩擦力的计算。在进行有限元分析时寻找接触对的方法通常采用增量搜寻或桶式分类搜寻。接触力的计算主要应用的是罚函数法,切向摩擦力的计算采用修正的库仑摩擦定律。
2数值模拟软件
经过多年的发展,利用冲压成形模拟技术和相关理论,人们已经可以对部分板材冲压加工过程进行准确模拟,并且人们开发了许多商业软件应用于生产实践中,通常软件的开发往往基于不同的原理,不同的软件反映了冲压成形分析中有限元方法的差异,例如按变形原理可以分为基于刚塑性变形的SHEET-3软件和基于弹塑性变形的Auto-Form、PAM-Stamp和Dyna-Form软件,按求解格式又可以分为基于静力隐式格式的Auto-Form软件和基于动力显式格式的PAM-Stamp和Dyna-Form软件。虽然基于不同的原理,但实践表明利用这些软件对板料成形过程进行模拟从而指导实际生产过程的方法是切实可行的。
但是由于汽车覆盖件本身的复杂性,覆盖件冲压成形的影响因素极其复杂,覆盖件冲压成形涉及的领域极广,所以对汽车覆盖件冲压成形问题的研究依然存在许多问题,例如仿真建模的合理性和准确性;材料屈服模型;计算效率和计算精度问题;回弹问题等。这些问题涉及复杂覆盖件成形模拟的关键部分,因此它的解决定会使汽车覆盖件成形的数值模拟产生质的飞跃,因此也成为人们关注的重点。
3结论
随着计算机技术和数值计算方法的发展,有限元数值模拟技术在汽车覆盖件成形工业中发挥着越来越重要的作用。利用它可以指导实际的冲压成形过程,可以实现新产品开发周期短、质量高、低成本的目标。目前板料数值成形技术在汽车覆盖件制造领域的应用越来越广泛,经比较和分析表明采用3参数Barlat屈服准则,单元类型为BT壳单元和求解格式为动力显式格式的有限元方法更适于汽车覆盖件冲压问题的分析。
参考文献
[1] 危熠平,王健,雷君相.汽车覆盖件冲压模具仿真设计[J].模具工业,2005,(10).
[2] 代洪庆,刘晓晶,闫巍,刘江涛.汽车覆盖件冲压成型的计算机仿真[J].机械工程师,2006,(5).
[3] 林忠钦.车身覆盖件冲压成形仿真[M].北京:机械工业出版社,2005.
[4] 王勖成,邵敏.有限单元法基本原理和数值方法[M].北京:清华大学出版社,1997.
关键词： 空心铝型材 有限元法 有限体积法 数值模拟
1.引言
目前，我国的铝型材加工企业在挤压工艺制定和挤压模具设计过程中，很大程度上依赖设计人员的经验，导致模具设计周期长、成本高，且产品质量得不到保证。因此研究铝型材挤压过程的机理和金属流动规律，进而为工艺和模具设计提供理论指导，就显得尤为重要。通过研究发现目前铝型材中的实心型材多采用有限元法和有限体积法，而针对空心型材的数值模拟研究很少。因此研究空心型材挤压过程的数值模拟方法对于合理设计挤压工艺和模具结构、减少试模时间和提高型材质量都具有十分重要的意义。本文以矩形管铝材的挤压为例，分别采用有限元法和有限体积法进行数值模拟，分析两种方法对于空心型材挤压模拟的可行性和模拟精度。
2.数值建模与模拟结果对比分析
矩形管尺寸60mm×20mm，壁厚3mm（见图1），坯料选择Al6063合金，尺寸170mm×50mm，初始温度为480℃，模具材料H13，预热温度为450℃。
对于矩形管挤压过程的有限元模拟，采用DEFORM-3D模拟软件，有限体积模拟则采用MSC Superforge软件。利用pro/E建立模具三维几何模型如图2所示。
（a）上模实体模型 （b）下模实体模型
在保证挤压进入稳定状态的前提下，有限元和有限体积模拟中均设置凸模压下量为38mm；有限元模拟时初始步长增量设置0.3mm，网格尺寸1.5mm，局部细划网格尺寸1mm；有限体积模拟时划分网格尺寸为2mm。选取摩擦因子为0.33，凸模速度为2mm/s。
图3为矩形管挤压过程有限元模拟结果，模拟在凸模压下量28mm处停止。从模拟结果可以看出，在挤压材料流入分流孔的过程中，模拟过程进展顺利，挤压件表面质量较好（见图3（a））。当变形材料流入焊合室以后，模拟计算变得比较困难，需要频繁的网格重划，所需模拟时间显著增加，模拟结果精度下降，但模拟仍然可以有效进行（见图3（b））。当变形材料进入工作带时，模拟过程变得异常困难，几乎每经过一步模拟就必须重新划分网格，体积损失严重；在焊缝处，由于频繁的网格重划，相邻分流孔内流出的金属无法焊合在一起，甚至出现体积缺失的情况，仿真结果已经严重失真，难以继续进行下去（见图3（c））。
对比图3（b）和图3（c）可以看出，由于有限元方法对此矩形管材的挤压模拟因为长宽不等而无法避免金属网格节点的自焊合问题，且该模型使用的是四个分流孔，所以其1/4模型总是存在两股料流在焊合室接触重新焊合的问题。分析原因是：①矩形管挤压为大变形挤压问题，有限元模拟过程中，由于网格畸变严重，需要频繁的网格再划分，每次网格划分都存在几何上的近似，即每次划分的新网格与旧网格总是存在一定的差异，当几个模拟步骤内划分次数较多时，几何形状就会发生较大差异。从模拟出错时的网格图4可以看出，焊缝处的三角平面上的每一条边并不是由两个三角平面共享。这就意味着在实际模拟中，表面上变形材料充满了焊合室完全焊合，实际还是两股料流，在后面的步骤里可见断裂，导致计算失真停止。②每次网格划分后，边界节点与模具边界的接触条件也相应地进行了重新确定，频繁的网格划分引起的几何形状的变化会导致部分边界节点并没有发生真正的接触，而边界节点的接触条件是导致模拟结果（特别是变形体几何形状的模拟结果）失真的主要因素。
此外，图3（c）表示凸模行程为28mm时的模拟结果，从模拟开始到图3（c）表示的位置，总计模拟时间50h，网格重新划分了92次。因此，无论是模拟精度还是运算效率，有限元法都不是模拟矩形管铝材挤压过程的有效算法。
图5为利用有限体积法模拟矩形管分流模的挤压结果。从结果可以看出，由于避免了有限元法中频繁的网格重划，有限体积法可以很好地解决有限元模拟矩形管出现的变形材料无法自接触问题，模拟进展顺利。
3.结语
本文分别采用有限元法和有限体积法对Al6063矩形管的挤压成型过程进行了数值模拟，比较了两种方法的模拟结果，发现在对矩形管成形过程进行模拟时，由于有限元模拟中频繁的网格重划导致网格几何外形和边界节点接触条件的失真，使得在焊缝处出现变形材料无法焊合的问题，最终模拟失败；有限体积模拟中由于无需网格重划，可成功模拟整个热挤压过程，为空心铝型材挤压模具结构的设计和工艺参数的优化提供了依据。
参考文献：
［1］周飞，苏丹，彭颖红等.有限体积法模拟铝型材挤压成形过程［J］.中国有色金属学报，2003，13（1）：65-70.
［2］于沪平，彭颖红，阮雪榆.平面分流焊合模成形过程的数值模拟［J］.锻压技术，1999，（5）：9-11.
［3］王锐，赵国群，吴向红等.大口径铝管平面分流模挤压过程数值模拟［J］.模具工业，2007（33）：17-21.
［4］马思群，孙彦彬，李荭.金属挤压成型的数值模拟技术研究［J］.大连铁道学院学报，2002，23（2）：8-11.
以裂缝为主要渗流通道的裂缝性油藏在碳酸盐岩油藏、低渗透性油藏中都占据着相当大的比例，而压裂开发又是其重要的开发方式。如何对大裂缝（包括天然大裂缝及人工压裂裂缝）进行科学而有效的模拟是影响开发效果预测的重要问题。
1 大裂缝油藏数值模拟处理方法
1.1 网格化表征法
网格化表征法主要采用裂缝网格化技术来显式地描述大裂缝的性质（包括走向，形态，开度、长度等），主要有密网格法，局部网格加密法，非结构网格法等。本文中以局部网格加密法为代表。
局部网格加密法可以显式地对裂缝进行建模和描述，在描述及显示驱替过程方面表现较好。但有可能因为局部裂缝孔隙体积过小而引起收敛性困难，特别是在裂缝中的重力分异过程描述水线突进时。另外对于复杂裂缝系统，局部网格加密法对于裂缝的描述非常复杂。
1.2 等效渗流特征描述法
等效渗流特征描述法主要采用流动能力的等效计算对裂缝的渗流能力进行等价。
等效级差法是该类方法的代表。主要做法是在原有网格系统的基础上，根据传导率等效原则，修改包含裂缝网格的渗透率，近似地等价裂缝渗流效果。该方法由于使用方便，只需要修改网格的渗透率，因此仍然是目前经常使用的近似裂缝模拟方法之一。但该方法在计算的过程中经常会出现即使将裂缝处的渗透率改的很大，仍见水时间晚的情况。
1.3 双重介质描述方法
双重介质模型把发育的互相连通的裂缝看成是一种连续介质，同时把被裂缝切割的岩块也看作一种连续介质。两个连续介质在空间上是重叠的，即每个几何点既属于裂缝连续介质也属于基质。裂缝和岩块中的流体按照一定规律进行交换。
它能够既能体现裂缝系统高渗高速流动的特性，同时还能反映渗吸效应、重力效应、分子扩散效应（气驱）等驱替机理，是目前裂缝性油藏描述中较完善的模型，但它更适用于中小裂缝且均匀分布的情况，在大裂缝非均质性较强的情况下对裂缝位置难以准确描述，表现水突进的直观性差，且计算收敛性较差。
2 单重介质传导型裂缝模型
Paul V. L.[1]等人提出了一种提出的一种在单重介质模型中加入传导型裂缝的新的处理方法。这种方法是针对裂缝尺寸与网格块相比较大，使用传统的方法已经无法描述的情况提出的，主要适用于基质作为主要储油介质，且裂缝之间距离较大（数十米）的油藏。
该方法中裂缝的处理方式与断层类似，主要考虑裂缝的传导性。具体应用通过修改网格的属性、生成拟相对渗透率曲线函数、修改井的生产指数来实现。
使用E c l i p s e数值模拟软件建立8×8×1的正方形概念模型，网格大小为22m×22m×9m。设计一注一采两口井位于对角线两端，定注采比为1，裂缝沟通注水井和生产井。假定裂缝的开度为0.01m，设计裂缝的有效渗透率为50MD，500MD，5000MD，计算对比各方案的含水上升情况。
模型对裂缝渗透率因素敏感，裂缝有效渗透率越大，生产水见水越早，含水上升越快。
3 各种方法综合对比与适应性分析3.1 概念模型建立
为适应网格化表征法，表现出基质储油，裂缝导流的机理，建立以下的概念模型：x方向平行于裂缝，y方向垂直于裂缝，在y方向将裂缝划分为一个网格，网格宽度为裂缝缝宽，x方向划分为多个网格。为减少数值困难，y方向网格划分时，网格大小在裂缝附近选小些，离裂缝越远，网格尺寸变大[2]。设置一口注水井，一口生产井，均完井于裂缝上。裂缝缝宽取0.01m，有效渗透率为5000md。
运用网格加密法、等效级差法、双重介质模型及单介质传导模型分别进行计算。3.2 计算结果分析
图1 四种方法计算含水上升规律
四种方法计算的含水上升规律如图1所示，单重介质传导型裂缝模型与局部网格加密法及双重介质模型计算结果基本一致，在计算结果上都能达到对真实裂缝的达西渗流规律的描述。而实际应用中常用的等效级差法见水晚且后期含水高，其中小网格的计算效果又好于大网格。
从计算效果、计算速度、收敛性、直观性等方面综合对比以上几种方法：
（1）局部网格加密法是对裂缝最直接的描述，且能显式地表现出油水的流动，计算结果最能反映真实的地下达西流动。但当裂缝尺度较小时，模型收敛性差；当裂缝条数多时难以操作。
（2）等效级差法优点是操作简单。但使用基质网格来等效裂缝增加了裂缝的储油能力，计算结果欠佳。但网格较小时也可以做一定程度的近似。
（3）单重介质传导型裂缝模型使用方便，且可以达到与局部网格加密法近似的计算效果，计算速度比双重介质模型要快。但在对裂缝的显示及油水运动描述方面存在不足。
（4）双重介质模型仍然是计算效果最好的模型，尤其适用于裂缝较多，难以逐条描述的情况，对大裂缝的计算效果也不比其它方式逊色。但在裂缝的显示及油水运动描述方面表现欠佳，且计算速度较慢。
上述几种方法仍存在着共同的缺点：
（1）仍然建立在达西流的基础上，因此对于大裂缝类管流的流动特征难以反映。
（2）对于裂缝形态的描述困难，对于不同方向的多条裂缝难以处理，更适用于简单的压裂裂缝。
4 结论
（1）在大裂缝的模拟中，单重介质传导型裂缝模型与局部网格加密法及双重介质模型可以反映真实裂缝的达西渗流规律，其中单重介质传导型裂缝模型计算综合性能最优，可操作性最强。
（2）以上几种方法在大裂缝高速非达西流描述、裂缝的形态描述方面还存在共同的缺陷，将是新一代裂缝性油藏数值模拟软件需要解决的问题。
参考文献
Abstract: In second-level location of 8 degree area, it has designed a prestressed concrete frame industry workshop wit single span 5 layer based on current regulations. At the same time, we analyze the earthquake response of the structure under major earthquake with three kinds of seismic wave. The results show that the industry workshop has a good seismic behavior in a given seismic conditions.
关键词：预应力混凝土框架；数值模拟；罕遇地震
Key words: prestressed concrete frame；numerical simulation；rare earthquake
中图分类号：TU37 文献标识码：A 文章编号：1006-4311（2011）28-0096-02
1 工程概述
该工程为单跨5层工业厂房，跨度15m，底层层高4.5m，其它层高4.2m，柱距7.5m。采用预应力混凝土框架结构，图1为结构的平面布置详图。设计标准为：8度抗震设防，二级抗震等级及II类场地[2-3]。选用的是普通II级钢筋，箍筋采用的是I级钢。梁柱采用C40混凝土。图2为预应力筋的布置，图3为结构配筋的示意图。
2 数值模型建立途径
2.1 模型基本假定 数值中，框架梁、柱结构符合平截面假定；钢筋和混凝土粘结可靠；结构构件受弯破坏先于受剪破坏；节点不破坏；不考虑混凝土的收缩和徐变影响。
参考文献：
[1]建筑结构设计规范（GB50011-2001）[S].北京建筑工业出版社，2001.
[2]耿耀明，黄鼎业.预应力混凝土框架结构抗震能力的静力弹塑性分析方法.工业建筑，2003，33（9）：29-32.
[3]混凝土结构设计规范（GB50011-2001）[S].北京建筑工业出版社，2001.
[4]熊学玉，李春祥，耿耀明等.大跨预应力混凝土框架结构的静力弹塑性分析.地震工程与工程振动，2004，24（1）：68-75.
Abstract： Lattice Boltzmann Method （LBM） has been used to simulate the seepage field of the fractured rock masss， which is born with the advantages of parallel characteristics， simple boundary conditions， easy to implement the program， clear physical image compared with the traditional numerical simulation. JRC numerical generation method is used to establish the fracture structure of two-dimensional rock mass with JRC of the joint roughness coefficient in a certain range. Based on the Lattice Boltzmann Method （LBM）， the the boundary conditions are that the inlet and outlet are set to the non-equilibrium extrapolation scheme， and the crack surface of the upper and lower rock mass is set to the standard rebound format. C++ programming can realize the Lattice Boltzmann Method simulation and verify the classic Poiseuille flow. In the rough fracture with a pressure difference of 0.002 and the width ratio of 400：21， the flow pattern of the rock mass seepage is linear sub-cubic flow， and with the increase of the JRC value， the deviation of the cubic value from the value of 3 is more obvious. In the case of JRC=14.1， the the fracture seepage can be represented by the sub-cubic seepage， cubic seepage and super-cubic seepage under different pressure difference and different fracture width.
关键词：格子Boltzmann方法（LBM）；数值模拟；吻合岩体单裂隙渗流；立方定律
Key words： Lattice Boltzmann Method （LBM）；numerical simulation；water flow in the single fractured rock mass with a close correlation；cubic law
中D分类号：TU45 文献标识码：A 文章编号：1006-4311（2017）10-0094-04
0 引言
岩体是指在一定范围内的自然地质体，它经历了漫长的自然历史过程，经受了各种地质作用，并在地应力作用下在其内部形成了各种地质构造行迹，如节理、褶皱、劈理、层理等，正是由于这些岩体裂隙的存在，为地下水的存储和流动提供了场地，由此导致的地质灾害及渗流问题在裂隙型油气田的开发、深基坑和隧道的开挖等方面得到重视并逐渐变成一个热点问题。因此，研究岩体渗流问题具有非常大的工程意义和工程应用价值。
从20世纪70年代开始W.Wittke用有限元法计算裂隙中的水流、Krizek等用有限差分法计算了各种裂隙网络系统内的流势形态。相应的，已有的计算流体力学方法也可以分为微观方法、介观方法和宏观方法三类[1]。介观方法中，流体被视为离散成一系列的流体粒子。最常见的介观模拟方法是格子Boltzmann方法。本文基于格子Boltzmann方法，在平直裂隙面下验证了该方法的可靠性，同时还分析研究了岩体节理裂隙粗糙度系数JRC、隙宽、压力差与渗流特性间的相互联系。
1 立方定律基本理论
由于理想的平板在自然界岩体裂隙中是不存在的，岩体裂隙往往是粗糙凹凸不平的，所以在光滑平板下推得的立方定理在粗糙裂隙渗流中不完全适用，所以要在原来的公式中做一些修正。依照上式我们已经知道隙宽的指数n=3时称之立方定理，为了让n=3命名为立方相应，本文将n3命名为超立方。前苏联学者Lomize[11]、Amadei及速宝玉[12]等用试验研究，认为处于层流状态的天然裂隙渗流依然满足立方定理，仅需要对系数进行修正。目前，在实际工程应用和科学研究中，使用该定理的比较多。
2 格子Boltzmann模型
格子Boltzmann方法（Lattice Boltzmann method）诞生至今已有28年，格子Boltzmann方法以介观层次为视角，建立将流体离散成流体粒子、物理区域离散成一系列的格子、时间离散成一系列的时步的模型。一个完整的格子Boltzmann模型一般由三部分组成：格子，即离散速度模型；平衡态分布函数；分布函数的演化方程。
2.1 格子Boltzmann方法的理论
本文选用的D2Q9模型中的离散速度如图2。
3 讨论
3.1 二维平板模型中的泊肃叶流的模拟验证
为了验证该模型的可靠性，本文用D2Q9做了一个简单的平板模型：长为200格子单位，宽为50格子单位。进口压力P-in=1.001；出口压力P-out=0.999；雷诺数Re=20；tau=1.2；niu=0.125。
为了更为直观地看到平板流中流场的信息，图4给出了计算收敛的数值解成像后的流场信息图。
为了使取值更加精确，本文取平板中间部位（x=100）上的横截面水平速度值和Poiseuille理论值进行对比，如图5所示，可以很明显地看到两者吻合度非常高。
3.2 压力差、裂隙开度与n值的关系
从图6可以看出当压力差值ΔP=0.31、0.41时，刚开始当隙宽为11时n值均大于3为超立方流，随着隙宽的增加n值逐渐减小由立方流转变为次立方渗流；当压力差ΔP=0.21、0.11时，n值随隙宽先增大后逐渐减小，但n值始终未超过3，所以均为次立方渗流；当压力差ΔP=0.01、0.001时，n值随隙宽的增加逐渐增加，最后逐渐趋于平缓，并且n值始终未超过3，渗流表现为次立方渗流特性。
4 结论
本文使用C++编程技术，数值模拟生成了吻合的岩体单裂隙粗糙面。在LB方法的基础上，采用D2Q9模型，模拟了平板流并且研究了在压力作用驱动下的单裂隙的渗流特性，得到了以下结论：
①LB方法模拟泊肃叶流数值解和理论解误差非常小，可以很好地模拟泊肃叶流。
②利用LB方法来模拟平板流是可靠的，LB方法适用于本文所要用来模拟的岩体裂隙渗流模型。
③可以将裂隙开度和压力差值作为次立方渗流和超立方渗流的判别标准，即当压力差值ΔP≥0.31且隙宽小于等于31时裂隙渗流为线性超立方渗流，当压力差值ΔP≥0.31且隙宽大于31时裂隙渗流为线性立方渗流或次立方渗流；当压力差值ΔP
参考文献：
[1]王志良，申林方，李邵军.岩体裂隙面粗糙度对其渗流特性的影响研究[J].岩土工程学报，2016：1-7.
[2]盛金昌，速宝玉.裂隙岩体渗流应力耦合研究综述[J].岩土力学，1998，19（2）：92-98.
[3]LIU J， ELSWORTH D， BRADYB H. Linking stress-dependent effective porosity and hydraulic conductivity fields to RMR[J]. International Journal of Rock Mechaincs and Mining Sciences，1999，36（5）：581-589.
[4]AMERICAN， BARTON N. Review of a new shear strength criterion for rock joints[J]. Engineering Geology， 1973， 7（4）：287-332.
[5]李化，黄润秋.岩石结构面粗糙度系数JRC定量确定方法研究[J].岩石力学与工程学报，2014，33（2）：3489-3497.
[6]Hyun-Sic Jang， Seong-Seung Kang， Bo-An Jang. Determination of Joint Roughness Coefficients Using Roughness Parameters[J].Rock Mech Rock Eng （2014） 47：2061-2073.
[7]Tse R，Cruden DM.Estimating joint roughnesscoefficients[J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts . 1979， 16（2）：303-307.
[8]Maerz NH， Franklin JA， Bennett CP.Joint roughness measurement using shadow profilometry[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts . 1990，27（5）： 329-343.
[9]Clarke，Keith putation of the FRACTAL DIMENSION OF TOPOGRAPHIC SURFACES USING THE TRIANGULAR PRISM SURFACE AREA METHOD [J]. Computers and Geosciences. 1986（12）：713-722.
[10]Turk N，Greig M.J，Dearman W.R，Amin F，F. Characterization of rock jointsurfaces by fractal dimension[C]//The28th U.S. Symposium on Rock Mechanics （USRMS）. Rotterdam：A.A. Balkema，1987.
[11]Lomize G M. Flow in Fractured Rock （in Russian ）[M].Gosemergaizdat， Moscaw， 1951： 127-129.
微分方程主要采用隐式方法,如:隐式RK方法,BDF方法,IRK方法等。而采用隐式方法将刚性方程离散化以后,其变为线性或非线性方程(组)的求解问题。目前,对线性或非线性方程(组)的求解,多采用Newton-Raphson迭代求解。但对于某些非线性方程组,由于方程之间的非线性化程度相差较大,采用Newton-Raphson迭代方法数值求解的结果并不理想。本文利用Brown算法求解此类非线性刚性系统,具有较高精度和较快迭代速度的优点,数值试验结果表明了该方法的有效性。
2 Brown算法
考虑多个实变量的非线性方程组
(2.1)
的数值求解问题,非线性方程组可以用向量形式表示:,其中,。
形如:的公式称为Newton-Raphson迭代公式。由于该方法是将,同时线性化,所以它并未考虑充分利用的具体结构。如果一个非线性的向量函数,其线性精度在各个分量,上的分布可能是不平衡的,有的分量是非线性函数,而有的分量是线性函数,同时非线性函数组中也有非线性程度高低的差别,在此情况下,利用Newton-Raphson迭代方法对所有分量采用完全相同的数值处理,不利于方法整体计算效率的提高。
针对以上情况,Brown于1969年提出了按分量函数方程,来形成迭代过程[3],其基本思想是对各分量逐个线性化并用其中每一个线性方程消去余下非线性方程中的一个变量,最后整个方程组就简化为一个仅含单个变量的非线性方程,应用一次单步Newton-Raphson迭代并结合逐一回代,即完成一次迭代过程[4]。
Brown算法的迭代步骤如下:
第一步,设为方程组(2.1)解的第次近似,函数在处近似用线性函数
替代,令,由此求出:
定义上式右端为。
第二步,对函数定义一个新函数Brown算法,且记,其中。类似地,用线性函数来近似替代。令,解出,
此时,为个变量的线性函数,并记此线性函数为。
第步,由线性函数,可得,利用Newton-Raphson迭代,求得,并由出发,利用逐一回代,即
(2.2)
从而可求出,至此完成了一次Brown迭代过程。
3 数值试验
考虑以下常微分方程组初值问题:
问题1
其中:;。
问题2
其中:;。
对于上述两问题,当时,可计算其右函数组的Jacobi矩阵的特征值,均有,其余特征值绝对值均不超过6,因此系统呈强刚性。此外,观察两问题中的右函数组,可以看出除最后一个函数是高度非线性化外,其余函数都是线性的。
对于上述两问题,采用隐式Euler方法离散方程组,并分别用Newton-Raphson迭代法与Brown迭代法求解,取步长,及相对误差界(表示迭代次数)控制每步迭代,最后得到数值解的最大绝对误差界,方程真解为:问题1,,,;问题2,,,。计算结果对比分析如表1所示。
表1 数值计算结果
4、结束语
对于实际问题中的刚性系统离散化后,如果非线性方程组的线性化程度不同,Brown迭代求解比Newton-Raphson迭代法具有较大的优势,另外需要指出的是在实际运算中,方程应预先进行排列,将线性方程放置在最前,再次为非线性化程度由低到高排列,可以有效的提高运算效率。
参考文献:
[1]李寿佛.刚性微分方程算法理论[M]. 长沙: 湖南科学技术出版社,1997.
[2]刘德贵,费景高.刚性大系统数字仿真方法[M].郑州: 河南科学技术出版社,1994.
关键词：水下航行体；变水域；数值模拟计算
0 引言
本文拟通过潜艇以相同的速度由无限水域进入不同阻塞比的有限水域和以不同的速度进入同一阻塞比的有限水域两种方法来分析在不同参数影响下潜艇进入变水域的阻力性能的规律（阻塞比β=潜艇面积/隧道面积）。通过几组参数计算结果的对比，揭示水下航行体进入变水域阻力性能与速度和阻塞比的关系，并通过此关系验证动网格技术对变域情况数值模拟的有效性，可作为后续导弹出筒等相似计算的方法依据。
1 基本理论与数值方法
1.1雷诺时均纳维叶-斯托克斯方程
雷诺时均纳维叶-斯托克斯方程和雷诺时均连续性方程：
―湍流动能生成项。其中的一些常数值如表1。
2 计算结果与分析
2.1 动网格模型
动网格区域与主体部分建模是分离的两个部分，而它们之间的数据交换是通过在动网格的外表面区域与主体部分与动网格交界内部区域之间完成的。
2.2 计算结果分析
经过数值模拟计算，可得到如图1数据。
从图1计算数据分析可以得出，潜艇表面压力在前体与后体变化比较剧烈，潜艇头部尖点处压力达到最大值，在前体与平行中体和后体与平行中体交接处变化最大，平行中体上压力值基本保持不变。在潜艇进入受限域的过程中，压力差值逐渐变大，在潜艇完全进入受限域也就是20s后压力差达到最大值并保持稳定。因此，后续计算数据主要以潜艇完全进入受限域后气动性能数据为主。
2.3 阻塞比的影响
阻塞比为潜艇横截面积与受限域横截面积的比值。当潜艇横截面积一定的时候，受限域横截面积越大，阻塞比越小。阻塞比对潜艇通过受限域时的水动力效应有很大的影响。为研究阻塞比影响，分别对阻塞比为0.08、0.12和0.2的数值模型进行了模拟计算，计算结果如图2。
在受限域内，随着阻塞比的增大，受限域内的压力值也增大。潜艇表面所受的压力也随之减小，在前体与后体不是很明显，在平行中体上受到阻塞比的影响相对较明显。在潜艇头部到达受限域至潜艇完全进入受限域的时候，随着阻塞比的增大，压力也增大。但是在潜艇完全进入受限域后情况却正好相反，阻塞比大的受限域潜艇表面最大压力值反而更小。随着阻塞比的增大，在进入受限域时潜艇阻力变化值也增大，完全进入受限域后，阻塞比对潜艇阻力的影响很小。
2.4 速度的影响
速度增大也会增大潜艇表面压力差值，压差值的增大会在潜艇表面产生比较大的扭矩，这对潜艇的结构安全性有很大的威胁。为研究速度对潜艇压力的影响，分别对速度为10m/s、15m/s和20m/s的潜艇进入受限域进行了数值模拟，得到结果如3所示。
潜艇速度对潜艇所受的压力最大值的影响特别大。随着潜艇的速度的增大，潜艇所受压力最大值的值也随之明显增大，并且潜艇受到的压力最大值的改变幅度也随之增大。也就是说，随着潜艇速度的增大，潜艇在受力最大值点的压力差值也随之增大。潜艇的速度值越大，潜艇在无限域中受到的阻力越大，在进入受限域过程中阻力差的值也越大。
3 结论
本文基于无粘、不可压缩流紊态流场的雷诺时均方程和k-？两方程紊流模型模拟了潜艇通过变水域时的水动力性能，建立了潜艇-变水域水动问题数值计算模型，应用有限体积法进行求解。对水下航行体在动网格计算法下得出的水动力性能进行了分析，得出了以下结论。
（1）利用CFD软件FLUENT建立了潜艇通过变水域的水动力学模型，并通过计算分析得出了潜艇由无限域进入受限域的水动力性能。
（2）随着阻塞比的增大，受限域内压力值也随之增大，但是潜艇表面压力值却是随之减小，在平行中体最明显；潜艇表面最大压力值的压力差值也增大，潜艇在进入受限域过程中总阻力值的最大值与阻力差值也随之增大，但是在无限域中和完全进入受限域后的总阻力值变化很小。
关键词：黄土隧道；数值模拟；沉降；变形
中图分类号：U455.4 文献标志码：B
0 引 言
随着公路建设的飞速发展，在路线车道数的选择中，以多车道代替双车道对于缓解交通压力更为有利，于是出现了很多大跨度隧道[15]。
相对于普通隧道的施工方法来说，大跨度隧道的施工因其断面大、跨度大、扁平率更小，所以隧道的受力更为复杂。目前，随着大跨度隧道逐渐增多[67]，对大跨度隧道施工方法的研究受到了业界人士的关注。主要的施工方法有两种，即三台阶七步开挖法和双侧壁导坑法，而对于这两种方法的特点与适用性的研究相对较少。为此，本文以神府高速墩梁大跨度黄土隧道工程为依托，对大跨度黄土隧道的两种施工方法进行对比研究，以确定哪种方法最优。
1 工程概况
墩梁隧道衬砌断面内轮廓采用三心圆方案。上行线起止桩号为RK25+705～RK27+120，全长1 415 m；下行线起止桩号为RK25+705～RK27+120，全长1 328 m，总计全长2 743 m。整个隧道均为V级围岩，隧道支护体系结构均为复合式衬砌，二次衬砌拱部厚度为06 m。深埋段开挖高度为1216 m，宽度为1726 m，面积为1652 m2。浅埋段开挖高度为1219 m，宽度为 1731 m，面积为1704 m2。设计时速为80 km·h-1，建筑限界有效净宽为1425 m，净高为52 m。
2 施工方法对比分析
2.1 三台阶七步开挖法
三台阶七步开挖法的施工步骤如图1所示，具体如下。
（1） 上部弧形导坑开挖。环形开挖上部弧形导坑应在拱部超前支护后进行，预留核心土宽度宜为隧道开挖宽度的1/3～1/2，核心土长度宜为3～5 m。根据初期支护钢架间距来确定开挖循环进尺，应控制在15 m以内，上台阶开挖矢跨比应大于03，开挖后立即初喷3～5 cm混凝土，并及时进行喷、网系统支护，架设钢架，在钢架拱脚以上60 cm高度处，按下倾角60°紧贴钢架两侧边沿打设锁脚锚杆，钢架与锁脚锚杆牢固焊接，最后复喷混凝土至设计厚度。
（2） 左、右侧阶开挖。根据初期支护钢架间距来确定开挖循环进尺，应控制在15 m以内，左、右台阶错开2～3 m，开挖高度一般为3～35 m。
（3） 左、右侧下台阶开挖。左、右侧下台阶开挖工序与左、右侧阶开挖工序一致。
（4） 上、中、下台阶预留核心土开挖。各台阶分别开挖预留的核心土，开挖进尺应与各台阶循环进尺一致。
（5） 隧底开挖。每循环开挖长度宜为2～3 m，开挖后及时施作仰拱初期支护，完成隧底开挖、支护循环后，及时分段施作长度为4～6 m的仰拱。
2.2 双侧壁导坑法
双侧壁导坑法的施工步骤如图2所示，具体如下。
（1） 左导洞（先行导洞）上台阶开挖。采用正台阶法开挖，每次开挖进尺0.5～0.75 m，开挖后立即初喷3～5 cm混凝土，及时架设钢筋网和型钢拱架，及时安设边墙锚杆和锁脚锚杆，随后喷射混凝土至设计厚度，以便形成较为稳定的支护体系。
（2） 左导洞下台阶开挖。导洞下半断面与上半断面前后错进15～20 m，开挖左导洞下台阶的临时支撑和左边墙（临时支撑比左边墙晚进一榀拱架），每次进尺1.0～1.5 m，这样可以避免同一断面的开挖面同时暴露。
（3） 右导洞上台阶开挖。采用正台阶法开挖，开挖完成后的施工工序与左导洞上台阶开挖工序一致。
（4） 右导洞下台阶开挖。分别开挖右导洞下台阶的临时支撑和右边墙（临时支撑比右边墙晚进一榀拱架），每次进尺1.0～1.5 m，这样可以避免同一断面的开挖面同时暴露。
（5） 主导洞上台阶开挖。主导洞拱部开挖前先测量画出开挖轮廓线，延开挖轮廓线打设超前小导管注浆加固，开挖后钢筋网和型钢拱架及时架设，边墙锚杆和锁脚锚杆及时安设，随后喷射混凝土至设计厚度，以便形成较为稳定的支护体系。
（6） 主导洞下台阶及仰拱开挖。开挖后型钢拱架及时架设，随后初喷混凝土至设计厚度，以便形成较为稳定的支护体系。
2.3 两种工法对比分析
（1） 在资源配置相同的情况下，三台阶七步开挖法施工每循环进行两榀，一榀0.75 m，26 h可以完成两个循环，每月最高进度指标为80 m，平均每月进度为76 m；双侧壁导坑法施工每循环0.6 m，每天可完成2个循环，每月最高进度指标为42 m，平均每月进度为36 m。
（2） 与双侧壁导坑法相比，三台阶七步流水法开挖施工技术可节省大量临时钢架施工支护费用，减少了投资。
（3） 根据现场实测数据，墩梁隧道在围岩级别相同的条件下，采用三台阶七步开挖法开挖拱部下沉值为90～120 mm，而采用双侧壁导坑法施工拱部下沉值仅为 30～50 mm。
3 数值模拟及分析
根据墩梁隧道设计情况，隧道总高1216 m， 净宽1726 m。将数值计算模型解析区域设为：竖向向下取隧道高度的5倍，左、右各取隧道直径的5倍，向上取至原地面线。在所取范围之外可认为不受开挖等施工因素的影响，即在这些边界处可忽略开挖等施工所产生的位移和应力。同时，保证模型不出现刚体转动及位移。
因隧道属于深埋段且无构造节理影响，故按初始自重应力场来考虑地应力。为了使得计算结果更可靠，在隧道内及其周围采取细密网络划分，采用4节点结构单元K线面映射网络划分模型网络。计算的主要内容为各个施工步隧道的拱部下沉与水平收敛。
3.1 三台阶七步开挖法数值模拟
对墩梁隧道采用三台阶七步开挖法进行数值模拟，其局部模型网格划分如图3所示。
（1） 净空收敛。初期支护封闭后隧道净空收敛计算结果如图4所示。
（2） 拱部下沉。初期支护封闭后隧道拱部下沉值计算结果如图5所示。
通过有限元数值模拟，在围岩级别等各种施工环境相同的条件下，采用双侧壁导坑法有效控制了隧道的下沉值与净空收敛值，而三台阶七步开挖法下沉值与净空收敛值均较大，说明在大跨度黄土隧道中采用双侧壁导坑法对于控制围岩变形更加有效。
4 结 语
（1） 侧壁导坑法是黄土隧道施工中最能有效控制沉降的工法，适用于松散易塌的软弱土层地段，这方面控制的重点是仰拱的及时封闭成环和左右侧壁永久支护拱架的对应，但由于其工序繁杂，故施工难度大、工期长、造价高。
（2） 三台阶七步开挖法主要适用于土体节理水平、稳定性较好的土层地段，其优点是工序流水作业能有效提高工效，加快施工进度，缺点是不能有效控制沉降。
（3） 有限元数值模拟结果表明，在围岩级别等各种施工环境相同的条件下，采用双侧壁导坑法有效控制了隧道的下沉值与净空收敛值，而三台阶七步开挖法下沉值与净空收敛值均较大，说明在大跨度黄土隧道中采用双侧壁导坑法对于控制围岩变形更加有效。
参考文献：
[1] 林镇洪.单边扩帮特大跨公路隧道开挖稳定性分析[J].铁道科学与工程学报，2009，6（4）：4650.
[2] 夏保祥，程崇国.三车道大断面公路隧道研究现状综述[J].地下空间，2002，22（4）：360366.
[3] 王梦恕，谭忠盛.中国隧道及地下工程修建技术[J].中国工程科学，2010，12（12）：410.
[4] 何开伟.超大跨隧洞洞口破碎围岩段施工[J].西部探矿工程，2007（2）：126128.
[5] 李国良.大跨黄土隧道设计与安全施工对策[J].现代隧道技术，2008，45（1）：5362.
