水动力-力学耦合模型（精选4篇）
水动力-力学耦合模型 篇1
1 引言
随着现代机械向高速、重载、高精度方向发展, 机械系统的结构逐渐多样化。采用构件刚性假设的传统多刚体系统动力学模型逐渐不适用于越来越复杂的机械系统, 其分析结果会产生一定的偏差以至于不能满足精度要求, 因此考虑系统中构件的弹性变形已成为多体系统动力学的重要研究方向之一[1,2]。
ANSYS是一款大型通用有限元分析软件, ADAMS主要用于多刚体运动学与动力学仿真, 将两款软件结合使用, 即ANSYS-ADAMS联合仿真分析是解决复杂机械系统动力学问题的一种有效的方法[3,4]。
本文基于柔性体理论, 以3-PUS并联机构为研究对象, 采用ANSYS-ADAMS联合仿真的分析方法, 对3-PUS并联机构的刚柔耦合模型进行了仿真分析, 以了解这一机构的动力学特性, 为进一步优化提供理论依据。
1 柔性体理论
1.1 柔性体系统建模理论基础
ADAMS用模态来表示物体的弹性, 这种弹性变形近似地用模态的线性叠加来表示, 刚体的位姿用反映刚体位置的笛卡尔坐标x= (x, y, z) 和反映刚体姿态的欧拉角u= (u, v, w) 表示, 模态坐标用q={q1, q2…qi} (i为模态坐标数) 表示。在完全约束条件下, 唯一能够确定质点系位置的独立变量可以用作广义坐标[5], 则柔性体的广义坐标表示为
根据选取的坐标系的不同, 柔性体建模可以分为3种方法, 包括浮动坐标系法、惯性坐标系法和随转坐标系法[6]。本文采用应用最为广泛的浮动坐标系法建立柔性体的动力学模型。
柔性体系统中建立静坐标系 (er) 和浮动坐标系 (eb) , 如图1所示。静坐标系是固定不变的, 作为参考坐标系, 浮动坐标系建立在柔性体上, 以其相对于静坐标系的平移和转动来表述柔性体的位姿。柔性体相对于静坐标系的运动可以分解为:刚体平移、刚体转动和柔性体的变形运动, 即柔性体上任意一点的位置的向量表示为
式中:ri为柔性体中一点在静坐标系中的向量, A为动坐标系到静坐标系的变换矩阵, rr为浮动坐标系原点在静坐标系中的向量, rb为柔性体中一点在浮动坐标系中的向量, up为相对变形量, up=φpζ, 其中φp为点p满足里兹基向量要求的假设变形模态矩阵, ζ为变形的广义坐标。
1.2 柔性体系统的动力学方程
柔性体的动能表达式为
式中, M (ξ) 为质量矩阵, 是一个3×3维的矩阵。柔性体的势能包括重力势能和弹性势能两部分, 其表达式为
式中:Wg (ξ) 为重力势能;2/1ξTξ为弹性势能;K为模态刚度矩阵, 是常值矩阵。
能量损耗函数为
式中:为模态坐标对时间的一阶导数;矩阵D为模态阻尼矩阵, 为常值对称阵。
为得到柔性体系统动力学方程, 将式 (3) ~式 (5) 代入下面的拉格朗日方程:
式中:L为拉格朗日项, L=T-W;Ψ为约束方程。得到柔性体系统的动力学方程为
2 3-PUS并联机构刚柔耦合模型仿真[7]
2.1 构建3-PUS并联机构的多刚体模型
本文研究的3-PUS并联机构包括动平台、静平台、3组连杆滑块和机架, 滑块和机架之间为滑动副, 滑块和连杆之间通过虎克铰连接, 连杆和动平台之间通过球铰相连, 通过电机传动使3个滑块沿导轨直线运动, 进而带动动平台做3个平移运动。将UG NX7.5中建立的3-PUS并联机构模型以Parasilid (*.x_t) 文件格式导入ADAMS中, 在ADAMS中将3-PUS并联机构的材料定义为45钢, 添加正确的运动副, 完成后并联机构的多刚体模型如图2。
2.2 构建3-PUS并联机构的刚柔耦合模型
在本文研究的3-PUS并联机构中, 其连杆的长度和直径的比值较大, 在工作中易发生弹性变形, 因此不能忽略连杆的柔性, 将ADAMS中3-PUS并联机构中原来的刚性连杆替换为柔性连杆并建立3-PUS并联机构的刚柔耦合模型再进行动力学仿真, 可得出较为精确的仿真结果。构建3-PUS并联机构刚柔耦合模型的方法如下:
1) 将UG NX 7.5中的3-PUS并联机构的三维装配模型中的连杆模型以Parasolid (*.x_t) 文件格式输出, ANSYS通过专门的接口读入Parasolid模型文件;
2) 连杆的材料为45钢, 在ANSYS中定义材料的弹性模量E=2.07×1011Pa, 泊松比PRXY=0.27, 密度ρ=7 801 kg/m3;单元类型选择Solid45和mass21;
3) 选用Solid45单元, 使用智能网格划分功能对连杆模型进行网格划分;
4) 根据ADAMS中3-PUS并联机构中构件的连接关系, 为ANSYS中的连杆模型添加Keypoints, 并使用mass21对Keypoints进行网格划分, 生成外部连接点;在ANSYS主菜单中选择“Preprocessor→Coupling/Ceqn→Rigid Region”创建刚性连接区域。
5) 在ANSYS的主菜单中选择“Solution→ADAMS Connection→Export to ADAMS”, 输出的连杆模型的模态中性文件, 如图3。
6) 在ADAMS中, 将3-PUS并联机构中的6个连杆替换为柔性体连杆, 替换后柔性体连杆会继承原来刚性连杆上的运动副、载荷和Marker点等特征, 图4为3-PUS并联机构的刚柔耦合模型。
2.3 3-PUS并联机构刚柔耦合模型的仿真
在ADAMS环境中为图4左侧滑块移动副添加滑移驱动, 使该滑块在前1 s以50 mm/s2的加速度做匀加速运动, 在1~5 s以50 mm/s的速度做匀速运动, 其他2个滑块保持静止, 驱动函数为IF (time-1:50*time, 50, 50) 。
表1为样机模型在ADAMS中的仿真参数设置, 仿真I的仿真时间为5s, 以观察仿真结果的整体趋势, 便于进行全面的分析;仿真II的仿真时间为0.1s, 以观察仿真模型在启动瞬间的运动特性, 是对仿真I的分析结果的补充和深化。
仿真运动中的并联机构如图5所示, 可以看出运动中连杆上Von Mises应力最大值为76.2 MPa, 而45钢的屈服极限为300 MPa, 可知该机构运行安全。
3 仿真结果与分析
多刚体模型与刚柔耦合模型仿真均可得到速度、位移和应力随时间的变化曲线图, 但刚柔耦合模型仿真的结果, 考虑了部件的弹性变形等实际物理特性, 因此反映的运动学和动力学特性更为准确。为说明刚柔耦合这一特点, 本文将3-PUS并联机构的刚柔耦合模型 (部件中含柔性体) 仿真结果和多刚体模型 (部件中均为刚体) 的仿真结果进行对比说明。
3.1 运动学仿真结果分析
仿真模型建好后, 按规划的路径驱动进行仿真, 仿真结束后, 在ADAMS/Postprocessor模块中观察仿真结果图 (图6、图7) , 图中Part2为并联机构的动平台。
图6 (a) 、图6 (b) 分别为仿真I中刚柔耦合模型和多刚体模型中动平台的速度和位移图。
图7 (a) 、图7 (b) 分别为仿真I中刚柔耦合模型和多刚体模型中动平台的X向速度和位移图 (动平台的Y向和Z向的速度和位移对比图与X向的速度和位移对比图与图7类似, 这里不再给出) 。
对比分析图6和图7, 可得出以下结论:
1) 比较两模型速度曲线可以看出, 速度曲线基本一致, 说明刚柔耦合模型的动平台的速度与多刚体模型的动平台的速度均值及整体变化趋势一致;但在仿真运动的开始前0.5 s内, 刚柔耦合模型中动平台的速度曲线有较为剧烈的波动现象。
2) 比较两模型位移曲线可以看出, 图6、图7中刚柔耦合模型和多刚体模型的位移曲线基本相同。
这里需要分析的问题是速度波动, 应该反映在位移曲线上 (位移是速度的积分) , 但仿真的结果两种模型的位移曲线却基本相同。原因可能是刚柔耦合模型在速度上的短时波动变化, 由于仿真步数的限制, 未能反映在更小时间单位上位移曲线存在波动的真实情况。
为了得到刚柔耦合模型中动平台的真实位移曲线和更清晰的速度曲线, 以确认仿真结果的正确性, 将仿真时间缩小为0.1 s、步长为20 000即仿真II, 以进一步观察和分析两者的动平台的速度和位移曲线的对比图 (如图8和图9) 。
对比分析图8和图9, 可得出以下结论:
1) 速度曲线表明:在仿真运动的初始时刻, 刚柔耦合模型中动平台的速度曲线呈周期约为0.01 s的简谐波, 并且其振幅逐渐衰减, 在t=0时其速度为0, 在t=0.001 8 s时刻速度就达到了54.022 5 mm/s, 而多刚体模型中动平台的速度t=0.001 8 s时速度仅为0.116 3 mm/s。
刚柔耦合模型仿真结果说明在实际运动的初始时刻, 动平台的速度以较高的频率和较大的振幅振荡, 并联机构动平台由静止到运动瞬间存在一定的振动, 其原因是杆件存在弹性变形。依据振动的特性与波形, 可以进一步地优化杆件结构尺寸以减小速度曲线的频率和振幅, 从而使并联机构的启动阶段更加平稳。
2) 位移曲线表明:在图8 (b) 中可以看出刚柔耦合模型动平台的位移曲线的斜率随速度的波动出现变化, 图10 (a) 中表示在前0.1 s内, 动平台的X向位移曲线的斜率随X向速度的波动也出现变化, 动平台的位移曲线也呈现一定频率的波形。
图6、图7中, 在0~5s时间内速度发生一定频率的波动, 但由于时间步长的原因, 不能在位移曲线发生相应的波动。而图8 (b) 、图9 (a) 在更小时间单位上, 刚柔耦合模型中动平台位移曲线证明了波动的存在, 从而反映了动平台的实际运动特性。
3.2动力学仿真结果分析
图10为仿真I中刚柔耦合模型和多刚体模型中球铰5处 (图4右前位置处) 沿X向、Y向、Z向受到的力及合力的对比图;图11为仿真II中球铰5处沿X向、Y向、Z向受到的力及合力对比图。
对比分析图10得知, 在仿真运动开始的0.5s, 刚柔耦合模型中球铰5处受到的力均以较高的频率和较大的振幅往复振荡, 之后球铰5处的力的变化趋于稳定, 总体上刚柔耦合模型和多刚体模型中球铰5处的力的变化趋势一致, 但是多刚体模型中球铰5的力学曲线并未反映出在运动的初始时刻球铰5处的力的振动现象, 在图11中可以看出刚柔耦合模型中球铰5处的受力曲线呈振幅逐渐衰减的简谐波形, 周期约为0.01 s, 在t=0.015 s处的力达到了最大值为2 886.735 4 N, 而刚体模型中球铰5在仿真运动的5s内最大力为1 531.419 7 N, 其原因在于在仿真运动初始时刻, 因惯性的存在, 刚柔耦合模型中的柔性连杆会发生弹性变形, 因此在运动初始时刻球铰5的应力变化更大。
动力学仿真结果进一步说明刚柔耦合模型可以更准确地反映部件受力情况, 从而为部件产品设计、优化和校核提供依据。
4 结论
1) 应用ADAMS等工程软件构建复杂机械系统的多刚体模型和刚柔耦合模型进行运动学和动力学仿真和求解, 为复杂机械系统的设计分析和优化提供了一种实用的方法。
2) 在建立虚拟样机时, 对含有长细比较大或在实际工况中易发生变形的部件的机械系统, 应考虑构建刚柔耦合的虚拟样机, 以获得更为真实有效的仿真结果。
3) 采用ANSYS和ADAMS联合仿真的分析方法, 对3-PUS的并联机构的刚柔耦合模型和多刚体模型进行了仿真, 分别将两者的结果进行了整体和局部的对比, 发现该并联机构在启动瞬间, 速度、位移、应力等均存在一定的波动, 刚柔耦合模型仿真结果反映了该并联机构真实的动态特性, 为并联机构的设计和优化提供了理论依据和方法指导。
摘要：以柔性体理论为基础, 使用多体动力学仿真软件A D A M S和有限元分析软件A N SY S建立了3-PU S并联机构的刚柔耦合模型, 对3-PU S并联机构的刚柔耦合模型和多刚体模型进行了仿真和分析。结果表明, 在仿真运动开始的瞬间, 刚柔耦合模型中动平台的速度、位移和铰接处的应力均有一定频率的波动现象, 反映了该并联机构实际的动态特性, 为并联机构的进一步的设计、优化和校核提供了依据。
关键词：并联机构,ADAMS,ANSYS,刚柔耦合模型
参考文献
[1]Shabana A A.Flexible multibody dynamics:review of past and recent developments[J].Journal of Multibodt System Dynamics, 1997 (2) :189-222.
[2]刘铸永, 洪嘉振.柔性多体系统动力学研究现状与展望[J].计算机力学学报, 2008, 25 (4) :411-415.
[3]张雷, 余克龙, 陈少钟, 等.基于ANSYS和ADAMS空间四连杆引纬机构的柔性动力学仿真[J].纺织学报, 2013, 34 (5) :116-120.
[4]张永德, 汪洋涛, 王沫楠, 等.基于ANSYS与ADAMS的柔性体联合仿真[J].系统仿真学报, 2008, 20 (17) :4501-4504.
[5]王爱勤, 冯振宇.也谈虚位移原理中的广义坐标的选择[J].大学物理, 2003, 22 (8) :13-14.
[6]Wasfy T M, Noor A K.Computational strategies for flexible multibody system[J].Appl Mech Rev, 2003, 56 (6) :553-613.
[7]李增刚.ADAMS入门详解与实例[M].北京:国防工业出版社, 2013.
水动力-力学耦合模型 篇2
烽火村乔木湾岩溶塌陷位于武汉市洪山区青菱乡烽火村二、三组,上、下倒口湖之间,东距京广铁路750m,西距107国道120m,距长江直线距离1.8km。1997-2005年先后共发生大小陷坑23处,且近几年有新的塌陷坑形成,陷坑平面形态呈椭圆形,烽火村塌陷坑直径一般<10m,深跨比≥0.3,认为是相互连通的扰动带和多个土洞发育而引起的塌陷。
垂向上分成三个发育段。
(1)上发育层。溶洞埋深20~40m,标高段一般为(0~20m),钻孔遇洞率15.71%,此埋深段为研究区溶洞垂向发育程度最强。几乎存在于岩溶强发育区及较强发育地区,分布较广泛,面状多发育为溶孔、溶槽的形式,此层溶洞发育规模一般都比较大,最大溶洞高度为7.8m。
上层发育带分为四段:20~25,25~30,30~35,35~40m;各埋深段的岩溶发育情况如图1所示。从图中可以看出20~25m埋深段是一个弱发育段,25~30m为岩溶最发育段。
(2)中发育层。埋深在40~65m,标高段一般为-20~45m,该段岩层溶洞较发育,所有的石灰岩区普遍存在,遍布青菱乡、毛坦港、阮家巷、烽火村等地区,溶洞规模不大,洞高一般在3m以下,全部分布在塌陷区,表明地表塌陷较深层次的溶洞也起着至关重要的作用。
但50~55m埋深段、60~65m埋深段钻孔遇洞率在3%以内,岩溶及溶洞发育相对较少,为岩溶弱发育段。
(3)下发育层。埋深约65 m以下(标高段一般为-45 以下),该段溶洞发育很不均匀,断层带溶洞发育规模较大,最高溶洞达7.08m,这与地层的构造有密切的关系,此标高段钻孔遇洞率为15.6%,分布在阮家巷及毛坦港附近,全部为重点塌陷区,区内深部岩溶应该较发育。
2 数学模型
2.1 地下水流动数学模型
根据Darcy定律和水均衡原理,忽略密度变化对地下水影响的三维流动数学模型为:
式中:h(x,y,z,t)为水位,L;Kxx、Kyy和Kzz分别为沿x、y、z主方向的渗透系数,L/T;w为单位体积单位时间内注入或抽出的水量,1/T;μs为单位储水系数,1/L;t为时间,T;h0(x,y,z)为流场初始水位,L;B1为一类边界条件;B2为二类边界条件;q(x,y,z,t)为二类边界上单位面积已知流量,L/T。
2.2 力学模型的简化
假定地层在土压力变化过程中引起的垂向总应力不变,土体侧向变形忽略,只发生垂向变形。土体的有效应力变化及垂向应变为:
有效应力变化为:
平均有效应力变化为:
垂向应变为:
式中:ΔP为压力增量;Δσx′、Δσy′和 Δσz′分别为x,y和z三个方向的有效应力变化;ΔσM′为平均有效应力,σM′=(σx′+σy′+σz′)/3;v为泊松比;K为体积模量;αp为Biot系数;εz为垂向应变。由此模型可以看出,孔隙水压力的增加引起水平有效应力减小;由于地表的自由变形,垂向有效应力只受孔隙水压力变化的影响;垂向应力应变满足线弹性的Hooke定律。
2.3 耦合方法
在初始流场时有相对应初始应力值,地下水位上升或下降都会引起土体有效应力的变化,利用GMS软件模拟出不同时段地下水的流场值,计算出不同时段的应力,不同时段之间应力差值是压力的增量为 ΔP,不同时段之间有效应力差值为Δσx′、Δσy′和Δσz′,再计算出不同时段土体的应力应变。饱和土体的有效应力表达式为:
2.4 土体破坏判别准则
土体是否产生应力破坏,采用库伦-莫尔准则:
首先需要根据应力状态得到应力莫尔圆,然后确定主应力(有效应力)方向,再得到最有可能发生剪切破坏的方向,最后采用破坏准则判定岩土体是否破坏。
式中:σ1、σ3最大、最小主应力,上式表明岩土体处于极限平衡状态,还有另外两种情况存在。可以得到如下判别关系式:
每一个单元体相应的最大应力及最小主应力为:
3 数值模型
3.1 模型边界条件及网格剖分
(1)侧向边界。研究区西临长江,可设此边界为河流边界,根据资料分析,长江历年的平均水位16.70m,因长江切割了砂层及基岩,且与承压含水层及岩溶水有水力联系,在汛期侧渗补给前缘孔隙水,通过越流形式补给岩溶水,所以砂层以下定为隔水边界,砂层以上至长江底部高程以下定为定流量边界,分段量化补给量。模拟区的东、南、北面边界定为定水头边界,根据观测孔的资料分析得到。
(2)底边界。选在灰岩埋深165 m处,此边界基本没有产生水力交换,可以作为隔水层边界处理。
(3)上边界。潜水层的自由水面作为系统的上边界,作为潜水面边界,系统与外界发生垂向水流交换。如:接受河流渗漏补给、农田灌溉补给、降雨入渗补给、潜水面的蒸发排泄等。各含水层的顶、底板高程根据钻孔资料确定。
通过对该研究区的地质结构模型、地下水水流特征及补、径、排关系分析,将盖层地下水流模型概化为非均质各向同性、三维非稳定流模型。基岩层定义为非均质各向异性、三维非稳定流模型。
把研究区共剖分网格单元6 315个(图2)。垂向上由各含水层底板埋深及厚度决定,总模拟深度180m。
3.2 模型参数分区及参数值
(1)渗透系数、弹性储水系数、给水度。渗透系数、弹性储水系数、给水度是正确建立水流模型的关键,在空间上具有很强的变异性,本次建立水流模型采用等效厚度法,结合钻孔、地质剖面图、观测孔资料、研究区前人的研究成果确定初始参数值,在模型识别阶段进一步调整,图3为各含水层参数分区图。灰岩及其他基岩的参数根据灰岩的埋藏及覆盖情况进行分区。
(2)力学参数初次选取。力学参数有压缩体积模量、回弹体积模量、天然重度、饱和重度、泊松比、黏聚力、摩擦角等。
3.3 模型识别
模型识别的目的是为了提高模型的精度,模型识别的过程是反复调整模型参数(包括边界条件)使计算结果尽可能与实际监测资料一致。模型识别资料选用2009年1月全区的统测水位作为模型识别的初始流场,2012年1月和7月的数据作为拟合检验。
(1)流场拟合。根据2009年1月为初始流场,计算出位于土洞上层即第三层的初始应变值及破坏分布如图4、5所示。
1月是武汉地区的枯水期、7 月是丰水期,长江水位的上涨、降雨集中在7月,2012年1月流场实测值与计算值比较(图6)、7月土体应变及破坏区的分布情况如图7、8所示。
(2)应力及变形拟合。2012年监测点土压力及土体变形拟合效果如图9、10所示。
(3)模拟结果分析。土体压力自2009年开始模拟计算,模拟数据表明,2009年1月第三层砂层处于破坏状态,这与2009年发生在烽火村的岩溶塌陷时间一致,2012年7月模拟的土体破坏区与地下水流场变化有着密切的关系,与实际的塌陷点位置相符。充分说明利用流场值模拟土体的变形及岩溶塌陷是可行的。
土体受力状况的变化值受地下水的影响较大,每年的7-8月份相对于全年的土压力值及有效应力值偏低,是由于地下水位明显升高,产生浮托力作用的结果。在模拟的年度内土压力值的震荡性较大,随着地下水位的变化而压实回弹,经过几次回合后,土体出现塑性变形,之后发展到破坏引起岩溶塌陷。
土体垂向变形与地下水的升幅有直接的联系,土体的沉降与回弹与每一层地下的流场走向相符,在没有考虑施加荷载时,土体沉降的变形值一般小于4mm。
土体破坏规律:地下水的活动和运移,将对土层产生潜蚀作用,从而形成土洞。此外,土洞形成后,其洞壁周围将产生应力集中现象,当地下水位发生变化时(上升或下降时),将进一步改变洞壁周围土体的应力状态,并致使洞体周边产生破坏最终导致塌陷。
(4)识别后的参数,通过模型识别,得到的模型参数如表1-5所示。
4 岩溶塌陷影响因素模拟评价
4.1 长江水位动态对岩溶塌陷的影响预测
实际调查监测表明,研究区内地下水水位动态受长江水位波动影响明显,地下水位动态与长江水位具有显著的相关性。因此,以长江水位变幅较大的2010年水位动态数据,在其他各项参数及源汇项不变的情况下,模拟预测了长江水位动态对岩溶塌陷的影响。
经过两次交替长江水位上升、下降,土体也经历了回弹、压缩作用,预测的第三层(砂层)应变及土体破坏情况如图11、12所示。结果表明,研究区地下水与长江水位具有很好的连通性,当长江水位在最高水位(年)与最底水位(低)交替出现时,地面塌陷点向长江方向偏移,这与地下水的排泄方向一致,岩溶塌陷很明显是由地下水渗流作用引起。
4.2 工程建设排水施工对岩溶塌陷的影响预测
研究区近几年拆迁还建、市政建设等工程大量增加,很多工程涉及到深基坑开挖,为保证施工安全,需要降低地下水。根据施工规范规定,地下水位必须降至基础以下1~1.5m,一般情况下基础的开挖深度10m左右,研究区地下水位需降至高程11m位置。为此,利用数值模型模拟了当某基坑开挖排水时得到的垂向应变及土体破坏情况如图13、14所示。模拟结果显示,土体的垂向应变发生明显变化,从降落漏斗最深处开始破坏。
4.3 大气降水对岩溶塌陷的影响预测
降雨是研究区地面塌陷发生的重要影响因素之一,实际调查表明,研究区内多起地面塌陷的触发因素都是降雨。降雨对地面塌陷的直接作用主要体现在两个方面:使上部土体饱水自重增加,物理力学性质降低;降雨入渗的渗流作用破坏土洞的稳定性。
根据武汉地区的降雨特征,在7月份流场的基础上,通过改变降雨入渗补给量,模拟3种降雨模式(表6),计算土压力的变化值分析降雨对土体变形及岩溶塌陷的影响。
在补给条件不变的情况下,以上几种降雨模式的土体变形及破坏情况如图15-20。
模拟结果表明,研究区降雨对于土洞稳定性有明显影响。三种模式下应力分布规律大致相同,由于区内地下水与长江有着较直接的排泄途径,短期的降雨没有造成水位大幅度上升,但使水压力增大导致土洞破坏。由以上模拟分析的变形位移图可知,相同时间下,降雨强度越大,土洞变形越大,同时降雨强度较小但持续时间较长时,对土洞的影响几乎和降雨强度较大的暴雨相同,也可成为塌陷的主要诱发因素。
5 结论
(1)基于地下水动力学理论和Terzaghi固结沉降理论,采用MODFLOW软件和自编Terzaghi模型计算程序,建立了青菱乡烽火村岩溶塌陷区地下水动力学- 岩土力学耦合数值模型。
(2)利用该模型模拟了长江水位动态、基坑排水、降雨等因素对岩溶塌陷的影响,结果表明:长江水位上升、下降交替出现时,土体发生破坏且塌陷点偏移长江方向;周边工程建设排水施工时,降落漏斗最深处沉降量4mm,对土洞的破坏性影响最大;降雨强度不大,但降雨时间长造成的土体变形大,最大处引起土体上浮1~2mm,大大降低了盖层的强度。
参考文献
[1] Lu Y R,Liu Q,Zhang F E,et al.Environmental characteristics of karst in China and their effect on engineering[J].Carbonates and Evaporites,2013,28(1-2):251-258.
[2] 郭殿权,刘道喜.武昌陆家街地面塌陷成因机理分析[J].工程勘察,1990,(6):11-13.
[3] He K Q,Zhang S Q,Wang F,et al.The karst collapses induced by environmental changes of the groundwater and their distribution rules in North China[J].Environmental Earth Sciences,2010,61(5):1 075-1 084.
[4] 胡亚波,刘广润,肖尚德,等.一种复合型岩溶地面塌陷的形成机理——以武汉市烽火村塌陷为例[J].地质科技情报,2007,26(1):96-100.
[5] 张人权,梁 杏,靳孟贵,等.水文地质学基础[M].6版.北京,地质出版社,2011.
[6] 方 云,林 彤,谭松林.土力学[M].武汉:中国地质大学出版社,2003.
[7] He K Q,Jia Y Y,Zhao M,et al.Comprehensive analysis and quantitative evaluation of the influencing factors of karst collapse in groundwater exploitation area of Shiliquan of Zaozhuang,China[J].Environmental Earth Sciences,2012,66(8):2 531-2 541.
[8] Zhao H J,Ma F S,Guo J,et al.Regularity and formation mechanism of large-scale abrupt karst collapse in southern China in the first half of 2010[J].Natural Hazards,2012,6(3):1 037-1 054.
[9] Liu X M,Chen C X.Prediction and evaluation of space collapse for typical covered karst[J].Disaster Advances,2010,3(4):120-126.
[10] He K Q,Jia Y Y,Wang B.et al.Comprehensive fuzzy evaluation model and evaluation of the karst collapse Susceptibility in Zaozhuang Region,China[J].Natural Hazards,2013,68(2):613-629.
水动力-力学耦合模型 篇3
需水量预测是实现水资源可持续利用的一项重要而复杂的热点研究课题[1],对一个城市的规划、建设和发展具有重要的意义[2]。本研究借鉴国内外相关研究成果,采用系统动力学模型[3](SD)原理和方法,建立济宁市需水量预测系统动力学模型[4],为实现水资源供需平衡和节水管理决策提供技术依据[5]。
1 研究区需水存在的问题
山东济宁位于山东省西南部,地跨黄、淮两大流域,黄河从北部穿境而过,南四湖处于济宁腹地,京杭大运河贯穿南北。全市当地多年平均水资源量22.78亿m3,人均水资源量558 m3,不足全国平均水平的1/4,且由于水资源时空分布不均,开发利用程度不一,资源型缺水、工程型缺水、水质型缺水并存,干旱年份缺水4~7亿m3,水资源短缺已成为经济社会发展的瓶颈制约要素。将SD 模型应用于该区域,对济宁市近10年的需水量进行预测,合理配置水资源是本模型要解决的关键问题[6]。
2 SD模型建模目的
通过模型建立,达到以下目的:①以社会需水量为核心,确定各系统要素之间的关系,直观地给出“水资源一社会一经济一环境”之间的相互耦合关系;②在不同的社会经济发展模式下,动态模拟水资源供求变化,定量地描述不同利用方案下社会需水量[7];③以水资源供应量为主要依据,对不同的社会经济发展模式进行可行性分析,从中寻找出最佳方案,为研究区社会可持续发展提供决策信息[8]。
3 SD模型建模步骤
SD模型建模步骤可归纳为:①对研究区过去十年的社会、经济发展和用水情况进行调查研究,确定系统边界,并根据实际情况确定系统因果关系;②根据系统因果关系与反馈关系画系统流程图;③根据系统流程图中各个变量间的关系,利用Vensim 提供的公式编辑器建立量化的系统模拟模型,书写动力学方程;④对模型的有效性进行检验;⑤确定现状年份,进行计算机仿真;⑥设计对比方案,寻找最优方案。
4 研究区系统动力学模拟
将SD模型应用于济宁市,对济宁市未来10年社会需水量进行预测,并通过用水参数调控,实现水资源供需平衡。
4.1 系统流程图
根据系统因果分析设计出系统动力学流程图,见图1。
4.2 模型参数拟定
利用2001-2005年历史统计数据,对模型参数进行拟定。由于模型复杂,参数较多,以下仅列举出人口、工业增加值、第三产业增加值3个参数的估计途径。
通过非线性多项式拟合分析,得到以下关系式。
(1)人口预测模型关系式为:
式中:y为总人口,万人;x为年份。
(2)工业增加值预测模型关系式为:
式中:y为工业增加值,亿元;x为年份。
(3)第三产业增加值预测模型关系式为:
式中:y为第三产业增加值,亿元;x为年份。
应用上述预测模型测算得到不同规划期的人口、工业、第三产业增长率指标。
4.3 模型有效性检验
模型投入实际应用和对其结果进行分析之前,必须检验模型的模拟运行结果和实际系统行为是否相符,以保证模型的有效性和可靠性。本文系统动力学模型的检验方法:①直观检验。根据专业知识和所掌握的系统动力学建模方法,直观地对模型的变量定义、边界选择、因果关系及系统方程进行检验,并做出正确的判断;②运行检验。建立系统动力学模型,并运用Vensim PLE软件进行模拟仿真,使用软件的模型检测和单位检测功能来检验模型方程的正确性、系统参数的合理性;③历史性检验。将模拟仿真部分结果与已有的历史数据进行对比,检验其符合的程度,判别模型的可信度。
前两项要求在建模过程中已经进行了考虑,这是因为系统动力学模型的检验证实工作贯穿于循环反复建模的过程中,济宁市需水预测SD模型在建模过程中满足前2项要求。现对第三项进行详细的模型有效性检验,即进行历史检验。被检验的变量有总人口、工业增加值、三产增加值、工业用水量、农业用水量、三产用水量、生活用水量和济宁市水资源需求总量等。由于模型参数过多,受篇幅所限,本研究以人口、工业增加值、三产增加值变量为例,进行历史检验。检验时间为2006-2009年,结果见表1。
由表1可知,模拟值与实际值比较,最大误差为6.64%,模型的模拟仿真结果与实际值拟合度好,满足SD建模精度要求,证明模型具有良好的可靠性。
4.4 决策方案分析
通过调整变量的取值,组合设计出常规发展方案、加大节水与再生水利用力度方案、调整产业结构方案和综合发展方案4种对比方案。利用模型对不同对比方案下的系统发展趋势进行模拟预测,分析预测结果,优选出能保障社会经济持续发展的水资源供需平衡推荐方案。
4.4.1 方案一:常规发展模式
该发展模式是按照当地“十二五”规划指导思想、目标及布局顺延发展,本模型将该方案作为基本方案,模拟结果见表2。
对以上模拟运行结果进行分析,得到以下结论:①从用水状况来看,农业用水占主导地位,占总用水量的70%左右,而同时农业的投入产出比相对要远低于其他产业,因此毛灌定额的调控是方案制定的必要指标;②水资源的供需差额对系统发展影响很大,在常规发展条件下,水资源可供量已远远不能满足社会、经济和环境的要求。因此,应采取各种可行性措施降低水资源的消耗量,以实现水资源的供需平衡;③第三产业增加值在未来10年增长迅速,而用水量变化趋势不大,工业产值增长迅速,但工业需水量在未来10年也高速增长,由此看来,在保证国内生产总值持续增长的趋势下,提高第三产业比例,限制高耗水工业增长速度具有很大的节水潜力。
4.4.2 方案二:加大节水与再生水利用力度方案
以水资源可供水量为约束条件,调整用水定额和加大污水处理率、再生水利用率。
将方案二的模拟运行结果与基本方案相比较,由图2可见,通过采用以节约用水为核心的管理方案,济宁市未来10年内的缺水量得到很大程度的改善。但由于济宁市水资源可供给量的有限性,方案实施后,还存在一定程度的缺水量。
4.4.3 方案三:调整产业结构方案
该方案在保证国内生产总值持续增长的前提下对产业结构进行适度调整,模型运行结果见图3。可见,通过产业结构调整,鼓励发展第三产业,限制高耗水工业的发展速度,可以缓解水资源供需紧张的状况,但若仅采用该措施,节水效果并不明显。
4.4.4 方案四:综合方案
该方案决策变量的取值是在以上各方案模拟实验的基础上形成的综合方案,即综合发展方案,主要变量取值见表3。模型运行结果见表4、图4所列。
注:表中部分参数依据山东济宁市水利发展“十二五”规划及国民经济发展纲要进行选取。
通过表4、图4对比分析可见,4种方案中,方案四到2020年基本达到了水资源的供需平衡,该方案使节水量大于以上任何一个方案,缺水量几乎降到零,达到济宁市水资源的供需平衡,相对基本方案将使济宁市2020年的缺水量由3.577亿m3降到0.067亿m3,到“十二五”规划末,缺水率降低到0.38%,到2020年,缺水率降到0.23%,所以选择方案四作为济宁市需水预测研究的最优方案。该方案全面考虑了节约用水,污水治
理回用,优化用水产业结构等措施,并将各种措施定量化处理,既缓解了水资源压力又加强了污水处理力度,同时促进了整体经济水平的提高,有利于济宁市社会、经济与环境的协调发展,从而使济宁市的经济社会发展获得良好的可持续性,同时水资源也得到了可持续利用。
5 结 语
针对研究区水资源、人口、社会经济的发展现状,以水资源供需平衡为核心,建立济宁市需水量预测系统动力学模型。模拟预测济宁市2010-2020年期间社会需水量,并通过状态变量调控,设计比选方案,实现水资源供需平衡。研究表明,建立的模型是有效的,得出实现济宁市社会经济可持续发展,需要在节水和治污方面加大投入力度,同时辅助产业结构调整。本模型将涉水措施的调控力度进行了定量化处理,为节水管理决策提供决策依据。
参考文献
[1]李杰兰,陈兴鹏,王雨,等.基于系统动力学的青海省可持续发展评价[J].资源科学,2009,31(9):1 624-1 630.
[2]达娃.西藏地区需水量预测分析[J].人民长江,2009,40(1):37-38.
[3]Forrester Jay W.World Dynamics[M].Cambridge,Mass.,TheMit Press,1971.
[4]刘俊良,臧景红,何延青.系统动力学模型用于城市需水量预测[J].中国给水排水,2005,21(3):31-34.
[5]严广乐.系统动力学:政策实验室[M].上海:科学技术出版社,1997:13-16.
[6]孙宗凤.系统动力学在水资源管理中的应用[J].水利水电技术,2005,36(6):14-16.
[7]袁汝华,耿小娟,邱德华.区域水资源供需的系统动力学仿真[J].水利经济,2007,25(4):7-10.
水动力-力学耦合模型 篇4
黄酮作为陈皮中一大类主要活性成分,种类繁多,主要包括橙皮苷( 陈皮苷、橘皮苷) ,新陈皮苷、川陈皮素( 蜜橘黄素) 、柚皮苷元( 柑橘素) 、二氢川陈皮素、红橘素等［3］,现代研究表明黄酮类化合物具有抗自由基、抗氧化、抗癌、防癌、抗菌、抗炎、抗病毒、抗过敏及抗糖尿病并发症等多种生理活性及药理作用,且无毒无害,对肿瘤、心血管病等多种疾病的治疗和预防有重要意义。由此可见,陈皮中总黄酮具有很好的开发前景和利用价值,有必要对陈皮中总黄酮的提取进行研究。参考本课题组已建立的黄酮类成分的提取动力学模型［6-10］,本实验对陈皮中总黄酮的水提取动力学过程进行了研究,建立陈皮总黄酮成分水煎煮提取的动力学模型,以扩大提取动力学方程的适用范围,建立黄酮类成分普适的提取预测模型的验证,丰富数据库资源,提高预测准确性。
1 仪器与药物
陈皮药材( 生产批号: 401003645,北京同仁堂有限责任公司) ,由北京中医药大学中药鉴定系刘春生教授鉴定; 橙皮苷对照品( 生产批号: 1807 /16335,纯度> 98% ,月旭材料科技上海有限公司) ;Tu180 型紫外分光光度计( 北京普析通用仪器有限公司) ; DZKW-4 型水浴锅( 北京中兴伟业仪器有限公司) ; 甲醇为分析纯( 北京化工厂) ,水为蒸馏水。
2 方法
2. 1 陈皮总黄酮测定方法的建立
2. 1. 1 对照品溶液和供试品溶液的制备精密称取橙皮苷标准品11. 78 mg到50 m L容量瓶中,用甲醇定容到刻度,得到0. 2356 mg /m L的橙皮苷标准品溶液。取陈皮15. 00 g,加300 m L蒸馏水回流提取60 分钟,滤过,得陈皮水提液作为供试品溶液。
2. 1. 2 测定波长的选择取上述2 种溶液,取适量蒸馏水稀释后用紫外可见分光光度计在紫外区( 190 ~ 600 nm) 进行扫描,结果,供试品溶液在28 ~285 nm处均有最大吸收,而橙皮苷标准品溶液在284 nm处有最大吸收,故选定284 nm作为测定波长［5］。
2. 1. 3 标准曲线的测定精密吸取上述橙皮苷标准品溶液0. 2、0. 4、0. 6、0. 8、1. 0、1. 2、1. 5 m L分别置于10 m L量瓶中,加蒸馏水稀释至刻度,摇匀,以相应蒸馏水溶液为空白对照,在284 nm处测定其吸光度。结果表明,在0. 004712 ~ 0. 03534 mg /m L范围内,吸光度( Y) 和橙皮苷质量浓度( X) 呈良好线性关系,其回归方程为: Y = 30. 877X - 0. 0072,r =0. 9991。
2. 1. 4精密度实验取浓度为0. 2356 mg / m L的橙皮苷溶液0. 8 m L,置10 m L量瓶中,加水稀释至刻度,重复进样6 次,计算其RSD为0. 0876% ,表明精密度良好。
2. 1. 5 稳定性实验精密吸取同一供试品溶液,分别于配置后10、20、30、40、50、70、90、120、150 分钟进行含量测定,测定其吸光度并计算RSD为0. 15% ,试验证明稳定性良好。
2. 1. 6重复性实验取15. 00 g陈皮药材,加水300 m L回流提取60 分钟,滤过,取提取液100 μL至10 m L量瓶中,加蒸馏水稀释至刻度,摇匀,作为供试品溶液,测定含量,以上实验重复进行6 次,计算相对标准偏差,结果RSD为1. 29% ,试验证明重复性良好。
2. 1. 7 回收率实验精密称取陈皮药材1. 00 g,置50 m L的圆底烧瓶中,再分别加入橙皮苷标准品1. 698 mg,加20 m L蒸馏水,回流提取30 分钟,提取液放冷,取续滤液60 μL至10 m L量瓶中,平行6份,加蒸馏水稀释至刻度,摇匀,进行含量测定。以上实验重复进行6 次,结果平均回收率99. 98% ,RSD为1. 77% 。
可以看出,所建立的陈皮总黄酮测定方法的各项方法学考察符合规定,因此根据此方法对其水提过程中的陈皮总黄酮含量进行考察,建立其提取动力学模型并进行验证。
2. 2 陈皮水提取动力学模型的研究
2. 2. 1基于FICK定律得到的中药材提取通用的动力学模型基于FICK定律以及希格比理论对提取过程的微观阐释,同时也根据实际应用的需求,选取重要的影响因素如提取时间、溶剂倍量、颗粒粒度作为实验参数,简化或转化复杂的非必求参数例如传质系数、药材密度、颗粒数目,以减小方程的复杂程度。但文献中应用到的方程大都存在方程形式复杂,需要求解的参数过多,且求解困难这一问题,因此欲建立与实际生产相适应的模型,就必须解决上述问题,通过一系列假设和数学推导最终得到中药材提取通用的动力学模型［6-7］:
它表示了煎煮出的有效成分浓度与提取时间、溶剂倍量及颗粒粒度之间的函数关系。CB为提取液中目标成分的质量浓度,单位为mg/m L。M为溶剂倍量,单位为m L/g。σ对于性状不同的饮片代表不同的含义,对于葛根这种近似为正方体形的药材,其代表饮片的一个边长,单位是μm。f为方程中的一个固定参数,与药材中目标成分的含量有关。R为干药材的吸溶剂率,由于本实验使用溶剂为水,即为药材的吸水率。n值与扩散系数有关,值为负。
2.2.2理论基础由于中药材饮片多数存在形状不规则的问题,故每味药材的粒径值即σ1是不易通过实验的方式确定的(除了少数形状比较规则的中药饮片如葛根等)。另外,在实际生产中,药材多以饮片投料,粒径这一因素是固定不变的。所以有必要对方程进行简化,将σ1这个参数由需要实验直接测定转变为最终固定参数a当中的一部分。这样减少药材粒径的直接测定能够增加模型的准确性,扩大模型的适用范围,并且和实际生产情况更加接近。因此根据合理的数学推导[8],得到了简化后的模型:
两边同时取对数可以得到:
通过对不同时间和溶剂倍量作为自变量[t1 /2/( M-R) ],提取液浓度CB为因变量。分别代入上述方程中,可以求得a和n值。从而得到提取动力学方程进行预测。
2. 3 陈皮模型方程建立试验方法
2.3.1固定参数f1和R的测定根据公式(1-1)的推导过程可知,。C10为溶液沸腾时溶液中溶质的浓度,C1b为沸腾时溶质药材的浓度,(C1b-C10)应当为药材有效成分的浓度,基于本实验,该值应为本批次陈皮药材中总黄酮的含量。以2010年版《中华人民共和国药典》(一部)陈皮项下含量测定中供试品溶液的处理方法[1],制备陈皮的总黄酮测定的供试品溶液,按照上述建立的含量测定方法进行测定,得到本批次陈皮总黄酮含量,并计算f1=0.089526,
分别取药材约2. 00、4. 00、8. 00、12. 00、16. 00、20. 00 g置于500 m L圆底烧瓶中,加入蒸馏水300 m L,加热回流60 分钟,滤过,取出称重。以药材湿重与干重的差值为Y对药材干重X做线性回归,回归方程的斜率即为中药材吸收溶剂率R。平行进行3 次,求算其平均值,求得陈皮的吸水率R为2. 8401 m L/g。
2. 3. 2 提取时间对提取率的影响采用圆底三口烧瓶,电热套加热进行水煎煮提取,为减少水分蒸发,安装了回流冷凝装置。固定溶剂倍量M =12 m L / g,在提取时间点为10 分钟、20 分钟、30 分钟、40 分钟、50 分钟、60 分钟时分别取出5 m L提取液后再补充5 m L。提取液放冷,待测; 固定溶剂倍量M = 16 m L/g,在提取时间点为10 分钟、20 分钟、30 分钟、40 分钟、50 分钟、60 分钟分别取出5 m L提取液后再补充5 m L。提取液放冷,待测。
2. 3. 3 溶剂倍量对提取的影响固定提取时间为20 分钟,分别于10、12、14、16、18 和20 溶剂倍量下,取出一定体积提取液,放冷,待测; 固定提取时间为40 分钟,分别于不同10、12、14、16、18 和20 溶剂倍量下,取出一定体积提取液,放冷,待测。
3 结果
3. 1 提取时间和溶剂倍量对提取的影响结果
采用上述建立的紫外分光光度法进行含量测定。对上述所有提取数据进行整理,求得ln[t1 /2/ ( M-R) ]和ln CB进行线性回归。结果见表1,图2。
3. 2 模型的检验结果
把图2 中的回归方程和( 1-3) 式中相应数值联立,1 /( 1 - n) = 0. 5902,1 /( 1 - n) lnαf1= 1. 2506。可以求得n = - 0. 6943,α = 92. 9600,代入( 1-2) 式中求得本次研究所建立的陈皮总黄酮水煎煮提取过程动力学模型: CB= ( 8. 3223 × t0． 5/ M - 2.8401)0． 5902。根据所建立的动力学方程式,可以计算出不同提取时间、不同溶剂倍量条件下陈皮水提液中总黄酮的浓度,计算结果见表2。
以上检验数据表明,本次研究所建立的动力学模型能比较好地与试验结果相吻合,各条件下的试验数值与计算数值的标准偏差基本能满足工业标准偏差要求( 小于10% ) 。
4 讨论
本实验陈皮总黄酮的测定方法精密度、稳定性、重复性、加样回收率均符合要求,稳定可行,可以用于陈皮总黄酮的测定。推导出的提取动力学模型对陈皮总黄酮的水提过程也有很好的适用性。但在实验中发现,超过一定范围后预测误差就会偏大,实验中超过80 分钟后,浓度增长很缓慢。分析原因可能是因为陈皮果皮质地柔软,总黄酮的成分容易煎出,因此陈皮煎煮过程中溶液能较快的达到饱和。在工业生产时,综合考虑到效率、能耗,一次提取时间不长,因此本动力学模型具有一定的参考价值［11］。
由实验结果可知,陈皮中总黄酮测定,通过简化后的模型,对两边同时取对数,得到了提取液浓度与提取时间和溶剂倍量之间的线性关系。然后在此基础上求解a和n更加方便,同时可以根据所有实验数据回归得到整体拟合度更高的回归方程,把两个因素之间共同作用的关系体现出来。相比课题组之前［6-8］研究中选择的求取n、a平均值更能反映整体趋势。