时间：2023-06-07 09:09:10
开篇：写作不仅是一种记录，更是一种创造，它让我们能够捕捉那些稍纵即逝的灵感，将它们永久地定格在纸上。下面是小编精心整理的12篇数值分析，希望这些内容能成为您创作过程中的良师益友，陪伴您不断探索和进步。
在计算机上运用数值分析解决实际问题的过程为：实际问题→数学模型→数值计算方法→程序设计→上级计算求出结果[1]。数值实验的设计应当充分体现这一过程，同时也应当充分体现数学建模思想. 而目前地方高校在数值分析课程的教学与实验中，普遍存在重理论、轻实践、纯数学十足的问题，部分即使避免这些问题，做到利用计算机进行可视化教学和算法编程实践，但学生也缺乏分析问题、解决问题的能力，因为学生的算法与程序多是百度而来，缺乏思考. 因此，数值分析课程的数值实验因当因时而异，因专业而异，紧跟时展，做到充分吸引学生的注意力，激发学生的兴趣，调动学生积极性，让学生积极主动的去想办法解决问题，而非被动的接受。这样的数值实验才能有效培养学生应用能力。
数值分析的内容广泛，包含插值拟合、数值积分与数值微分、数值代数、微分方程数值解法、非线性方程与方程组的数值解法，最优化等等。下面以插值拟合、微分方程数值解法、非线性方程数值解法、最优化的相关实际问题为例，研究数值实验的设计。
2 结束语
通过以上案例发现：如果知道解决实际问题的数值分析方法，利用MATLAB能很方便的求出实际问题的结果。但是对于多数学生而言，由实际问题无法得到数学模型，更无法知道相应的数值分析方法。因此在学生有一定数值分析基础后，才能引入这些实际问题的数值实验；同时选取的数值实验必须为学生深入研究预留了空间，因为在讲解相关的知识背景，详细分析问题，建立模型后，对学生进行的是分层次指导解决这些问题（即或利用MATLAB函数命令简单编程计算，或利用MATLAB的专用工具箱计算，或设计算法流程利用MATLAB编程计算）；这样每一个学生都能参与到实验中来，每一个学生在实验中都有收获；最后撰写实验报告，阐明实验的目的、要求、过程、收获等.如此，通过这些数值实验，学生既自己动手解决实际问题，培养了其综合应用能力，又让其充分体会到数值分析的魅力，为大学生数学建模竞赛和后续专业学习打下坚实基础.当然，随着技术的进步，不同专业的数值分析应用实验需要教师不停的探索和研究。
参考文献：
[1] 徐翠薇，孙绳武.计算方法引论[M]. 3版.北京：高等教育出版社，2007.
[2] 吕国英.算法设计与分析[M]. 2版.北京：清华大学出版社，2009.
[3] 张德丰等. MATLAB数值分析[M]. 2版.北京：机械工业出版社，2012.
[4] 杜廷松.关于《数值分析》课程教学改革研究的综述和思考[J]. 大学数学，2007，23（2）：8-15.
[5] 张涛，. 数值实验在《数值分析》教学中的重要性分析[J]. 长江大学学报：自然科学版，2009，6（1）：371-372.
1计算方法
1．1桨叶／机身非定常面元
1．1．1面元基本原理除物面附近及尾流区外，旋翼流场可假设为无粘、无旋、不可压。在惯性坐标系下，连续方程可表示成速度势的函数［17］，即式中SB与SW分别为物面（桨叶或机身）和尾迹涡面，n为物面外法线单位矢量，r＝（x，y，z）为空间点位置。
1．1．2边界条件物面边界条件要求相对于物面的法向速度为0，远场边界条件要求物体对流体的扰动在无限远处为0，即假设物体表面由N个面元组成，尾迹涡面由Nw个面元组成，采用等强度四边形偶极子面元，则式（2）可表示成如下：
1．1．3面元压力旋翼流场确定之后，可根据非定常Bernoulli方程，通过速度势和物面速度计算压力分布。非定常项／t可通过求解物体表面速度势得到。对于机身，非定项主要来源于桨叶和旋翼尾迹的影响。桨叶影响可通过速度势直接求解，而尾迹影响为尾迹对面元的诱导速度与尾迹自身速度之积［19］。
1．2时间步进自由尾迹为求解桨叶和机身面元强度分布，在解式（5）或式（8）之前需计算旋翼尾迹。本文采用时间步进自由尾迹［12－14］。时间步进自由尾迹基于不可压假设，并把旋翼尾迹漩涡简化为直线涡线。旋翼涡量场可由三维不可压粘性Navier－Stokes方程描述，表示成速度－涡量采用有限差分近似时间和空间导数求解式（14）。涡线位置由时间步进格式求解得到。文中采用二阶精度的预估－修正格式（PC2B）［12－14］。
1．3旋翼桨叶运动方程旋翼尾迹和桨叶面元汇／偶极子分布与桨叶的挥舞运动方程紧密相连，因此在描述旋翼尾迹时需求解桨叶的挥舞运动。根据桨叶挥舞铰力矩为零建立刚性桨叶挥舞运动方程。桨叶挥舞运动可表示成一组常微分方程，并采用四阶Runge－Kutta求解［14］。
1．4桨叶面元／尾迹耦合为计算旋翼尾迹的畸变效应，采用全展涡线代替偶极子面元尾迹。旋翼桨叶由非定常面元构成，桨叶脱出的尾随涡由尾随偶极子面元构成，旋翼尾迹则由连接于尾随涡的全展涡线构成，并从桨叶尾随偶极子面元中脱出（如图1）。基于桨叶后缘Kutta条件及尾迹偶极子面元强度与涡线涡量强度等价原则，建立桨叶面元与尾迹之间的联系。偶极子面元与涡线等价原则可表示成如下。在各时间步，旋翼尾迹涡线强度由桨叶面元强度决定，同时桨叶面元的汇／偶极子强度又与旋翼尾迹涡线有关，由此确保桨叶非定常面元与旋翼尾迹的紧密耦合。
1．5旋翼尾迹／机身干扰低速前飞状态下，机身浸润在旋翼尾迹中，因此旋翼尾迹涡线将靠近机身表面。由于机身的阻塞效应，旋翼尾迹涡线靠近机身表面的速度减小，而切线速度增加，此时机身非定常压力主要来源于旋翼尾迹。由式（11）可知，非定常项由尾迹涡线移动速度和尾迹诱导速度构成，因此旋翼尾迹几何特性对旋翼尾迹／机身干扰影响显著。由于机身表面载荷与旋翼尾迹几何密切相关，因此旋翼尾迹涡线靠近机身表面的运动特性就显得非常重要［22］。为满足机身表面无穿透条件，并模拟涡线靠近机身表面的加速现象，文中采用涡线镜面法。与二维点涡镜面类似［22］，尾迹涡线由两点直线构成，因此可通过涡线中点的矢量镜面得到镜像涡线，镜像涡线涡量为Γ′＝－Γ（如图2）。在各时间步，通过桨叶和机身非定常面元同步求解，得到桨叶和机身的非定常气动力，而后推进旋翼尾迹，由此计算旋翼／机身非定常气动干扰。
2计算结果与分析
为验证本文旋翼／机身非定常气动干扰分析方法的准确性，文中将计算前飞状态的Maryland、ROB－IN旋翼／机身干扰下的机身非定常压力分布，并与可得到的实验值、CFD计算结果对比验证。随后分析前飞速度、旋翼与机身高度对非定常气动干扰的影响。
2．1Maryland旋翼／机身干扰本算例为前飞状态下的Maryland旋翼／机身干扰试验［23］，旋翼系统由4片直径为1．65m的矩形桨叶铰接构成，桨叶线性负扭为－12°，翼型为NASARC310和RC410，弦长为0．0635m，旋翼转速为1860rpm，机身长度为1．94m，机身最大截面直径为0．254m，机身尾梁与机身截面直径之比为1∶2．5。桨毂中心与机身重心高度为0．24m。机身压力传感器分布如图。旋翼／机身干扰下的各传感器非定压力随桨叶方位角变化历程如图4。从图4中可以看出，本文计算方法计算得到Maryland旋翼／机身干扰下的非定常压力时间变化历程与实验测量值吻合较好。图4中各传感器非定常压力随方位角的变化均表现为4Ω周期波动，此倍频与旋翼桨叶片数相同，由此说明桨叶通过机身上方所产生的显著非定常干扰效应。机身头部（传感器1）非定常压力呈现出类正弦波动，主要原因为传感器1在旋翼下方，受桨叶通过性影响显著。传感器9、10在旋翼下方之外的尾梁，受到桨叶通过性影响减小，而主要受到旋翼尾迹与尾梁干扰影响，因此非定常压力呈现锯齿形状。传感器9比传感器10更靠近旋翼，因此传感器9的非定常压力受到桨叶通过性影响更显著，表现的类正弦特性更显著。传感器11、12在尾梁左右两侧，主要受到旋翼尾迹／机身干扰影响，表现出锯齿形状。从图4（b）、（c）中可以看出，传感器9的压力峰值相位超前于传感器10，主要原因为尾梁传感器9比传感器10更靠前，旋翼尾迹将先靠近传感器9。但随着诱导速度的向下作用，尾迹距传感器10的距离更小，因此负压峰值更大，旋翼尾迹／尾梁干扰更显著。从图4（d）、（e）中可以看出，尾梁左侧传感12的非定常负压峰值大于尾梁右侧传感器11，主要原因为旋翼右旋，旋翼尾迹贴近尾梁的左侧，因此对传感器12的干扰作用大于右侧的传感器11。
2．2ROBIN旋翼／机身干扰本算例为前飞状态下的ROBIN旋翼机身干扰试验［24］。旋翼系统为2MRTS（2－MeterRotorTestSystem）缩比旋翼［25］，机身为流线型机身。2MRTS旋翼由4片矩形桨叶铰接构成，桨叶半径为0．861m，弦长为0．0663m，线性负扭为－8．0°，翼型为NACA0012翼型，旋翼转速为2000r／min，前进比为μ＝0．151。机身长度为1．999m，桨毂与机身重心垂直距离为0．322m，旋翼周期变距由风洞试验据得到［25］。各片桨叶由弦向60段和展向20段面元组成，旋翼系统共由4800个面元构成，ROBIN机身由10842个面元组成。机身头部、发动机舱、尾梁、机身左右两侧压力传感器分布如图5所示。ROBIN旋翼／机身干扰下的旋翼尾迹如图6所示，从图中可以看出，旋翼左右两边形成比较明显的桨尖涡。由于机身的排斥作用，旋翼中间尾迹向上、左右两侧移动，但尾梁后段仍然浸入在旋翼尾迹中，因此将产生显著的旋翼／机身干扰。桨叶／机身压力分布如图7所示，机身头部和发动机舱前、后缘部分产生较大压力。前飞状态下，由于旋翼前行边和后行边桨叶相对来流的非对称，因此需通过周期变距改变桨叶的桨距以保证整机左右平衡，并由此导致前行边和后行边桨叶气动环境、桨叶脱出涡量不一致，从而引起旋翼尾迹的非对称。由于旋翼尾迹非对称和桨叶位置的变化，导致机身前后、左右两侧压强非对称，由此产生时变载荷。ROBIN机头顶部、发动机顶部、尾梁顶部、机身两侧非定常压力随桨叶方位角变化如图8。从图8中可以看出，本文计算方法计算得到ROBIN旋翼／机身干扰下的各传感器非定常压力间变化历程与实验值［24］和CFD计算结果［4，26］比较吻合，由此验证本文计算方法的可靠性。由于旋翼系统采用4片桨叶，因此图8中各传感器非定常压力随方位角变化历程均呈现4Ω的周期特性。从图5可以看出，各传感器均位于旋翼下方，因此各传感器非定常压力主要表现为桨叶通过性影响。为反映机身各处压力的变化特性，机身各传感器非定常压力峰值相位和幅值如表1．从表1可以看出，机身头部传感器6非定常压力峰值出现在桨叶通过机身后，而发动机舱传感器22和尾梁顶部传感器15非定常压力峰值出现在桨叶未通过机身前，主要原因为机身头部距离桨尖平面距离更大，阻塞效应较小，且存在阻塞滞后，并由此导致机身头部非定压力幅值小于尾梁顶部。由于传感器22处于桨根下方（图5），因此受桨叶通过性的影响小于传感器15。由于桨叶右旋转，导致机身左侧流场阻塞，机身右侧流场扩展，由此导致左侧传感器13的压力幅值大于右侧传感器19，且右侧峰值相位滞后与左侧。由于机身头部、尾梁、机身左、右侧非定常压力幅值与相位的差异导致机身力与力矩的非对称，由此产生4Ω周期激励载荷。
2．3前飞速度对旋翼／机身干扰影响以ROBIN旋翼／机身干扰为基本算例，分析前飞速度对机身非定常压力的影响。从图9中可以看出，随着前飞速度的增加，旋翼载荷增加，桨叶通过性对机身非定常压力影响增加，由此导致机身头部、发动机舱、机身右侧压力幅值均增加，且随着前飞速度增加，压力幅值增加速率增加。但尾梁传感器15压力幅值随前飞速度的增加而先加后减小，原因为前飞速度较小时，尾迹／尾梁干扰显著，速度增加导致尾迹对尾梁的诱导非定常项影响增加，因此压力幅值增加，但随前飞速度的继续增加，尾迹距尾梁的距离增加，因此对尾梁的影响减小。机身与发动舱连接处传感器19的压力幅值随前飞速度增加而减小，原因为传感器在前行桨叶下方，前飞速度增加，为满足配平条件，需减小前行桨叶桨距，前行桨叶桨根载荷减小，因此旋翼桨叶通过性影响减小。
2．4旋翼与机身距离对旋翼／机身干扰影响以ROBIN旋翼在前进比为0．15状态下为基本算例，旋翼与机身高度分别增加5％、10％、15％、20％、30％后机身各部分压力幅值的变化如下。从图11中可以看出，随着旋翼与机身距离的增加，桨叶通过性对机身非定常压力影响减弱，由此导致机身头部、尾梁顶部、机身左右侧各处压力幅值均减小。随着旋翼与机身距离的增加，压力幅值减小速率逐渐减小，距离增加20％，传感器幅值减小为参考值的80％以下。但发动机舱顶部传感器幅值随旋翼与机身距离的增加而先减小后增加，主要原因为旋翼与机身距离的增加，桨叶通过性影响减小，因此压力幅值先减小；旋翼与机身距离的继续增加，方位角为270°处桨叶的桨尖涡将贴近发动机舱顶部，由此导致压力幅值增加。
3结论
关键词: 数值分析 数学建模 Matlab
数值分析又称计算方法,是一门与计算机使用密切结合的实用性很强的一门课程,重点研究如何运用数值计算方法去处理实际工程问题,因此数值分析在科学研究、工程建设和经济建设等很多方面有着广泛的应用。在信息科学和计算机技术飞速发展的今天,这门课程中的数值方法更显得极其重要,但是对多数学校来说,还没有引起对这门课足够的重视,而且在数值分析的教学过程中都存在很多不足。不少学者也讨论过我国高校中数值分析课程的教学情况,其中存在一些普遍问题,例如学生理论学习模式化、实践能力不够、缺乏应用性,学习过程中学生感觉到枯燥或者学习效果不佳,学校软、硬件设施无法满足学生的上机实习等。如何更好地开展这门课程的教学工作,对于我们来说是一个巨大的挑战。下面我们来谈谈在教学过程中遇到的几个问题。
1.理论基础知识扎实,同时采用启发式教学
课程中的很多公式是推导出来的,推导过程比较烦琐,得到的公式也比较冗长,而且比较难记,对于已经复杂并且很冗长的数值公式,还需要进一步进行抽象的理论分析,包括算法的收敛性如何,数值算法是否稳定并进行误差分析,以及分析算法的空间和时间复杂性等,同时还涉及如微积分、线性代数、常微分方程等。过多地强调数学理论证明,大多数的学生觉得这门课很难,学得很枯燥,也感觉不到乐趣,从而越来越厌烦学习这门课程。
因此,我们要将“因材施教”的理念落到实处。方法的讲授应该尽量地从实例中提出问题,引导学生去思考如何运用数学知识去构造解决的方法,然后给出相应的数学理论。并且,给出一种方法,可以换位思考,激发学生思考是否能用另外的已学方法来求解。这样不仅能复习已学的知识,而且能巩固各种知识之间的联系,还可以启发学生把学过的知识学以致用,真正了解学习带来的乐趣。
2.将数学建模的思想融入到教学过程中
数值分析是对实际问题的数值模拟方法的设计、分析与软件实现的理论基础。要解决具体的实际问题,首先需要建立起适当的数学模型,将实际问题的解决归结为相应的数学问题的求解,然后对所归结的数学问题建立相应的数值方法。这样就可以以实例启发学生弄清为什么要进行数值分析、应该如何引进数值方法进行分析,建立一种数值分析的方法后,哪些问题是值得且必须研究的。例如在汽车、飞机等的外形设计过程中,利用样条技术设计的外形越来越光滑、美观。学生了解了样条插值的实际应用背景后就会对样条插值的理论更感兴趣,也会更有动力来学。
将数学建模的思想融入到数值分析教学过程中,要求我们必须有一个合适的切入点,不能用数学建模课的内容过多占有数值分析课的教学,因此精选只涉及相应数值分析理论和方法而又能体现数学建模思想的内容,既能吸引学生又是学生以后可能碰到的案例,将其融入到数值分析课程中是十分重要的。下面具体举两个例子,插值方法可以引入人口增长的模型和设计公路平面曲线的问题,常微分方程的差分方法可以引入导弹追踪和估计水塔的流量问题,方程求根的迭代法可以引入一般战争模型,线性方程组的解法可以引入投入产出模型和小行星轨道问题等。
3.结合Matlab进行实践教学
在结合多媒体教学的过程中,尽量地在讲解数学模型的过程中,无论是问题的引入还是算法的讲解和实现,以及结果尽可能地转化成图形等一些可视的结果展示给学生,以激发学生的学习兴趣,引人入胜,Matlab软件的可视化功能能够实现这一点。
在计算机技术飞速发达的今天,只要有效地把教学过程和相关的计算机技术结合起来,就能够做到减轻教师教和学生学的负担,优化学习环境,实现高效教学。在一些数值分析教材中一些常用的算法都已经有了现成的程序,因此在授课的过程中,对这些算法进行展示时,要让学生从中学会如何将一个算法转变成一段程序。鼓励学生自己根据算法写出程序流程图,然后使用Matlab语言将其转变成程序,将自己所得程序与课本中的结果进行比较分析,这个过程有助于学生更好地理解算法,增强学生动手实践的自信心。
4.结语
数值分析是研究数学模型的数值计算方法。随着电子计算机的迅速发展、普及,以及新型数值软件的不断开发,数值分析的理论和方法无论是在高科技领域还是在传统学科领域,其作用和影响都越来越大,实际上它已成为科学工作者和工程技术人员必备的知识和工具。
对于理工科的本科学生而言,它的理论和实践知识对学生的要求都比较高。因此要让学生学好这门课程,需要在教学中采用一些技巧性的教学方法,比如采用启发式的教学方法,融入数学建模的思想,以及结合Matlab进行实践教学等。这样可以调动学生主动学习的积极性,提高学生的综合素质,使学生真正学好这门课程。
参考文献:
[1]赵景军,吴勃英.关于数值分析教学的几点探讨[J].大学数学,2005,21(3):28-30.
[2]孙亮.数值分析方法课程的特点与思想[J].工科数学,2002,18(1):84-86.
Abstract： Based on the bubble dynamic equation under the consideration of liquid surface tension, viscosity and radiative resistance, this essay adopted numerical simulations to investigate single cavitation bubble dynamics with different kinds of acoustic driving frequencies and bubble initial radiuses.
关键词： 超声空化；空化气泡；数值分析
Key words： ultrasonic cavitation；cavitation bubble；numerical analysis
中图分类号：G30文献标识码：A文章编号：1006-4311（2011）01-0196-02
0引言
随着科学技术的发展，超声已在众多领域得到了广泛的应用，在这些应用中，超声空化是引发各种物理、化学和生物效应的主要机理，这些效应与瞬态空化气泡崩溃时所产生的高温高压等现象有关。在研究单一空化气泡动力学过程的方法中，数值分析是除理论和实验方法之外的一种研究方法，至少有两方面原因表明它是必要的。首先，由于气泡运动过程中高度的非线性，使得从理论上建立能够精确描述空化过程的方程实际上是不可能的，其次，微米级大小的空化气泡半径和持续时间为微秒至纳秒级的气泡运动周期使得实验测量也难以进行。本文基于考虑了液体表面张力、液体粘滞性和辐射阻尼的气泡运动方程，采用数值分析中Runge-Kutta方法研究在不同声场频率和气泡初始半径条件下单一空化气泡的运动过程。
1气泡动态的数值分析
考虑了液体表面张力、液体粘滞性和辐射阻尼的单一空化气泡运动方程：
R+=P+-P-P-
-+P+-P（1）
方程（1）是二阶常微分方程，使用Runge-Kutta方法求解时应当用替换法化为形如y′=f（x，y）的一阶微分方程组，再加以使用。方程（1）可化为下列一阶方程组：
y=R′y′=f（t，R，R′）=-（R′）+P+-P-P--R′+P+-PR?佐=R，y?佐=R′?佐=0
若设时间步长为h，则使用Runge-Kutta方法求每个离散时间点上对应气泡半径大小Rn的递推公式为：
R=R+hR+k+k+kR=R+k+2k+2k+kk=hft，R，Rk=hft+，R+R，R+k=hft+，R+R+k，R+k=hft+h，R+hR+k，R+k
设声场激励波形为P=Psin（ωt），若声压P
f=P+-（2）
同时假设气泡运动过程为等温过程，即方程（1）中的泡内气体多方指数n=1。计算时与液体相关的各参数取值分别为：液体密度ρ=1000kg/m3，液体表面张力系数σ=0.076N/m，液体粘滞系数μ=0.001kg/（m•s），液体中声速c=1481m/s，液体中的静压P=1.013×105Pa。
下面研究单频声场激励下，声场频率f和气泡初始半径R对气泡动态的影响。具体为计算f和R在不同的取值条件下，空化气泡半径随时间的变化曲线，即R（t）曲线。
1.1 声场频率对气泡动态的影响给定气泡初始半径R，通过改变声场频率f的大小，讨论fa的变化对气泡动态的影响。设声场激励为P=-Psin（2πft），P=5.0×105Pa，R=0.6μm，由公式（2）得到其自然共振频率f=7.55MHz。若取f=20MHz，则反映气泡动态的R（t）曲线为图1所示，气泡在多个声波周期内做复杂振荡。
若取f=7MHz，则R（t）曲线为图2所示，气泡在一个声周期内即趋向崩溃。需要指出，这里所说的趋向崩溃是指气泡半径在增大到最大值后急剧向R=0趋近。
气泡在不同声场频率激励下其趋向崩溃的程度是不同的。图3显示了声场频率从4MHz变化到10MHz时，气泡在趋向崩溃时的半径与其初始半径之比R/R0的变化趋势。
通常认为:当超声波频率与气泡的自然共振频率相等时，超声波与气泡之间就达到了最有效的能量耦合，气泡将迅速崩溃。但数值计算的结果表明，当f=f时，气泡半径在通常所指的崩溃阶段趋向0，但不为0；当f小于f至一定限度时，气泡半径才在10-5数量级上趋向0，这时可以认为气泡彻底崩溃；当f大于f至一定限度时，气泡可在多个声波周期内稳定振荡，且振荡波形复杂无规律。
同时大量实验研究表明，随着频率升高，声空化过程变得难以发生。对这种现象的定性解释为：频率增高，声波膨胀相的时间相应变短（如f=20kHz，其膨胀时间为25μs；如f=20MHz，其膨胀时间为25ns），气泡核来不及增长到可产生效应的空化气泡，或者即便空化气泡可以形成，但由于压缩相时间也短，空化气泡可能来不及收缩至发生崩溃。为使在较高频率下产生空化，可以提高声强，即空化阈值将随频率升高而增大。此外，从声波的传播特性可知，频率升高，声波的传播衰减将增大，这也使得空化强度减弱以及可能发生空化的区域减小。
1.2 气泡初始半径对气泡动态的影响给定声场频率f，计算气泡取不同初始半径R时的动态曲线R（t）。设P=P=1.013×105Pa，f=20KHz，这是超声工业清洗及声化学中较常使用的频率，按照公式（2）与其对应的共振气泡半径为R=150μm。图4为R从60μm变化到160μm气泡趋向崩溃时R/R的变化趋势。同样可以看出，当声场频率一定时，在该频率上自然共振的空化气泡并没有彻底崩溃，而是那些半径小于自然共振半径至一定限度的气泡才趋向彻底崩溃，半径大于该自然共振半径的气泡将持续振荡若干周期而不崩溃。
一般情况下，当时，计算得出下述结论：对初始半径大于共振半径的气泡，将发生复杂的持续振荡，一般不会崩溃；对初始半径小于共振半径的气泡，随着声压负压相的到来而不断增大，当声压正压相到来时，气泡先因惯性继续生长到最大半径，然后迅速收缩，直到崩溃。
2结论
本文采用Runge-Kutta数值分析方法研究在不同声场频率和气泡初始半径条件下单一空化气泡的运动过程，数值分析结果表明，当给定气泡初始半径大小时，声场频率在小于气泡自然共振频率至一定限度时，气泡将迅速崩溃，而大于该共振频率时，气泡将持续振荡而不崩溃，即随着声场频率升高，声空化将难以发生；当给定声场频率时，只有其半径小于与该频率对应的气泡自然共振半径至一定限度的气泡才会彻底崩溃，半径大于该自然共振半径的气泡将做持续振荡。
参考文献：
[1]冯若.超声手册[M].南京：南京大学出版社，2001.
[2]钱梦J，程茜，葛曹燕.单泡声致发光中气泡的动力学特征――振子模型[J].声学学报，2002，27（4）：289-294．
[3]李信真，车刚明，欧阳洁，封建湖.计算方法[M].西安：西北工业大学出版社，2000.
关键词地基承载力，塑性力学上限，最优化方法。
1前 言
地基承载力、土压力和边坡稳定是土力学的3个重要领域。这 3个问题都基于共同的极限平衡分析原理 ,可以采用相同的分析方法。但是 ,在长期的实践中 ,这 3个领域各形成了自己的体系。在地基承载力领域 ,目前常用的计算方法仍然是基于 Prandtl解的各种经验修正公式。应用塑性力学上下限原理 ,在建立地基承载力、土压力和边坡稳定分析统一的理论和方法方面作了大量的工作 ,但其有关的研究一直是在变分原理基础上进行的 ,因此 ,难以扩展到具有复杂边界和分层土体的实际工程问题中。曾提出一个基于滑楔破坏模式的分析方法 ,其普遍适用性还有待进一步论证。显然 ,只有开发数值分析的方法 ,方可使大部分实际问题方便地获得解答。
近期 ,笔者在二维领域应用塑性力学上限定理进行边坡稳定的理论研究[4 ,5]。该方法从变形协调出发 ,对于一个设定的滑裂面和斜分条模式 ,建立协调的速度场 ,根据外力功和内能耗散相平衡的原理确定相应的安全系数或加载系数 ,然后应用最优化方法 ,确定对应于最小安全系数的那个临界滑裂面和斜分条模式 (以下简称能量法 )。这一方法在精确地确定边坡稳定安全系数方面获得了成功。由于地基实际上是一个坡度为零的边坡 ,将该成果推广到地基承载力 ,是一个十分具有理论和实用价值的课题。
2极限分析法的理论基础和计算步骤
2.1上限定理的基本命题
在边坡稳定和地基承载力分析领域 ,对上限定理的描述可以用下面的命题表达(图 1):
在塑性区 Ω*,给出一个机动可能的应变场εij* ,并在滑裂面 Γ* 上给出一个相应的速度场 V*,那么,按照下式计算获得的外荷载T* 将比一个包含有真实的塑性区 Ω和真实的滑裂面 Γ的临界荷载 T大或与其相等。
∫Ωσi*jεi*j dΩ +∫vdDs* = WV* + T* V*(1)
上式左边的第一、第二项分别为塑性区内和滑裂面的内能耗散；W为塑性区土体重。因此 ,在诸多协调的位移场中给出最小的 T* 的那个一定最接近真实的临界荷载 T。
在地基承载力问题中 ,通常定义加载系数 η* 为
η* = (T * - T0)/ T0 (2)
其中T0为地基的实际承受的外荷载 ,那么上限图定理的命题具体化为寻找一个使 η* 获得最小值 η的应变场和速度场。
2．2计算内能耗散
如果材料遵守莫尔-库仑破坏准则和相关联的流动法则，则可确认沿滑面的速度V与滑面夹角为摩擦角φ，单位面积内能耗散为（图2）：
d D = (c cosφ u sinφ)v(3)
其中c为凝聚力；u为孔隙压力；V为滑块沿滑面的在单位荷载增量下产生的相对位移 ,通常称变形速率。
2.3计算多块体破坏模式协调的速度场
对某一边坡的塑性区 ,将其用一系列倾斜的线分成若干楔块 ,每一楔块都视为刚体 ,其变形速率为 V。图 3示出 3个块体的系统。 V与滑面夹角为 φ,与右边相邻块体的相对速度为 Vj ,V j与该两块体交界面的夹角为 φj。内能耗散发生于该楔块的底面和楔块间的界面 ,在刚体内为零。
根据位移协调要求 ,可以得到
1 (4)
1 (5)
其中Vl 和 Vr 分别为左侧和右侧条块的速度 ;θj =π/2-δ+ φj, θl=π+αl-φl,θr =π+αr -φr;α为底面与 x轴正向夹角 ;δ为侧面与 y轴正向夹角 ;θ为速度与正 x轴的夹角。如果将条块的宽度取为无限小 (图 4),还可通过积分获得滑面上坐标为 x的条块的绝对速度和相对速度。
V0为左端点（x=x0）的速度。在滑裂面上第k个α或φ 发生突变。上标l和代表该不连续点左和右的物理量。计算从第一个界面开始，到第 n-1个界面终止。这样，对滑面上横坐标为x的任意一点，其条块绝对速度V和条块侧面的相对速度Vj 都可表达为滑面左端点x=x0处的速度V0的函数。将获得各条块的绝对速度和相对速度代入式 (3)再代入式 (1),其中式 (1)左侧第一项可通过将Vj替代式（3）中 V获得。消去左右侧 V0 ,就可求解加载系数 η*如下。
定义
其中d W =土条重量 ;T0x, T0y分别为 T0在 x和 y轴的分量 ;L为土条侧面长度 ;η′为水平地震力系数。式 (10)最后一项计及滑面上 (n -1)个 α或φ的不连续点相应的界面上的内能耗散。由式 (2)定义的加载系数可通过下式计算 :
η* = G/Gb(12)
2.4求解临界滑动模式
陈祖煜和邵长明曾详细介绍了对传统的极限平衡法计算最小安全系数和临界滑裂面的数值分析方法。最优化方法为使用计算机搜索临界滑动模式创造了条件 ,这些研究成果可以方便地推广到本文介绍的极限分析方法中。所不同的是 ,滑动模式和垂直条分法相比 ,增加了一个土条界面倾角 δ。每一条块的 δ也将成为自由度。最优化方法将最终找到相应最小加载系数的滑裂面和斜分条模式。具体计算步骤通过下节 [例 2 ]介绍。
3验证
为了验证上述推导的正确性 ,下面通过两个例题进行分析探讨。
[例 1 ]对具有垂直表面荷载的例题 (图 5),索科洛夫斯基 ( Sokolovski , 1954)给出的临界垂直荷载 q的计算公式如下 :
的 q=c ctgφexp-2x)tgφ-1
其中 χ为边坡斜面相对水平线的夹角。相应的临界滑裂面由三段组成 ,AB , CD为直线 ,分别与边坡线和坡顶线夹角为 μ。
μ=
BC为一对数螺旋线 ,其左右边界线 BO和 CO分别与边坡线和 y轴线夹角为 μ。
当边坡处于极限状态时 ,加载系数 η=0。对 AB ,BC和 CD段分别进行积分 ,按式 (13)确
定的 q将使按式 (10)确定的 G为零。
这一实例说明 ,本文提出的上限定理的命题可通过解析解获得印证。
[例 2 ]某一坡角为 35°的均质边坡 ,其水平顶面上作用一均布荷载 ,荷载方向相对铅直
线夹角为 δ′(图 6)。根据索科洛夫斯基 (1954)提出的滑移线方法 ,此题的理论破坏面由直线 AB , CD和对数螺旋线 BC组成 , CD和 CO分别相对铅直线夹角 μ+ρ和 -(μ+ρ)。其中 μ = 45°-φ/ 2 ;ρ为大主应力相对铅直线的夹角。
主要参数 :c = 720kPa ,φ = 37°,χ = 35°,δ′= 24°,理论解提供的解答是 q =6 . 228 MPa , ρ= 28. 4°。理论的滑裂面和土条侧面示于图 6线 4。滑裂面通过联结 4个点的样条函数形成。对设定的初始滑裂面 1和相应的斜分条模式使用式 (12)求得 η3 =0127。从这个滑裂面开始 ,进行最优化方法计算最终得临界滑裂面和条间界面 (滑裂面 3 ,虚线 ),相应 η=01019。滑裂面 2是优化计算过程中通过随机搜索获得的滑裂面。如果用 5个点来模拟该滑裂面 ,则可得到 η=010028。与理论解相比 ,无论是最小加载系数 ,还是临界滑裂面和临界条间界面均十分接近。
通过 [例 2 ]说明 ,应用最优化方法可以自动找到相应最小加载系数 η的临界滑裂面和相应的斜分条模式。
4能量法在地基极限承载力计算中的推广
4．1传统的承载力计算方法
地基极限承载力的计算包括两个方面，一方面是允许位移的校核 ,另一方面是极限承载力的计算。对于后者， Prandtl于 1920年根据塑性力学理论导出了刚性基础压入无重量土中滑裂面的形状及其相应的极限承载力计算公式。由于数学上的严格解答在大部分的实际问题中是不可能得到的 , Terzaghi ,Meyerhof ,Vesic等众多学者在 Prandtl解的基础上对承载力理论进行了研究和发展 ,最终形成地基极限承载力的近似解答。这一近似解答的一般表达式为
Nγ为一半经验数据 ,可从地基规范承载力表中查取或用半经验公式 (表 1)计算 ;B为基础宽度 ;D为基础埋深 ;γ为土容量 ; qu为地基的极限承载力 ,即 T在单位宽度上的强度。
4.2能量法在地基极限承载力计算中的应用
选取宽度 B = 17m的条形基础进行计算分析 ,相应参数为 :c = 144. 5kPa ,γ=0. 0kN/ m3。这个例子针对土的不同内摩擦角 φ值进行计算。对于具有理论解的实例 ,使用式 (13)获得的 qu应保证使用式 (12)获得的 η的最小值为零。图 7示出 φ =0°和φ = 20°两种情况。使用同样的初始破坏模式如图 7 (a),应用最优化方法获得的临界破坏模式分别如图 7 ( b)和 7 (c)所示 ,η分别为 01004和 01008。计算机在搜索最小 η值时 ,准确地将中部的土条侧面收敛于地基的左侧点 ,由此将滑裂面分为 3个区域 :荷载作用下面的三角形区域 ,对应于理论上的主动 Rankine区 ;条间界面一端收敛于一点所对应的放射形区域 ,对应于理论上的 Prandtl区;放射形区域另一端的另一个三角形区域对应理论上的被动 Rankine区。 φ=0°时 ,滑裂面形状接近于圆弧 ;φ≠0°时 ,滑裂面形状为两段直线接一段对数螺旋线。这就说明 ,采用建立在塑性力学上限 解基础上的地基极限承载力数值分析方法直接获得了理论上严格的地基极限承载力解答。
表 2将一系列φ值的计算成果与理论解对比，可见成果的准确性相当稳定。所得的η值与理论值的误差均在1%左右，而且自动搜索得到的临界滑裂面形状也与理论解一致。
表2 上限数值解和理论解成果对比
注：qn和 qp分别为根据数值解和理论解获得的极限承载力 ;η= ( qn/ qp -1);q的单位为 kN/ m。
表 3比较了容重不为零情况下表 1所列的各种经验方法的准确性 ,并示于图 8。可见 W1F Chen的公式计算结果与采用 Prandtl的经验公式求得的结果最为接近 ,但是在 φ值超过 25°后 ,η值为负 ,意味着高估了地基的承载力。而 Mayerhof以及 Terzaghi的方法则偏于保守。
表 3对 γ≠0情况各种不同的经验公式和数值解成果对比
注 :下标 n,v, m,t,c分别代表数值解、采用VesicMeyerhofTerzaghi 和 W. F. Chen方法的计算成果，q的单位为 kN/ m。，
在有容重且有埋深的条件下 ,式 (15)的经验公式将基础两侧埋置深度以内的土重以连续均匀分布的荷载考虑 ,未能计及这部分土体的抗剪作用 ,因此不可避免地存在着一定的误差。
[关键词]压缩式封隔器；扩张式封隔器；封隔器胶筒；有限元分析
1.封隔器的应用现状
压缩式封隔器的工作原理：当胶筒承受轴向载荷时，封隔器胶筒将产生径向大变形，使胶筒与套管之间产生接触压力，藉此封隔环空，隔绝产层，以控制产（注）液，保护套管。胶筒所用材料为橡胶，橡胶最显著的特性就是它的超弹性，即在很小的力的作用下就能产生很大的变形，正是这一性质使得橡胶材料被制成密封元件并在工程中得到广泛应用。
但是，国内还没有专门用于深水井下作业的相关封隔工具，因此，研究封隔器胶筒的材料以及研究用于恶劣环境下的井下措施工具以及技术，是国内各大研究单位以及相关高校的当务之急。
2.胶筒材料本构关系
2.1 基本理论
橡胶材料在较短时间内及恒定的环境温度下通常被视为各向同性不可压缩材料，其应变能密度函数W是变形张量不变量的函数，即：
（1）
其中：
（2）
（3）
（4）
式中：λ1、λ2和λ3是3个主伸长率。
Rivlin用唯象法从纯数学角度出发导出应变能密度函数的最一般的形式为：
（5）
式中：Cijk为Rivlin材料常数。
橡胶材料可以近似认为是不可压缩的，因此I3=1，式（5）可简化为：
（6）
取式（6）的前两项得到Mooney-Rivlin模型：
（7）
其中C10和C01为Rivlin材料常数，均为正定常数。
2.2 Mooney-Rivlin模型材料常数的确定
对两参数Mooney-Rivlin模型，Mooney-Rivlin材料模型的应力-应变方程可表示为：
（8）
这个方程可用σ/2（λ-λ-2）对1/λ作图。在1/λ=1时，相应值为C10+C01并且直线的斜率为C01，初始剪切模量与材料常数的关系是G=2（C10+C01）。如果材料被假定为不可压缩的，那么初始拉伸模量是E=6（C10+C01），G=2（C10+C01）和杨氏模量的值近似等于3倍的剪切模量。
利用橡胶材料的硬度与C10和C01存在一定的关系，根据橡胶材料硬度与弹性模量的实验数据，经拟合可得：
（9）
利用硬度计可测得橡胶试样的硬度HA，通过式（9）可算出橡胶弹性模量E，根据公式E=6（C10+C01）可以得到C10+C01的值，通过公式d=（1-2μ）/（C10+C01）可以测出模型参数d的值。利用经验公式C01=0.5C10求解下面的方程就可得到C10。
E=6（C10+C01），因此E=6（C10+0.5C10），则。
由实验测得在常温25℃和高温150℃下胶筒的邵氏A硬度分别为94和91，通过公式（8）、（9）即可计算出的胶筒的Mooney-Rivlin模型材料参数。
3.封隔器胶筒有限元数值模型的建立
封隔器采用两胶筒形式，上、下胶筒、中心管和套管结构。上下胶筒力学参数和几何尺寸完全相同，两胶筒由隔环隔开。由于封隔器中心管、胶筒、套管以及胶筒所受载荷均为轴对称分布，故取过轴线的剖面建立有限元计算模型， 胶筒采用四节点四边形超弹单元PLANE182，中心管、套管和隔环采用四节点四边形等参数单元PLANE42，用面面接触单元CONTA171和TARGE169模拟胶筒与套管、中心管和刚性隔环间的接触，并考虑接触时的摩擦作用，通过实验测定胶筒与钢材之间的摩擦因数为0.3，钢与钢之间的摩擦因数为0.1，中心管和套管上下两端采用固定约束，橡胶筒在上端隔环的作用下压缩膨胀，与套管接触，随着坐封压力不断增加，胶筒与套管的接触压力也随之增大，从而起到密封套管内上下环空压差的作用。
建立胶筒有限元模型，对胶筒进行非线性接触分析。胶筒具有不可压缩性，因此取胶筒泊松比为0.4996。
4.数值模拟结果分析
4.1 坐封过程中接触压力分布规律
为模拟封隔器坐封过程中胶筒与套管内壁接触情况，稳定温度150℃不变，分步向胶筒上端隔环施加轴向荷载，观察胶筒的变形情况。
由于胶筒与套管内壁接触压力和摩擦阻力的存在，2个胶筒所受到的轴向压缩载荷由上至下依次降低，所以尽管上下胶筒的材料和几何形状完全相同，但是下胶筒的变形、与套管的接触压力均小于上胶筒，坐封完成以后，上胶筒的接触压力比较大，下胶筒的接触压力比较小。显然，随着封隔器坐封载荷的增加，上胶筒与套管的接触压力明显高于下胶筒的，所以上胶筒在密封套管环空上下压差中起到关键作用。
4.2 坐封过程中轴向载荷与胶筒变形的分布规律
根据ANSYS分析输出结果可以得到在150℃时不同坐封力下胶筒的轴向压缩位移和轴向载荷的关系。稳定温度在150℃时，依次施加10t、12t、13t、15t、17t、18t坐封力时胶筒的轴向压缩位移分别为23.095mm、24.319mm、24.564mm、24.905mm、25.086mm、25.186mm。
实验发现，随着轴向载荷的逐步增大，胶筒轴向压缩位移也渐渐增大，当轴向载荷小于12t时，胶筒轴向压缩位移随轴向载荷增大的梯度比较大，当轴向载荷大于12t时，胶筒轴向压缩位移增大梯度慢慢减小。因此，当轴向载荷过大时，胶筒轴向压缩位移的变化不是很大。
5.结论
从数值模拟可以看出，封隔器的两个胶筒与套管的接触压力分布为：下胶筒接触压力较小，上胶筒的接触压力较大，随着轴向载荷的增大，轴向压缩量也增大，同时也可以发现，随着坐封力的增大，胶筒与套管接触长度逐渐增加，胶筒外表柱面部份径向变形受限制，胶筒内表面变形如外表一样向外鼓，当载荷增加时胶筒被压扁并在最后压实。并且随着轴向载荷的逐步增大，胶筒轴向压缩位移一开始是渐渐增大，增大的梯度也比较大。当轴向载荷继续增加时，胶筒轴向压缩位移增大梯度会慢慢减小。因此，当轴向载荷过大时，胶筒轴向压缩位移的变化不是很大，为了达到良好的密封，应该计算合适的坐封载荷。
参考文献
[1]曾宪平.对压缩式封隔器胶筒变形分析[J].石油钻采工艺，1983，13（4）：65-68.
[2]杨秀娟，杨恒林.液压封隔器胶筒座封过程数值分析[J].石油大学学报，2003，05-0084-04.
[3]刘永辉.封隔器胶筒性能有限元分析及结构改进[D]西南石油大学，2006
关键词：混凝土 ;随机骨料模型; 映射网格;细观数值模拟
中图分类号：TU528.01
文献标识码：B
文章编号：1008-0422（2008）08-0163-03
1前言
混凝土作为重要的建筑材料已有百余年历史，人们从不同的结构层次来研究混凝土材料的力学性能。（1）宏观层次。一般认为特征体积相当于3～4 倍的最大骨料体积。当小于特征体积时，当研究对象大于特征体积时，材料假定为均质。（2）细观层次。从分子尺度到宏观尺度之间，其结构单元尺度变化范围在10-4cm至几cm。在这个层次上，混凝土被认为是一种由骨料、水泥砂浆及二者之间的界面组成的三相复合材料[1]。(3) 微观层次。材料的结构单元尺度在原子、分子量级，即从小于10-7cm～10-4cm 。材料由晶体结构及分子结构组成。随着高性能计算机及数值算法的迅速发展，用数值试验来模拟混凝土材料破坏过程的研究已成为热点，相对廉价的细观数值分析有望弥补或替代部分材料试验[2]
2骨料尺寸
选取300mm×300mm×300mm的立方体试件，并将它们取为二级配混凝土，其中大石(直径40mm~80mm)、中石(直径20mm~40mm)进行数值模拟。按富勒理想级配曲线：
（1）
其中D为筛孔直径；Dmax为骨料最大粒径；P为能通过筛孔D的骨料的质量百分比。
取Dmax为80mm，求得大石、中石的质量百分比：
由表1可知，粒径为40mm~80mm的骨料质量百分比为100%-70.7%=29.3%，粒径为20mm~40mm的骨料质量百分比为70.7%-50%=20.7%。骨料容重取2800Kg/m3，试件的石料用量为1227Kg/m3。适当选择球体骨料的个数与直径，使之满足级配关系并容易随机填充，本文取大骨料为大石，个数为27,直径为62.58mm。小骨料为中石，个数为189,直径为29.16mm。
3三维有限元网格的划分
随机骨料模型的网格剖分，长期以来是三维细观力学分析的难题之一。分网技术简单地说分为三种：六面体映射或拉伸网分、四面体全自动网分、混合网分。（1）商用软件中，TRUEGRID、ANSYS ICEM CFD、AVL_FAME等可以完成六面体映射网格划分。划分随机骨料模型时，主要思想是：先生成均匀的有限元网格覆盖整个分析区域，根据蒙特卡罗方法和骨料从大到小的投放原则，随机生成一套虚拟的骨料域。根据单元节点是否全部落在、部分落在和不落在骨料域，判别单元是否属于骨料单元、界面单元和砂浆单元。这种映射网格技术能够解决复杂形状骨料的三维模拟，全部是六面体网格，有限元的单元精度高，但界面区域不够光滑。在计算上，考虑到界面的厚度很小，要采用较小的单元而在全域上采用一样大小的单元，导致单元个数太多，计算量较大，在微机上计算难度大。（2）多数的CAE软件都提供四面体网格全自动划分功能，配合局部加密技术划分单元，可以建立随机骨料网格模型，计算量可以接受，但四面体单元的精度较低，空间填充效率不高，且视觉效果稍差。（3）CUBIT、MARC、HARPOON等软件提供混合分网功能，产生的网格主要为六面体，边界复杂的区域用五面体、棱柱体或四面体填充。分网效率高，但边界区域通常有应力集中，需要特别关注的地方却用的是低精度网格。
本文提出先分块映射网分再随机装配的方案，可以兼顾单元精度、局部加密和计算量的要求。首先，将混凝土试件看成由三种块体组合而成：（1）包含大骨料、界面和砂浆的块体，简称大骨料块。如图1所示。（2）包含小骨料、界面和砂浆的块体，简称小骨料块。如图2所示。（3）只包含砂浆的块体，简称砂浆块。如图3所示。分别对这三种块体用映射网分。
以上是球体骨料的情形，类似地，可作出椭球体骨料的映射网格。如图4所示。
注意到界面的材料刚度与骨料的刚度有较大的差别，为了减少应力集中，将骨料的网格作了适当的加密。
三种块体都是正方体，使大骨料块的边长为66.7mm，小骨料块和砂浆块边长33.3mm，界面单元厚度为1mm，骨料和界面单元全部位于块体的内部。大骨料块的边长为小骨料块和砂浆块边长的两倍，表面上的单元尺寸相同，显然这三种块体可以无缝、无重叠地填充三维空间，而且将重合的节点合并后，各块之间形成共节点连接的整体有限元网格，如图5所示。随机装配过程利用了MARC软件的PYTHON语言前处理接口，调用随机函数以及单元的复制、移动功能。整体有限元共有252424个六面体单元，258877个节点，其中砂浆单元139968个，如图6所示。骨料单元82922个，如图7所示。界面单元28512个，如图8所示。
4材料参数
Drucker -Prager 屈服准则在MSC.Marc有限元软件中称其为Linear Mohr-Coulomb准则，屈服条件为[4]：
为应力张量的第一不变量； 为偏应力张量的第二不变量。
MSC.Marc提供了可模拟开裂和应变软化的损伤模型。各主要参数如表2所示。
5结果分析
按位移控制加载，铅垂方向分级加拉伸位移，ε表示混凝土试件在加载方向的总应变，图10~图14显示的是等效开裂应变云图，相当于显示出了裂纹分布的区域。
图10显示当伸长总应变为ε=0.663×10-4时在大骨料的界面上开始出现裂纹，裂纹方向为垂直于加载方向。图11显示当伸长总应变为ε=0.773×10-4时在小骨料的界面上也开始出现裂纹，而且大骨料界面上的裂纹继续增加，裂纹方向为垂直于加载方向。图12显示当伸长总应变为ε=0.8673×10-4时，大骨料界面上的裂纹继续增多，在小骨料的界面上裂纹没有明显增多。
图13显示当伸长总应变为ε=0.970×10-4时，在小骨料的界面上裂纹有明显增多。图14显示当伸长总应变为ε=1.07×10-4时，在大、小骨料的界面上裂纹足够多了。虽然裂纹之间还没有贯穿相连，但材料非线性的求解得不到下一载荷步的收敛解答，可认为材料已达到单向拉伸的极限载荷。
6结语
6.1本文从细观力学角度，使用分块映射网格划分法结合各子块网格随机组合的思路，建立了混凝土三相细观有限元模型，模拟了混凝土试件在单项拉伸时裂纹的产生、发展过程，结果与试验结果定性相符。
6.2没能模拟出应变软化现象，综合实验、理论及数值方法对混凝土应变软化现象的研究还需大量的工作。
参考文献:
[1] 刘光廷，王宗敏. 用随机骨料模型数值模拟混凝土材料的断裂[J].清华大学学报，1996，36(1):84~89.
[2] 马怀发，陈厚群，黎保琨.混凝土细观力学研究进展及评述[J].中国水利水电科学研究院学报，2004，2(2) :124~130.
[3] 彭一江，黎保馄，刘斌.碾压混凝土细观结构力学性能的数值模拟仁[J ].水利学报，2001，6(10):19~22.
关键词： 数值分析 课程教学 教学改革
“数值分析”又称“数值计算方法”，是研究如何利用计算工具（如计算机）求解数学问题数值解的一门科学，因此数值分析随着计算工具的发展而发展。它与传统的数学课程有明显的区别，它与概念、逻辑推理及各类问题求解的数学技巧相对而言不如先前数学课程那么重要，对具体问题求解方法的构造及利用数学理论对各类算法的收敛性和稳定性分析是课程的研究核心。根据这门课程的特点及授课对象的差异，本文提出了课程教学改革的思路。
一、更新观念、改革课程内容
在理、工科数学类科目教学体系中，“数值分析”是一门重要的课程，理工科学生在各自专业学习中面临各种专业问题需要建立数学模型求解，复杂的数学模型往往得不到解析解，那么如何求解？数值分析正好为他们提供了解决思路，因此能够激发学生进一步学习数学、应用数学的意识和兴趣，同时也能培养学生的创新思维和创新能力。
《数值分析》课程教学内容多，公式复杂，推导繁琐，而学时少，教师不可能对每个知识点讲解透彻，且算法实现对编程能力的要求较高，因此不少学生由于“掌握不了”容易厌学[1]。笔者认为根据学生不同专业、不同院校可以选择符合自己特色的教材。该课程授课对象主要有两类：一类是理科学生，他们学习的目标是不仅要会“使用”算法，还要会“创造”算法。另一类是工科学生，他们学习这门课程的目的主要是“使用”算法。由于这两个专业群体学习的目标不同，教材的内容、体系及侧重点应有差别。理科生学习时需要“研究”算法，学得相对精细点，学时不够，因此在学习内容的选取上可以删去一些在后续课程中将会学到的内容，比如常微分方程数值解。针对学生编程基础薄弱，可选用基于MATLAB的数值分析教材[2]，由于MATLAB编程简洁、数据处理方便，具有强大的数值计算功能，即使不具备较强编程能力的学生，学习这种教材也不至于为算法实现而发愁，为提高数值分析课程的教学质量创造了良好的条件。笔者认为基于MATLAB的数值分析教材符合创新性实践教学的要求。当然对于编程能力强的部分理科学生来说，学习C语言等高级语言的教材有利于加深对算法思想的理解。
二、课程教学方法和教学手段的创新
学生是学习的主体，教师是教学活动的引导者、组织者，因此应充分发挥学生的主体作用。传统的教学模式是教师在课堂上讲授，学生被动接受，“满堂灌”的教学模式，教师往往费力不讨好。数值分析的教学目的不仅让学生掌握算法的基本思想、基本方法和基本原理，还通过这门课程的学习提高学生数学素养，注重培养学生举一反三的能力。因此，在教学内容的讲解上要有所侧重，以主带次，在有限的学时里讲清每一个主题，突出讲授典型的、具有代表性并能体现其思想方法的常用算法和理论；对那些原理相近的内容不求面面俱到，可以精心设计教学提纲加以引导和提示，并且开发网络教学平台，把教学资源于网络平台，让学生利用这一平台以小组为单位课后讨论自学完成。数值分析课程教学可以贯彻“少而精”的原则，达到创新教学的目的。
课堂教学应激发学生的兴趣，在教学方式上，力求讲清每一个主题的实际背景、目的和算法设计的出发点，通过对实际应用背景的描述，激发学生的学习欲望。不少学生在前期数学学习中感到困惑：数学到底有什么用？以实际问题为出发点，让学生真切地体会到学习数值分析课程的意义所在。其次，充分调动学生的积极性和主动性，启发学生进行独立思考，带领学生分析算法的适用范围、优缺点，启发学生对算法进行改进。让学生明白任何一种算法都有局限性，有值得改进的地方，在课程学习中激发学生的科研兴趣。针对某一个知识点，要求学生课前查阅资料认真预习，让一些学生像老师备课一样精心准备，并上台讲解。学生经过查阅资料、独立思考这一主动学习过程提高兴趣，感受到成功的喜悦，增强自信，由对这门课程最初的“害怕”逐步变为“喜爱”，也通过这一知识点的学习掌握学习方法。我们在信息与计算科学专业学生中进行这方面的尝试，实践证明通过这种以点带面的学习，大大调动了他们的主动性和积极性。
三、重视实践环节的教学和改革
传统的课程教学重理论、轻实践，而《数值分析》是一门实践性很强的课程，学习这门课程时，应该将重点放在掌握数值分析的基本理论和思想及与结合计算机解决实际问题上，所以必须重视《数值分析》的实践教学。
实践教学环节可以分为两大类：一类是课堂实验，另一类是与工程实际背景相关联的数学问题的求解，比如数学建模竞赛、与学生专业领域相关的问题求解。对于课堂实验，教师应编写实验指导书，鼓励学生分组讨论，让学生通过实验进一步体会算法思想，掌握算法精髓。比如学习用Newton迭代法解非线性方程时，要求初值比较接近真实解，这时引导学生实验时选择不同的初值进行比较，并归纳结论；由于实际问题是不知道真实解，难以满足初值要求，因此对算法进行改进，提出了下山牛顿法，这一过程让学生在实验中自己感受，印象会深刻得多。对于工程实际问题的求解，可以引导学生分析全国大学生数学建模竞赛历届试题，并将数值分析课程中的逼近思想，如插值法、最小二乘法、曲线拟合等方法用于求解中，我们必须将数学建模的思想融入数值分析的教学中。对于学生专业领域中的数学问题可引导学生分组讨论，自己设计算法方案。通过实践环节的教学，让学生真正体会到学有所用，同时激发他们的学习热情和科研兴趣。
参考文献：
关键词：数值分析课程设计；实践教学；教学改革；教学方法
中图分类号：G420 文献标志码：A 文章编号：1002-2589（2013）24-0291-02
《数值分析》是一门研究数值求解各类数学问题的课程，是综合性大学数学专业与某些理工科专业的必修课程，也是面向全校学生开设的通识选修课。在《数值分析》课程设计的实践教学中，学生经常抱怨时间短，内容多且复杂，难学且不会应用。事实上，《数值分析》既有数学课程在理论上的抽象性和严谨性，又有解决实际问题的实用性和实践性，它作为科学计算的基础理论与基本方法，已经广泛地应用于物理学、力学、计算机应用、航空航天、土木工程、机械工程、风险投资和经济管理领域。如何改变学生的观点、提高学生的兴趣、提高《数值分析》课程设计的教学质量，已成为数学专业课程教育教学改革的焦点之一。笔者经过几年的《数值分析》课程的教学，分析了《数值分析》课程设计实践教学中存在的一些问题，并根据学生专业的特点提出几点教学改革建议。
一、实践课程中存在的问题
《数值分析》课程设计实践教学，课程难学与本课程教学过程中的一些特点有关。具体说有如下几点。
（一）任课教师的专业限制
《数值分析》课程是信息与计算科学专业的专业必修课，它的任课教师是数学专业教师。而《数值分析》课程设计培养的是厚基础、宽口径、应用型本科人才，属于应用教育。在教学过程中，任课教师常常强调专业和基础，忽视实践与应用，这使得学生要花许多时间学习，才能学会该课程，如果不这样做，很难达到老师的要求。这样学生的实践练习时间就很少，更难达到用所学知识来解决实际问题，造成教与学脱节、学与应用脱节。另外由于数学专业教师缺乏工科背景，教学中多是只给出概念、公式、定理，罗列出方法和计算，很难做到时时处处与实际应用需要结合，这样会让学生感觉该课程就是一门枯燥的数学课，很难做到学以致用，逐渐地丧失学习兴趣。
（二）学生的专业特点限制
《数值分析》课程是面向数学系信息与计算科学专业开设的专业必修课，而信息与计算科学专业的学生数学基础好，计算机基础薄弱，程序操作能力差，课程设计实践教学的目的是各种数值算法的计算机实现，学生本就觉得这门课程理论枯燥难学，由于不熟练编程序，很难有兴趣自己操作实现各种算法，违背了课程设计的目的，使得这部分内容多是老师讲解演练各种算法，学生很少动手操作，教学效果较差。
（三）数值分析课程设计实践教学时间短，内容多
《数值分析》是一门与计算机密切结合的课程，该课程的任务是要根据计算机的特点给实际问题提供切实可行的算法，因而，课程设计是该课程教学中一个不可缺少的环节。然而，由于学时较少，在保证理论的限制下，课程设计实践教学只有一周时间。而且《数值分析》课程设计多安排在教学的最后一周，学生在学期最后几天往往会因很多事情影响学习，实际授课和实践的时间就更少了，在这样短的时间内完成许多算法的实践是很难做到的，只能选择几个重点的算法由老师讲解学生练习，学生只掌握了几种算法的原理，掌握效果也并不好，很难用所学知识去解决实际问题，这与《数值分析》课程设计解决实际问题的实用性和实践性相违背。
（四）数值分析课程设计教学模式单一
数值分析课程设计教学目标是教学生将所学的理论知识转化为实践能力，采用的教学模式就是上机讲解和练习。连续五天的这种教学模式会让学生觉得形式单调、内容枯燥、难学，在没有老师在身边监督下，有的学生干脆不听不练习，玩起了游戏，老师过来检查时关掉游戏，老师离开又继续玩，反反复复地看着老师，几天的实践课程就这样荒废，很难达到预期的教学效果。
二、教学改革的几点建议
针对上述的几点问题，笔者在《数值分析》的教学中提出以下几点改革建议。
（一）明确教学目的
必须明确我校人才培养的目标不是专业人才，而是应用型本科人才，培养操作能力强、宽口径、复合型的人才。针对这种人才培养模式的特点，我们确立《数值分析》的教学目的是：通过这门课程的学习，使学生了解数值计算的基本原理，掌握一般的数值计算技能，增强学生实践能力，使学生能够在实际需要涉及无法求解的数学问题时，会使用数值计算的思想和方法，为以后的实际操作和研究工作打下良好的基础。
（二）充实和精炼教学内容
《数值分析》是一门严谨、完整的学科体系，内容丰富。连同课程设计部分它的学时是64学时，笔者在近几年的授课过程中感觉到，在少学时的情况下，如果所有内容都讲到，往往会顾此失彼，不能突出重点，学生学下来后并没有在脑海里对所学知识留下深刻印象。如果将实际应用较多的内容重点强调，相关的方法应用详细讲解，并在上机实验上重点练习，其他知识点弱化讲解，应用较少却很繁琐的内容精简掉，这会使得与实际应用相关的重点内容有充足的时间讲解，实用性很强的方法有动手操作的机会，不仅减轻了学生的学习负担，还激发了学生思考用这些方法解决相关问题的兴趣，增强了学生解决实际问题的能力。
（三）不断改善实践课教学模式和教学方法
《数值分析》研究的是各类数学问题的数值算法，对于实用性较强的算法一定要多给学生动手操作的机会，加深学生对这些算法的理解，熟练它们的应用。数值试验的题目选择一些与实际应用相关或生活相关的案例，让学生将分析问题、建立模型、选择合适的算法、编写程序、分析结果这一计算主线完整地进行下来，以缩小学会该课程与解决实际问题的距离，增强学生的能力和信心。数学软件可以指导学生使用简单且方便调试的软件，如Matlab等，这可以兼顾到学生计算机语言基础薄弱，避免调试程序时语法错误频出而丧失信心。在教学方法上，做以下改革：
1.采用在解决问题中引出理论内容的教学法
每节课前都提出一些与学生专业和生活联系紧密的问题，和学生一起探讨这些问题的解决方法，越简单越好，从实际问题的解决方法入手，运用归纳、分析的手段引出相应的数值计算理论，之后从理论上研究解决问题的思想和方法，分析方法的优点、缺点以及所能解决问题的类型，进而给出解决实际问题的程序。例如在介绍数值积分时，课前不妨假设要购买一块不规则形状的养鱼池，按面积付费，判断所付价格是否合理。以其中两点确定的方向为横坐标轴，选择坐标轴上的点后，用船行驶的时间乘以船速，得到不同点的纵坐标，每小块面积用矩形面积近似，这些近似值加起来得到养鱼池面积的近似，将这种思想和方法提炼出来就是数值积分的思想和方法。结合生活中的实例，学生既能掌握《数值分析》中各部分内容，又能了解其在生活中的应用，进而将遇到的问题分类，学会用不同的数值方法解决相应的问题。
2.采用优劣分析式教学方法
这里所涉及的优劣分析包含数值解法与解析解法优劣分析、同一问题不同解法优劣分析等。
如讲到数值积分时，选几个实际生活中遇到的问题为例，如购买性状不规则的养鱼池，需要估计面积，而在边界曲线函数无法知道的情况下，不能直接应用计算定积分的牛顿-莱布尼兹公式来计算面积，而若被积函数是一张数据表、用牛顿-莱布尼兹公式积分后的原函数过于复杂等，都无法直接应用牛顿-莱布尼兹公式计算，而数值积分公式则可以克服这些弊端，强调了数值积分的实际应用价值，给出具体的数值算法，比较数值算法的计算量。又如，在介绍线性方程组的数值解法时，指出线性代数中的线性方程组的求解方法只对阶数较低的方程组适用，阶数稍高就会导致计算量很大，列举一个五阶线性方程组和一个十阶线性方程组，观察它们的计算量，而源自工程实践中的线性方程组多是高阶线性方程组，如土木工程中的有限元方法，最后都会归结为解高阶线性方程组，如果用传统的解法会导致计算量相当大，而给出数值算法，比较数值算法与传统算法的计算量，这样，数值解法的优点可以吸引学生的学习兴趣。
启发学生思考身边问题，发现问题、解决问题。课上，老师讲授之外，师生共同围绕学生提出的那些与生活相关的问题进行讨论，提出解决问题的方法并在程序上实现问题的解决，对大家提出的算法进行点评，这种授课方式能够激发学生学习的热情。
（四）改革考核方法
以往的实践课考核方式是给出几道题目，给学生一定的时间，由他们自己设计算法，给出问题的解决方案，提交课程设计报告，根据这些报告给定成绩。有的同学平时不练习，在最后抄袭别人的实践报告。为了杜绝这类现象，在平时上课时，将学生分组，每次课给每组指定任务，大家分组讨论，每组都给出问题的解决思路，每个组员都要给出补充，根据这些课上表现综合给出考核成绩，这既可以监督学生浪费课上时间玩游戏，又可以保证所给成绩的公平性，调动学生积极性。
三、结束语
教学内容和教学方法的改革，是提高教学质量的关键，是高校提高自身的竞争力适应社会发展需求的必然趋势。讲授《数值分析》课程设计，必须做到与实际问题相结合，根据实际应用取舍教学内容，安排试验学时，培养学生应用所学知识解决实际问题的能力，拓展学生思维，为培养基础全面、操作能力强、综合素质好的应用型本科人才奠定基础。
参考文献：
[1]李庆扬，王能超，易大义.数值分析：第5版[M].北京：清华大学出版社，2008.
[2]周生田，李维国.工科硕士研究生数值分析课程建设与教学改革[J].石油教育，2009，（1）：57-59.
关键词：柴油机；螺旋进气道；设计参数；缸内涡流比；数值分析
中图分类号：TK422 文献标识码：A
Numerical Analysis of the Effects of the Design Parameters of a Diesel Helical Intake Port on the Swirl Ratio
HAN Zhiyu1,2，WANG Yong1，CHEN Zheng1，DENG Peng1
(1.Research Center for Advanced Powertrain Technology，Hunan Univ，Changsha，Hunan 410082，China; 2.State
Key Laboratory of Advanced Design and Manufacturing for Vehicle Body，Hunan Univ，Changsha，Hunan 410082，China)
Abstract:To study the effects of the design parameters of a helical intake port on the incylinder swirl ratio for an automotive diesel engine, the transit port and incylinder flows were simulated with the Converge CFD software. By changing the values of three main design parameters of the intake port, their effects on the incylinder swirl ratio were simulated and studied. The results have shown that the incylinder swirl ratio achieved the maximum values when the spiral chamber height, the vortex shell cutoff amount and the helical intake port angle were changed to the values of 13 mm, 1.27 mm and 0 degree, respectively. The results also indicated that the effects of the parameters on the swirl ratio were not additive when multiple parameters were varied at the same time. This study has shown that the numerical analysis of multiple design parameters of a helical intake port can be done with the use of the Converge software.
Key words:diesel engine；helical intake port；design parameters；incylinder swirl ratio；numerical analysis
随着石油能源的不断减少以及环境问题的日益加重，针对车用柴油发动机的燃烧和排放提出的要求也越来越严格.为了满足日益严格的法规要求，必须合理的组织缸内燃烧过程，以提高燃料的燃烧效率和降低排放.
柴油机螺旋进气道结构参数对缸内涡流的产生有至关重要的作用，合理的气道结构参数不仅能保证在进气过程中产生较强的涡流，还能保证在压缩过程以及压缩上止点附近缸内都具有较强的涡流.通过与喷油时刻的匹配，能保证在喷油燃烧时，缸内具有最优的涡流比，促进燃料与空气的混合，提高缸内空气的利用率，以达到改善发动机动力性、经济性以及排放性的目的［1-2］.本文采用最新的三维CFD计算软件Converge进行数值研究，分析了柴油机螺旋进气道的3个主要结构参数对缸内涡流比的影响，实现了螺旋进气道多结构参数的数值研究.
1 基本情况介绍
1.1 发动机基本参数
建立了一台高速柴油发动机的气道及燃烧三维模型，如图1所示.气道三维模型是通过三维激光坐标扫描仪扫描气道砂芯获取点云图，并导入三维建模软件中得到［3］.该柴油机基本参数如表1所示.
图1 气道及燃烧室三维模型
Fig.1 3D model of transit port and combustion chamber
表1 柴油发动机基本参数
Tab.1 Basic parameters of the diesel engine
型号
排量/L
压缩比
燃烧室
形状
进气
方式
标定功率
/kW
最大转矩
/(N・m)
直列
4缸
1.6
18.1
直喷
SymbolwA@
燃烧室
增压
中冷
80(4 000
r/min)
230(2 000
r/min)
湖南大学学报(自然科学版)2012年
第4期韩志玉等：螺旋进气道设计参数对涡流比影响的数值分析
1.2 三维模拟计算软件介绍
本文采用了最先进的三维计算流体力学软件Converge进行CAE研究.其不同于传统CFD软件之处，在于它的计算网格是在计算过程中自动生成，无需用户提前画出计算网格，从而为用户节省了大量网格制作时间.在使用Converge时，用户只需将制作的STL格式的几何文件导入Converge的前处理软件中进行简单的几何检查和边界设置即可［4］.
2 涡流比与进气流量系数计算
2.1 涡流比与进气流量系数
2.1.1 涡流比的计算
涡流比SR的计算公式为［4-5］：
SR=ΩflowΩcrankshaft. （1）
式中：Ωflow为缸内气体的角速度，r/min；Ωcrankshaft为发动机的曲轴转速，r/min.本文利用Converge计算软件可直接计算出缸内涡流比的大小.
2.1.2 进气流量系数
Ricardo无量纲流量系数C［6-7］F：
CF=QnAVV0,（2）
V0=2・ΔPρ, （3）
AV=π・d2V4. （4）
式中：Q为试验测得的实际空气体积流量，m3/s；n为进气门数目；V0为理论进气速度或速度头，m/s；AV为气门座内截面积，m2；dV为气门座内径，m；ΔP为进气道压差，Pa；ρ为气门座处气体的密度，kg/m3.本文通过Converge软件得出计算的空气体积流量Q，然后代入式（2），得出计算的流量系数CF.
3 边界条件设定以及模型标定
3.1 边界条件的设定
在发动机台架上进行2000转倒拖试验并测得相关试验数据.模拟发动机的工作范围设定为进气阀开启前10°至压缩上止点后20°，则对应的曲轴转角范围为-386°~20°.根据倒拖试验测定的缸压、温度等，设定-386°时刻发动机的边界条件如表2和表3所示.
表2 边界条件1
Tab.2 Boundary condition 1
进气温度
/K
进气压力
/kPa
缸内温度
/K
缸内压力
/kPa
排气温度
/K
排气压力
/kPa
325.8
124.2
352
245.2
342.8
196.5
表3 边界条件2
Tab.3 Boundary condition 2
活塞温
度/K
气缸温
度/K
缸盖温
度/K
进气阀
温度/K
进气道
壁面温
度/K
排气阀
温度/K
排气道
壁面温
度/K
380
353
353
353
343
370
370
3.2 模型的标定
将数值计算的发动机缸压与倒拖试验所得到的实验缸压进行对比，如图2所示.
由图2可以看出：在开始时刻计算缸压和实验缸压存在微小的波动，在之后的时刻，计算缸压和实验缸压十分吻合.在压缩上止点时刻，计算缸压略高于实验缸压，计算的最大缸压（5.812 1 MPa）比实验的最大缸压（5.718 4 MPa）高1.6%.因此，可以认为：数值模型能很好地反映实际发动机的工作情况，可以使用此数值模型进行下一步研究.
曲轴转角/° ATDC图2 计算缸压与试验缸压对比
Fig.2 Comparison between calculation and
experimented results
4 调整结构参数的数值模拟研究
4.1 结构参数介绍
将螺旋进气道的三维模型在CAD软件中直接导出二维工程图并进行适当简化，得出决定螺旋进气道结构形状的主要结构参数，如图3所示.
图3 螺旋进气道主要结构参数
Fig.3 Main structure parameters
of the helical intake port
结构参数包括：涡壳切割量η、气道入口截面积S1、气道最小截面积S2、螺旋底坡角β1、螺旋坡角β2（气道最小截面的法线与水平面的夹角）、螺旋室高度H、螺旋室直径D、螺旋进气道高度μ、涡壳半径R1、R2、螺旋进气道转角θ、螺旋进气道偏心距e以及涡壳转角β3.本文主要研究螺旋室高度H、涡壳切割量η以及螺旋进气道转角θ改变对缸内涡流比的影响.
4.2 参数H的影响
如图4所示，H为螺旋进气道的螺旋室高度，通过三维造型软件改变H的取值，分别取H1=7 mm，H2=7.5 mm，H3=10.5 mm，H4=13 mm以及H5=13.5 mm.计算得出H不同取值时的缸内涡流比，如图5所示.对比H不同取值的进气量，如图6所示.由图5可以看出，在整个压缩过程中，缸内涡流比呈先减小后在上止点附近快速增大，达到一个峰值以后再减小的过程［8］.原因是：在压缩过程中，随着活塞的上移，进气初期形成的缸内涡流受压而减弱，从而使涡流比降低；当接近上止点时，缸内气体被强制压入燃烧室，气体旋转半径减小并产生挤流，从而使缸内涡流增强，涡流比增大并在上止点附近达到最大值；随着活塞下移，缸内气体从燃烧室内流出，气体旋转半径再次增大，缸内涡流减弱，涡流比减小.
图4 螺旋进气道结构参数H
Fig.4 Helical intake port structure parameter H
曲轴转角/°ATDC图5 H不同取值时的涡流比
Fig.5 Swirl ratio of different values of parameter H
曲轴转角/°ATDC图6 H不同取值的进气量
Fig.6 Air inlet of different values of parameter H
同时，可以看出随着H取值的增加，缸内涡流比整体上成先减小后增大再减小的规律，并在H=13 mm时，缸内涡流比获得最大值.可以得出：H与缸内涡流比之间存在非线性的相关关系，存在一个最佳H取值，使涡流比具有最大值.
继续分析图5可知，改变H的取值，对-20°和0°时刻的涡流比影响较大，而对20°时刻缸内涡流比的影响不明显.因此，-20°和0°时刻的缸内涡流比对结构参数的变化较敏感，结构参数的微小变动，都会影响此时的缸内涡流比，而20°时刻的缸内涡流经过了上止点的压缩过程，对结构参数的变化已经不那么敏感了.
由图6可知，虽然H的取值不断增加，但进气量的变化却不明显.在-120°时刻，缸内进气量最大值与最小值之间相差0.72%，因此，可以认为：结构参数H的改变对进气量没有影响.
4.3 参数η的影响
根据实际试验中，涡壳切割量对缸内涡流比有一定的影响，因此，取η为涡壳切割量，如图7所示.通过修改使其取值分别为η0=0 mm，η1=1.27 mm，η2=2.54 mm.
图7 结构参数η
Fig.7 Structure parameter η
计算得出η不同取值时的缸内涡流比，如图8所示.对比η不同取值的进气量，如图9所示.
曲轴转角/°ATDC图8 η不同取值的涡流比
Fig.8 Swirl ratio of different values of parameter η
由图8可以发现，随着η的增大，缸内涡流比先增大后减小.在η=1.27 mm时，涡流比在-180°~20°曲轴转角之间都具有最大值，且其最大涡流比（1.114 8，-4°）比η=0 mm的最大涡流比（0.950 6，-3°）高17.3%；同时，在-20°和0°时刻，η=1.27 mm的涡流比比η=0 mm时高17.9%和17%.随着η值的继续增大，缸内涡流比降低，并低于η=0 mm的涡流比.
曲轴转角/°ATDC图9 η不同取值的进气量
Fig.9 Air inlet of different values of parameter η
由图9可知，虽然η的取值不断增加，但进气量的变化却不明显.在-120°时刻，缸内进气量最大值与最小值之间相差0.58%，因此，可以认为：结构参数η的改变对进气量没有影响.
4.4 参数θ的影响
以进气阀的轴线为轴，进气阀的轴心与进气道入口截面的中心连线，绕此轴旋转的角度为θ，如图10所示.设气道初始位置的θ角度为0°，且顺时针旋转为+，逆时针旋转为-.通过三维建模软件修改θ的取值，分别为-15°，-10°，-5°，+5°，+10°以及+15°.
计算得出缸内涡流比，如图11所示.对比参数θ不同取值的进气量，如图12所示.由图11可以看出，无论是顺时针旋转还是逆时针旋转，缸内涡流都降低，且顺时针旋转降低的程度要远小于逆时针旋转降低的程度，同时随着旋转角度的增大，涡流比下降的速率也逐渐增大.根据不同的θ取值，可以得出初始位置时螺旋进气道所处的θ角是使缸内涡流比为最大值的最佳角度.
图10 结构参数θ
Fig.10 Structure parameter θ
曲轴转角/°ATDC曲轴转角/°ATDC
图11 θ不同取值的涡流比
Fig.11 Swirl ratio of different values of parameter θ
曲轴转角/°ATDC图12 θ不同取值的进气量
Fig.12 Air inlet of different values of parameter θ
从图12可以看出，无论θ的取值如何变化，进气量的变化却不明显.在-120°时刻，缸内进气量最大值与最小值相差0.6%，因此，可以认为：结构参数θ的改变对进气量没有影响.
5 调整多个结构参数的研究
为研究多个参数改变时，对缸内涡流比的影响，分别同时改变2个参数和3个参数，设置5个Case的模拟计算，表4所示，为每个Case中各参数的取值.
表4 各Case参数取值
Tab.4 The parameters values of different Cases
名称
H/mm
η/mm
θ/（°）
Case1
10.5
Case2
13
1.27
Case3
12
1.27
Case4
13
1.27
+10
Case5
13
1.27
-10
其中，Case1的参数取值为所研究发动机螺旋进气道结构的实测值；Case2的参数取值是将各参数调节到使涡流比为最大值的参数取值；Case3的参数取值是保持Case2中η和θ的取值不变，调整H的取值；Case4和Case5的参数取值是保持Case2中H和η的取值不变，调整θ的取值.
按表4所示，修改结构参数后导入计算，得到5个Case的涡流比，如图13所示.对比5个Case的进气量，如图14所示.
曲轴转角/°ATDC图13 5个Case的涡流比
Fig.13 Swirl ratio of five Cases
曲轴转角/°ATDC图14 5个Case的进气量
Fig.14 Air inlet of five Cases
比较图14中的Case1和Case2可以发现，虽然分别单独调节H=13 mm或η=1.27 mm时能使涡流比达到最大值，但当同时调整H和η时，并没有使涡流比得到提高.因此，可以认为：各结构参数对涡流比的影响效果不是简单的相互叠加，而是存在相互制约、相互抵消的作用.在此基础上，对比分析Case2和Case3，在保持η和θ取值不变的情况下，适当调整H的取值时，发现Case3的涡流比得到一定幅度的提高，但整体上仍小于Case1的涡流比.同样，对比分析，Case2和Case4和Case5可以发现，在保持H和η取值不变的情况下，调整θ的取值，涡流比较Case2有增大也有减小，但变化幅度没有调整参数H明显.
通过上面的分析，可以得出：缸内涡流比是各结构参数相互作用的结果，各结构参数之间存在相互制约的作用.欲使缸内涡流比得到优化，必须准确了解各结构参数与涡流比的相关关系，通过分析相关关系可以得出优化缸内涡流比的各参数取值.
分析图14可以看出，在各种情况下进气量的变化都不大.在-120°时刻，缸内进气量最大值与最小值相差1.1%，因此，可以认为：在误差允许的范围内，本文中所讨论的结构参数改变对进气量没有影响.
6 缸内流动情况的分析
分析修改结构参数后对缸内流动情况的影响，沿两进气阀的轴线做一竖直截面，取进气阀在最大升程时刻，且右侧为螺旋进气道，获得缸内流动速度分布图，如图15所示.
(a)Case1(b)Case 5图15 缸内流动速度分布
Fig.15 Distribution of incylinder flows speed
分析图15可以看出：在15(a)图中，气缸中心处产生了一个较大涡流，并进一步使进气气流在两气门之间向左侧流动，加速了缸内涡流的产生；而15(b)图中，由于缸内产生的涡流较小，进气气流进入缸内后直接流向燃烧室，缸内不能很好地组织涡流.分析原因，可以得出：由于Case1和Case5的螺旋进气道在结构设计参数H，η和θ的取值上存在一定的差异，从而使Case1和Case5的进气道在进气过程中产生涡流的能力发生改变.Case1中各结构设计参数的取值优于Case5，从而使Case1的缸内涡流比高于Case5.因此，合理的设计螺旋进气道的结构参数能很好的改善发动机的缸内流动情况，提高缸内涡流比，最终达到优化发动机缸内燃烧过程的目的.
7 结 论
1）本文讨论的螺旋进气道结构参数的改变对进气量的影响较小，而对缸内涡流比的影响较大.当分别单独修改螺旋室高度、涡壳切割量以及螺旋进气道转角为13 mm，1.27 mm以及0°时，缸内涡流比得到大幅度的提高，并取得最大值.
2）分析了发动机缸内涡流比的动态变化情况，研究了对燃烧有意义的时间段（-20°，0°和20°）的缸内涡流比，其结果比稳流试验台的结果更反映发动机的实际工作过程，更真实合理.
3）同时改变多个结构参数对缸内涡流比的影响比较复杂，各参数之间的作用效果不是简单的相互叠加的，而是相互制约的.
4）采用三维计算流体力学软件Converge，可以较方便地实现螺旋进气道多结构参数的数值研究和优化，对气道的设计有重要的指导意义.
致谢
感谢Convergent Science, Inc.提供Converge软件以及相应的技术支持.
参考文献
［1］ 周龙保．内燃机学［M］．3版．北京：机械工业出版社，2011：261-266．
ZHOU Longbao．Internal combustion engine theory［M］．3nd ed．Beijing：China Machine Press，2011：261-266．(In Chinese)
［2］ 蒋德明．内燃机燃烧与排放学［M］．西安：西安交通大学出版社，2001：512-538．
JIANG Deming．Internal combustion engine combustion and emissions［M］．Xi`an：Xi`an Jiao Tong University Press，2001：512-538.(In Chinese)
［3］ 王天友，刘大名，沈捷，等．四气门柴油机进气道开发［J］．内燃机工程，2008，29（2）：51-55．
WANG Tanyou，LIU Daming，SHEN Jie，et al．Development of Intake ports of foulvalve diesel engine［J］．Chinese Internal Combustion Engine Engineering，2008，29（2）：51-55.(In Chinese)
［4］ RICHARDS K J，SENECAL PK, POMRANING E.CONVERGE (Version 1.3)［M］.Middleton：ConvergentScience，Inc,2008．
［5］ 吴克启，舒朝晖．高等流体力学［M］．北京：中国电力出版社，2009：13-26．
WU Keqi，SHU Chaohui．Advanced fluid mechanics［M］．Beijing：China Electric Power Press，2009：13-26．(In Chinese)
［6］ 周俊杰，徐国权．FLUENT工程技术与实例分析［M］．北京：中国水利水电出版社，2010：14-17．
ZHOU Junjie，XU Guoquan．FLUENT engineering techniques and examples analysis［J］．Beijing：China WaterPower Press，2010：14-17(In Chinese)．
［7］ 邱卓丹，沈捷．直喷式柴油机螺旋进气道性能试验及评价方法［J］．内燃机工程，2005，26（3）：27-30．
关键词: 改善; 模拟；气流组织; 分析; 比较;
中图分类号： F407.474文献标识码： A
1 机舱物理模型
机舱内设备及管线众多，结构复杂，建立模型时必须进行简化来使模型适合数值模拟计算，模型简化的基本原则是: 对空间气流场产生作用较小且温度影响不大的设备管线进行移除; 对结构形状复杂，简化后对模拟结果影响较小的设备规则化; 对厚度边界做无厚度壁面处理; 忽略风管，在风口段用立方体代替，底面为风口;抽风围阱对流场影响小，以风口面代替等;该船舶机舱模型空间大小为: 长 21m，宽 13m，高 10m，主要设备有汽轮机、锅炉、涡轮增压机组、高温除氧器、等离子过滤器、部分烟气管道等；该舱底层为铺板层，沿高度方向分别为第 1、2 格栅层，顶部为甲板层，风口均布置在第1、2 格栅层下方，为防止送风直吹发热表面，造成较大热应力，风口均采用矩形风口，风向垂直向下，尺寸大小为0.4m×0.2m 的风口55个，尺寸0.4m×0.28m 的风口 4 个 ( 侧 45°方向增加局部绕流) 。机舱简化后几何模型如图 1 所示。
图 1 机舱几何模型
2 流体数学模型及边界条件
2.1 流体数学模型
CFD 技术经过长时间的发展，在计算流体流动传热传质方面已经非常成熟，本文通过 Fluent软件CFD技术模拟机舱内气流组织分布，机舱内通风为湍流对流换热问题，本文选用图1模型，为获得较好的模拟结果，提高模拟效率，对模拟进行如下假设和简化:
1) 机舱内流体为不可压流体，密度符合假设;
2) 舱内流体流动及传热视为稳态过程;
3) 各壁面的辐射传热不计;
2.2 边界条件设置
利用开局让棋法进行几何建模及划分网格，因机舱内结构复杂，模型网格采用非结构化网格，对风口及重要发热面进行局部加密，网格数量为 88 万.
舱内通风是采用典型的置换通风与自然通风相结合的形式，边界条件设置如下:
1) 送、回风口边界条件: 送风及机械抽风口为速度入口边界，送风方向与出风口面垂直；
2) 自然排风口边界条件: 通风时保证舱内正压，因此自然排风口采用压力出口边界；
3) 固体壁面边界条件: 舱内发热固体壁面采用定热流密度，壁面厚度为0.
3 设计方案及模拟结果分析
3.1 设计方案
机舱内送风口位置相对固定，采用下送上回的通风方式，回风进入抽风围阱后通过管道与各风机相通；本文设计锅炉燃烧空气量占总送风量50% 左右，从舱内吸入会造成内部各区域气流压力梯度过大，因此采用独立系统从外界抽进燃烧空气。根据《ISO8861-1998造船--柴油机船舶机舱通风设计要求和计算基准》进行热负荷及通风量估算，得到该通风系统排除舱内余热所需风量为160000m3/h，舱内设计温度54℃，送风温度39℃，风机开度 100%时，送风口的设置送风速度为9.4m/s，考虑到送风量对气流组织及舱内温度的影响，设定送风开度来进行分析，分别为 100%，85%，70%.机舱环境的优化可通过改变送抽风口位置形状、面积及送风量等来实现，为获得最佳气流组织形式，该舱通过增大抽风围阱及隔屏来增加舱内空气扰流，防止大漩涡及通风死角，理论上，机械抽风口分布越多，抽风围阱尺寸越大，会有效改善舱内气流组织，但是机舱内空间狭小，因风口而增加的风管会减少舱内可用空间; 且过多地增加隔屏会给工作人员管理维护带来不便，经过分析给定 4 种调整方案来进行研究，见表 1
3.2 模拟结果分析
3.2.1 气流组织分析
计算收敛后，选取典型截面进行分析，由图 1可以看出，两锅炉间至主机右侧面区域散热面积最大，流场也最复杂，因此截取y=﹣0.5m 面进行气流组织分析，y=﹣0.5m 截面如图 2所示.
在假定机舱送风口位置固定情况下，要消除漩涡和通风死角，一般采用下列方法: 一是设置障碍物，障碍物对室内气流组织有显著影响; 二是布置空气射流通风系统，由于气流至底层后从下至上流动，障碍物迎风面不存在滞留区，同时考虑经济成本问题，对主机右侧安装隔屏作为障碍物来消除此漩涡，隔屏位置如图 2 所示"
4 结 语
模拟结果表明，利用CFD技术对机舱气流组织进行模拟分析是一种高效简洁的方法，获得的结果是可信的，且通过模拟分析可知:
1) 合理的布置机械抽风口能有效地改进舱内流场，在舱内增加隔屏作为障碍物并非完全是负面影响，在不影响舱内设备运行和人员工作情况下，增加隔屏能有效改变舱内气流组织
2) 自上而下的通风方式能使各工作面温度得到较好的控制，人员活动区域温度符合设计要求
3) 风机开度直接影响舱内温度，通过改变风机开度能节省通风资源，机舱内环境的优化有很多方式，可通过改变风口数量大小、形式、位置等，这些还需要进一步的研究和补充
参考文献:
［1］ 向立平，王汉青空调客车内气流组织与污染物浓度场数值模拟［c］中南大学学报( 自然科学版) ，
关键词：拱形支护结构；空间特性；内力；变形
Based on the MATLAB numerical analysis of arched retaining structure is analyzed
Li Xin
（HeBei Ruizhi Transportation Technology Consulting Co., Ltd
HeBei Provincial Engineering Eesearch Center for Road Structure and Material
Shijiazhuang050091）
Abstrct: Arched retaining structure of deep foundation pit can give full play to its space effect, the uniform stress, the cross section only compressive stress, which can effectively control the deformation and ensure the stability of foundation pit. Engineering on the arch supporting system, using the shell theory to establish the differential equations, and then calculate the internal force and deformation of retaining structure, due to considering the spatial characteristics of this kind of supporting structure, so this method is more in line with the actual working situation.
Key words: arch supporting structure; spatial features; the internal force; deformation
中图分类号： TV551 文献标识码： A
1前言
拱结构拥有良好的力学性能，可以将上部荷载沿着拱圈传递至拱脚，从而在拱圈内拉应力值很小，甚至不出现拉应力[1~2]。拱形支护结构在基坑开挖支护中应用较为广泛，这种结构体系与常规的直线形支护结构相比，具有较多优点，工程实践已表明采用这种支护形式后可大大减低基坑支护的造价，经济效益相当可观[4~5]。
2结构计算模式
拱形支护结构通常由几段曲率不同的拱组成，在每个拱与其他结构连接处（拱脚）受力最大，故拱脚处往往得到加强，每一段拱从空间上看都可认为是直线边支承的壳体结构，它所受的荷载为土压力，方向为壳体的法向，见图1所示。圆柱壳的内力分布形式见图2所示。在基坑工程中，圆拱使用较多，现以圆柱壳作为支护结构的计算模型。对于只承受法向荷载的圆柱壳，其基本微分方程为式（1）：
（1）
其中：分别为圆柱壳的母线、环向坐标；分别为圆柱壳在母线、环向及法向上的位移；E为壳体材料的杨氏模量；为壳体材料的泊松比；R为圆柱壳的半径；t为圆柱壳厚度；P为壳体所受的法向荷载（即土压力）。
图l拱形支护结构计算模型
图2（a） 图2（b）
图2拱形支护结构计算模型
引入位移函数F（），则可将壳体的位移表示成式（2）：
（2）
将式（2）代人式（l），前两式满足，由第三式可得：
（3）
式中，
对于拱形支护结构，把拱脚处假设为简支端，则位移函数F（）可表示为单三角级数形式：
（4）
式中，m为大于1的整数。
将法向荷载P（）也展开成三角级数形式为
（5）
将式（5）及式（4）代人式（3），比较方程两边的系数，可得
（6）
式中，。
对于拱形支护结构，被动区的土压力可以用弹簧模拟，即，其中为水平基床系数；w为支护结构的位移。
主动区的土压力采用主动土压力，因此对于拱形支护结构其主动土压力在母线方向上是线性变化的，而在环向上保持不变。
壳体支护结构上作用的荷载为：P=，其中h为基坑深度；c，d为与土压力有关的量。
支护结构在坑底以上部分的微分方程为：
（7）
壳体的内力由位移函数表达为：
根据基坑开挖面处（a=h）支护结构的连续性，有：
位移连续：
力的连续：
其中上标为u的符号表示开挖面以上的量，上标为d的符号表示开挖面以下的量。根据基坑顶部和底部的边界条件，得：
边界为固支时：；
边界为间支时：；
边界为自由端时：，。
由连续条件及边界条件即可求出式（4）的位移函数，根据式（2）和位移函数表达式可计算出支护结构各点的位移及内力。
3工程实例
某深基坑为局部扇形结构，开挖深度6.8m，采用由单排650钻孔灌注桩组成的拱形支护结构，圆拱半径为25m，圆心角为120o，在两边拱脚处各布置10根650钻孔灌注桩，桩侧加打水泥搅拌桩以加固周围土体，形成支撑体系。支护结构平面布置及剖面见图3所示。该深基坑主要土层自上而下分布情况为：杂填土，厚度1.5m；粘土，厚度1.0m；淤质粘土，厚度11.0m；粉质粘土，厚度3.2m；黄褐色粘土，厚度4.5m。为了监测基坑附近土体的侧向位移，在拱顶处设置测斜管进行现场实时实测。
图3支护结构平面及剖面
根据土层分布及地下水位情况，计算得到土水压力分布情况。利用MATLAB数值分析，得到拱顶处的位移（w）计算值与实测结果，见图4所示；拱顶弯矩（M）计算结果，见图5所示。监测结果表明，尽管为单排桩组成的自立式支护结构，但由于采用了拱的结构形式，使整个支护结构的侧向位移和内力均较小。
图4拱顶处侧向位移计算域实测值比较
图5拱顶处弯矩计算值
由图3可知，拱顶处的侧向位移较小，最大值仅为3.6cm，侧向位移在0~6m范围内随深度的增加而增加，6m处的侧向位移最大；6m以下范围内侧向位移随深度的增加而减小。由图4~图5可知，拱顶处的弯矩绝对值也较小，最大值仅120kN.m，0~5m范围内负弯矩绝对值随深度增加而增大，5m以下范围负弯矩绝对值随深度增加而减小，10~15m之间出现反弯点，反弯点以下范围内正弯矩随深度增加而增大。
4结语
通过上述分析，得到以下几个方面的结论与建议：
（1）拱形支护结构是一种合理的支护形式，在本文介绍的非闭合情况下，拱顶处的侧向位移及内力均较小，同时又充分发挥材料的抗压性能，比常规直线形支护结构具有较大的优越性。
（2）拱形支护结构可以把支护形式由单拱向多拱方向发展，推动拱形支护结构向更大、更高结构层次发展和应用。在材料选用上，可以向更为广泛的钢筋混凝土灌注桩、地下连续墙、钢桩支护等支护方式上拓展拱形结构，有效降低工程成本、节约工期。
（3）由于拱形支护结构固有的外形特点，使得拱形支护结构在施工场地狭小的工程在使用上受到限制。
参考文献：
[1] 陆瑞明,赵锡宏.挡土拱圈的内力和变形计算[J].岩土工程学报,1997,19(6).
[2] 张祖闻,陈新余.高层建筑地下结构及基坑支护[M].宇航出版社,1994.
[3] 吴连元.板壳理论[M].上海交通大学出版社.1989,12~15
[4] 梁春江.连拱式基坑支护与实践[J].施工技术.2002.
[5] 房集明.连拱形基坑支护体系[J].岩土工程技术.2000
[6] 龚晓南等.深基坑工程设计施工手册[M].中国建筑工业出版社.1998
[7] 张士乔等.基于拱梁法原理的深基坑拱形围护结构分析[J].土木工程学报.2001
[8] 黄生根.连拱式水泥土支护结构的应用[J].岩土工程技术.1998
作者详细信息：李欣，女，（1986―），助理工程师；
