非稳态模拟

2024-07-09

非稳态模拟(精选八篇)

非稳态模拟 篇1

1 温度场计算理论

根据固体传热学理论的傅立叶定律, 对无内热源系统, 瞬态温度场问题中第三类边界条件的热量平衡控制微分方程为[3]:

式中, c为混凝土比热;ρ为混凝土密度, λ为混凝土的导热系数, x, y, z为结构计算域Ω内的坐标;T为结构中鄣x, y, z鄣点处的混凝土在时刻的温度;f鄣x, y, z鄣为结构初始温度场;n为结构边界面的法线方向;Tf为随时间变化的环境温度, 本文按环境温度为恒定值考虑;α为A边界的对流换热系数;A为包围三维区域Ω的封闭曲面。

对 (1) 式的非稳态导热问题, 可用变分原理来解决。令时间变量t暂时固定, 先考虑在某具体瞬态下对泛函变分, 然后再考虑t的变化把鄣T/鄣t用差分展开。由变分原理可以证明, 在各种温度边界条件下, 结构的瞬态温度场控制方程等价于泛函:

进行有限元分析时, 将求解域离散, 任意t时刻单元内的温度分布函数取为:

式中, N=鄣N1, N2, …, Nm鄣, m为单元节点数, Ni鄣i=1, 2, …, m鄣为单元形函数, Te=鄣T1, T2, …, Tm鄣T为单元节点温度向量。将 (3) 式带入 (2) 式, 再对单元节点温度向量在单元计算域内Ω变分, 即可导出求解温度场的有限元控制方程:

式中:

节点温度Te为时间的函数, 采用向后差分来离散时间域, 有:

式中, △t为时间步长。将 (9) 式带入 (4) 式, 有

从 (10) 式可逐步计算时间离散点上的各空间离散点处的温度, 从而可得到非稳态温度场的数值解。

2 温度场计算模型

选择某核电站方案论证阶段12000m2逆流式自然通风冷却塔为计算模型, 其主要几何参数为:冷却塔总高, 170.0m;进风口高度, 11.6m;壳体底部斜率, 0.30;底部直径131.08m;喉部直径75.188m;喉部高度, 127.5m;壳体最大厚度1200mm;壳体最小厚度220mm;塔顶厚度400mm。

采用空间六面体温度单元对上述壳体结构进行网格剖分, 在厚度方向控制单元数为4个。按照上述网格剖分原则, 冷却塔壳体结构的总单元数为17496个, 节点数为23760个。冷却塔塔筒非稳态导热计算三维有限元模型见图1。

3 温度场计算条件及计算结果

3.1 计算参数

为方便与《火力发电厂水工设计规范》的稳态温度场筒壁温差值进行比较, 混凝土的热力学参数与规范建议值取为一致, 即:c=960.0J/ (kg·℃) ;ρ=2500kg/m3;

环境温度按《火力发电厂水工设计规范》附录K中, 环梁有挡水设施, 大气温度为-15℃地区, 采用母管系统的冷却塔为计算环境温度条件。

3.2 筒体三维非稳态导热计算结果

筒体在冷却塔运行1.0h、12.0h、1.0d、2.0d、3.0d、5.0d、7.0d、10.0d、15.0d、时的温度场云图见图2~图10 (各温度均为开氏温度) , 高度为20m、50m、100m、150处的塔筒外部壁面温度变化曲线见图11。

计算结果表明, 冷却塔双曲等厚部位温度场最先达到稳定导热状态, 随着时间推移, 冷却塔塔筒外部温度随壳体厚度增加而依次趋于稳定。其中, 20m高处达到稳态导热状态的时间为293.0h, 50m高处达到稳态导热状态的时间为53.0h, 双曲等厚部位达到稳定导热状态的时间为43.0h (100m高处、150m高处) 。图11中各曲线的变化趋势表明, 冷却塔塔筒厚度越大, 达到稳态导热的时间就越长, 且达到稳态导热时的壁面温度越低。

3.3 筒体一维、三维非稳态导热计算及与规范值的比较

《火力发电厂水工设计规范》、《工业循环水冷却设计规范》中关于塔筒筒壁内外表面温差的计算公式如下:

式中:α0、αi为筒壁外、内面向空气的放热系数;h为筒壁厚度 (m) ;λh为混凝土的热导率[规范单位为W/ (m2·℃) , 应为W/ (m·℃) , 建议升版或修编时改正];△tb为筒壁内外表面温度差 (℃) ;△t为筒壁内外空气温度差, △t=ti-t0 (℃) ;Kch为传热系数 (W/m2·℃) 。

鉴于《火力发电厂水工设计规范》、《工业循环水冷却设计规范》中关于冬季壁面温差的计算公式均为一维推导结果, 因此, 本文亦对不同厚度的壳体进行了一维非稳态导热计算。塔筒不同厚度、不同计算维数、不同计算方法的筒壁温差及达到稳态传导时间见表1。为了直观体现混凝土非稳态导热情况, 选取500mm厚混凝土壳体的导热计算结果为代表, 给出了一维非稳态计算的内表面温度-时间变化曲线、外表面温度-时间变化曲线、内外壁温差-时间变化曲线, 分别见图12、图13、图14。

注:上述计算的△t=ti-t0均为30.0℃, 三维计算值均来自本文3.2节。

计算结果表明:

⑴数值解和解析解完全一致, 证明数值解具有极高的精确度;

⑵一维计算和三维计算均表明, 达到稳态传热的时间与塔筒厚度正相关;

⑶一维计算结果表明, 达到最大温差的时间小于达到稳态导热的时间, 其原因是内壁面升温较快, 能迅速达到壁面最大温差, 然后内外壁面再同步缓慢升温;

⑷三维计算结果表明, 达到最大温差的时间和达到稳态导热的时间差别不大, 其原因主要是一维计算未考虑塔筒内部的热传导问题。

4 结论及建议

综合前面的计算及论述, 得出如下结论及建议:

⑴冷却塔开始运行时, 塔筒主要以非稳态对流换热方式传递热量。其传热过程为: (1) 冷却塔开始运行; (2) 热蒸汽和塔筒内壁对流换热, 塔筒内表面温度升高; (3) 随着时间推移, 温度变化波及的范围逐渐扩大; (4) 塔筒自内向外, 各截面温度依次升高, 温度变化逐层传导到塔筒的外表面; (5) 塔筒外侧表面温度升高; (6) 塔筒外表面与塔外空气对流换热; (7) 随时间推移, 壁面温差逐渐增高, 接近恒定状态, 此时, 非稳态导热过程结束, 进入稳态导热状态;

⑵《火力发电厂水工设计规范》和《工业循环水冷却设计规范》给出了冬季壁面温差计算公式, 而这一公式是基于一维稳态导热推导出来的, 不能反映塔筒的非稳态导热过程;

⑶塔筒非稳态导热数值计算表明, 达到稳态传热的时间与塔筒厚度正相关;

⑷塔筒非稳态导热一维数值计算表明, 达到最大温差的时间小于达到稳态导热的时间, 其原因是内壁面升温较快, 能迅速达到壁面最大温差, 然后内外壁面再同步缓慢升温;塔筒非稳态导热三维数值计算表明, 达到最大温差的时间和达到稳态导热的时间差别不大, 其原因主要是一维计算未考虑塔筒内部的热传导问题。

摘要:本文考虑冷却塔初始运行的实际情况, 根据瞬态温度场计算理论, 对冷却塔塔筒非稳态导热问题进行了三维数值模拟, 得出了塔筒冬季非稳态导热的时间特性。并对筒体不同壁厚情况下一维、三维非稳态导热计算及规范值计算三者进行了比较, 给出了三者结果的差异。

关键词:冷却塔,非稳态导热,温度应力

参考文献

[1]DL/T5339-2006, 火力发电厂水工设计规范[S].

[2]GB/T50102-2003工业循环水冷却设计规范[S].

四氧化二氮低温非稳态导热模型研究 篇2

四氧化二氮在低温条件下容易凝结为固体而对航天发射产生难以预见的影响,本文依据非稳态传热理论对低温环境中不同管径圆柱管道内的四氧化二氮导热情况进行了研究,建立了圆柱管道四氧化二氮低温非稳态导热模型,经验证与实际情况拟合较好.

作 者:焦天恕 信聪 齐彦兴 Jiao Tianshu Xin Cong Qi Yanxing 作者单位:焦天恕,信聪,Jiao Tianshu,Xin Cong(中国酒泉卫星发射中心,酒泉,732750)

齐彦兴,Qi Yanxing(中科院兰州化学物理研究所,兰州,730000)

非稳态模拟 篇3

关键词:储油罐,温降,FLUENT,数值模拟

一、物理模型的建立

本文利用FLUENT软件对大型储油罐的温降进行了模拟计算。图1为储油罐的二位几何模型, 对此物理模型的相关假设如下:①认为罐内横向温度温度分布均匀, 所以只对纵向二维温度场进行计算;②对计算区域进行简化, 取过油罐中心轴的径向切面为研究对象, 变三维问题为二维问题;③忽略太阳与油罐间的虑辐射换热;④忽略原油内部物理化学因素而产生的内热源, 忽略油罐边缘以及底部原油和石蜡的凝固潜热;⑤假设外界环境温度近似为随时间变化的周期性函数;⑥将储油罐浮顶近似处理为单层钢板, 同时也将罐壁、罐底假设为均匀的等厚钢板。

其中各符号的具体意义为:

To——为稳罐后储油罐内原油的初始平均温度, 为常数;

Td——为大地恒温层处土壤的温度, 为常数;

T1——为大地温度场的初始温度分布;

fw——为储油罐壁 (保温层外侧) 环境的温度函数;

ft——为储油罐顶部环境温度的分布函数;

H1——为浮顶储油罐内原油液位的高度 (10.1米) ;

H2——为储油罐罐底钢板到大地恒温层的深度 (H2=10米) ;

R——为储油罐径向平面半个平面的宽度;

δ1——为保温层厚度, 保温层平均厚度 (δ1=80毫米) ;

δ2——为储油罐罐壁平均厚度 (δ2=10毫米) ;

δ3——为储油罐顶板厚度 (δ3=5毫米) ;

δ4——储油罐底板的平均厚度 (δ4=10毫米) ;

二、数学模型

控制方程的直角坐标张量符号形式如下所示

质量守恒方程:

动量方程:

其中p为静压, τij为应力张量, gi和Fi分别为i方向上的重力体积力和外部体积力, 式中:

能量方程:

其中E为内能, keff是有效热传导系数, Jj’是组分j’的扩散流量。而且hj’为同分子物质的焓值, Sh包括了化学反应热以及其它用户定义的体积热源项。

式中: 是组分j’的质量分数。

湍动能k方程为:

湍动能耗散ε方程为:

在标准k-ε模型中, C1ε、C2ε和C3ε为经验常数;σk和σε分别是与湍动能k和耗散率ε对应的Prandtl数;根据Launder等的推荐以及后来的验证, 有下式:C1ε=1.44, C2ε=1.92, Cμ=0.09, σk=1.0, σε=1.3。

三、边界条件

1、储油罐顶部与大气直接接触, 它们之间的热交换主要形式是对流换热。公式表示为:

其中:λ——为钢的导热系数, W/ (m·℃) ;Tw——钢表面温度, ℃;

ft——大气温度函数, ℃;Γ——边界;

h——大气与钢板间的对流换热系数, W/ (m2·℃) 。

2、储油罐底部基座为大地, 取大地温度恒温层位置H2=10米, 近似认为基座左边界上为绝热。应用公式表示为:

其中:q——为热流密度q=0, W/m2;λ——大地的导热系数, W/ (m·℃) ;Γ——边界。

3、储油罐基座底边, 由罐底向下10米深处为 (H2=10米) 恒温。用公式表示为:

Td——储油罐底向下10米处至无穷远处温度, 为常数[***]为2.2℃, Γ——边界。

4、同理可知罐壁为第三类边界条件, 对称轴边界上为绝热边界条件。

四、结论

1、高温区并非集中在油罐中心, 而是相对在靠近罐壁附近的原油温度最高。

2、在自然对流过程中, 原油温降的流动分为两种形式, 一是沿着原油上液面—对称轴—罐底—侧壁面这一路进行, 这是主要形式;二是在靠近罐壁处原油沿着罐壁向下形相对一形式逆时针形态的涡流。

3、对比分析现场实测数据与应用FLUENT模拟计算的结果数据, 实验值与数值计算值的相对误差在3%以内, 符合工程上的要求。

参考文献

[1]王明吉、张勇、曹文:《立式浮顶原油储罐温降规律》, 《油田地面工程》, 2005年第24期。

[2]大庆油田有限公司储运销售分公司:《储油罐温度分布和温降规律研究评定材料之八》, 储运销售分公司档案室, 2001年。

城市基因的非稳态现象初引 篇4

关键词:稳态,城市基因,时间,记忆

“2013年4月10日, 埋有唐代高僧玄奘大师灵骨的西安兴教寺, 正面临大规模拆迁。当地政府给出的拆迁原因是“丝绸之路联合申遗的需要”, 联合申遗名单上只有兴教寺塔, 不包括寺庙内其他建筑。寺院会被作为纯粹的文物建筑而圈定, 本身及周遭全体作为景区开发, 拆迁僧寮、斋堂”。

——2013年4月12日凤凰网

1兴教寺与“城”

唐朝都城长安 , 是我国古代规模最为宏大的国都之一。其作为唐都存在的近三百年历史中, 除了城内纵横阡陌的街坊、恢弘的宫殿、精巧的寺观, 仅都城近郊就有许多风景秀丽的游览胜地。譬如都城西郊是以沣河、昆明池、定昆池等运河组成的独特水系风景区;东郊是一大片以皇家离宫和皇室贵族庄园为主体 , 以高官贵族庄园为补充的旅游风景区;南郊是主要的文化风景区, 大的有韦曲、樊川、辋川等园林景区 , 小的如辋川景区内套有的小块景区——王维的“辋川别业”; 与城南遥相对应的是隐士庄园区终南山。“城南”地区是唐代长安重要的文化景区, 诸水环绕、川源相接 , 它不但有韦、杜两姓贵族世代聚居, 也有大量的陵墓遗址 , 还是文人学士游赏吟咏的个人别业集中区。此外 , 还有大量的寺庙庄园存在 , 是宗教活动的集中地。

兴教寺, 即位于唐长安城“城南”地区——少陵塬南畔, 其南面遥望秦岭终南山, 俯视樊川, 潏河从樊川流过。唐代非常重视佛教, 在樊川之北的少陵塬边和樊川之南的中南山坡 , 有不少的寺院和古佛塔。而兴教寺更为唐代樊川八大寺之首, 且兴教寺是东土法相唯识宗祖庭之一, 是各方释子向往的禅林 (图1.1) 。就是在当下, 也是东亚各国佛教信徒游历西安时必至之地。

在唐代法相唯识宗是一支皇家宗学, 只有皇室和高级士大夫才有机会接触。但到了五代至宋元明清, 关中政治经济和文化都渐次衰落, 陕西佛教也是如此。尤其明清时代, 陕西佛教基本上呈现抱残守缺的状态。清同治年间兴教寺毁于兵火, 延至民国初, 只留塔孤立在少陵塬畔。民国二十六年, 日寇发动侵华战争, 许多老法师发表抗日声明, 呼吁人民和僧伽应积极参加抗战。在八年抗日战争期间, 兴教寺曾多次举行过护国息灾法会等佛事活动, 祈祷抗日战争胜利, 追悼阵亡将士和死难同胞。如此而塑造了百年来兴教寺背后的爱国爱教的热忱。

如今兴教寺虽远离城市, 但分布在城市外围山林或周边地带的有寺庙、道观形成的宗教文化空间, 从空间距离上看似乎远离城市的空间, 实际上和城市空间又有着不可分割的关系, 是城市景观的一部分。宗教本身, 无论佛道, 总是世人心灵的归属之地, 或是信仰皈依, 或是体闲流连。寺庙道观是与城市的人们生活紧密相连的空间场所, 与城市建立了一种融入关系 (图1.2) 。

2城市营造的内在——特质基因

穆拉托里在《城市生命与历史》中, 谈到描述城市发生和转变的规律并不是“自然”的, 而是从精确的文化行为结果中浮现出来的;现代城市营造不能抛弃其“内在”和“遗传性”知识, 还有对建筑变化内在逻辑的了解。穆拉托里所谈到的这种城市遗传性内在就是城市的历史、文化等等, 都是根植于城市的文化遗传基因。而形态上的“遗传”更为明显, 从唐长安城城市景观到现今西安城景观, 其纵横阡陌的道路形态、隐含的轴线特征依然在延续 (图2.1) , 被视为形态基因。

而基因, 本是一个生物学的概念, 它是指遗传信息的载体, 可以通过复制把遗传信息传递给下一代, 从而使后代表现出与亲代相同的性状。基因作用的表现也离不开内在环境和外在环境的影响。城市作为一个动态过程, 不断地更新、变化, 保持相对稳定的运动状态。当城市各方面的结构形态稳定下来, 相互间不可避免地发生影响, 逐渐形成该城市所特有的内在信息, 并代代相传, 形成城市文化基因、形态基因。实际上, 在城市基因的传承过程中, 一方面, 城市基因凭借其自身的秉性和位势, 不断地进行传承或传播, 保持其特有的个性;另一方面, 城市基因在传承或传播的过程中, 为了适应环境的变化, 又往往会产生一定的变异, 从而获得更好的传承或传播形式。

兴教寺所形成的宗教氛围、爱国热忱恰如一种不断传承的城市文化基因, 其文化、宗教氛围在历史长河中不断形成、演变。兴教寺是东土法相唯识宗, 两代三祖的塔院, 不论当时还是现在, 各个时期的法师及其弟子都坚持对其宗派发玄探微, 并且始终以玄奘法师为精神内核, 实践着佛教的构想。日寇侵华战争时期, 更是发表抗日声明, 谴责日寇侵华罪行, 呼吁人民和僧伽应积极参加抗战等等。

3 城市特质基因的“稳态”特性

城市的每个要素都在不断运动, 保持着可变又相对稳定的状态, 所以城市依然记录着过去, 也包容着现在种种, 事件与事件环环相扣, 连续地记载着时间。这种可变而又保持相对恒定的运动状态, 称之为“稳态”。“稳态”一词最初是运用在生物学当中, 19世纪法国著名生理学家Claud Bernard指出 : 动物越高级 , “内环境”越完备 , 摆脱外部条件制约的能力就越大 , “内环境的稳定性乃是自由和独立生命的条件”。后来 , 美国生理学家W.B.Cannon (1926) 于二十年代末期就提出了稳态 (homeostasis) 这个概念, 它是对内环境恒定概念的引伸和发展。W.B.Cannon认为稳态是指具有开放性的生命系统中协调一致的运动形态。

卡尔维诺如此描述过:“在欧林达, 你若拿着放大镜仔细寻找, 你就能在什么地方看见针头大的一个点, 稍加放大, 就能看见里面的屋顶、天线、天窗、花园和水池, 悬挂在街道上方的横幅……这个点不会停留在那里不变的, 过上一年, 它会变得有半个柠檬那么大, 然后像一只牛杆菌那样大, 然后像一只汤盘那样大;然后会成为自然大小的城市……旧城墙和旧城区一起扩展, 而新的城区也在边缘成长, 而且变细了一些, 以便给从里向外积压过来的更新地城区让位, 依次一环接一环, 直达城市核心:一座全新的欧林达。虽然缩小了尺寸, 但保持了最初的和后来所有从中衍生出来的欧林达的特征与活力……[1]”

所以不论物质的还是非物质的城市基因在长期的演化中达到历时性与共时性的统一——这就是城市基因的稳态。它是城市内在属性的延续, 是一种环环相扣的传承。从玄奘法师宣扬唯识宗, 五代至宋元明清的衰落, 民国时期参加抗战、办学堂到成立藏学院, 再到近代西安佛教复兴运动, 兴教寺始终坚持与维护宗教文化与爱国精神, 并努力传承。民国二十一年到二十三年, 由于“倡修慈恩塔院, 盖以佛法救正人心, 拨乱反治之本原也”的政治目的, 大兴土木重修兴教寺, 完成了现在中间大殿和法堂奉佛, 东边经楼奉法, 西边塔院供养祖师的稳定格局。

4 城市记忆

“城市记忆”是城市形成、变迁和发展中具有保存价值的历史记录。它是城市发展过程中长期积淀的结果, 同时当代城市人正不断地建构着新的城市记忆。城市记忆与城市基因一样, 在历时与共时中延续。城市历史痕迹的保留, 能够创造出所在场所的历史延续感和事件感。我们对记忆的推进, 不能仅仅是抽象存在, 只有保留这样的城市要素, 找到时间的载体, 并且不变其地点, 才不会泯灭对历史的回忆, 也才会持续地记录事件。

城市事件转瞬即逝, 承载过这些事件的建筑或是环境却会长期存在, 使之成为记忆历史的物质载体。然而城市的记忆不仅仅是空间上的存在, 时间轴上更应该连续记录, 记忆无论是个人的还是集体的, 都是通过叙事而形成的。只有连续的历史事件的出现更决定着场所意义的表达。由此, “城市成为了一个实在的剧场, 凝聚了事件和情感。每一次事件都包含了历史的记忆和未来的潜在记忆, 其本身也构成事件[2]”。

兴教寺, 在每一个历史层中, 都在不断叠加新的记忆, 也有不同的物质载体。兴教寺的宗教价值远远大于再创造的人文旅游价值。所以拆迁僧寮、斋堂, 毁掉的不仅仅是建筑, 更是宗教文化记忆的延续。缺乏了佛教精神的底蕴的寺院, 又怎么传递城市基因?

也正是因为城市记忆的存在, 才使城市的独特性凸显, 城市记忆构成了每座城市独特的“文化资本”。在其中, 我们分享着记忆, 同时延续着记忆, 让城市的每个片段都在时间轴上完整地被记录着。然而城市在不断拆迁过程中, 其物质载体不在、记忆不在, 城市基因也无法在历时与共时中经久、传承。

5城市“失忆”——城市基因的非稳态现象

我们的城市在不断变迁中, 差异不断消失, 每座城市的记忆也在一点一点泯灭, 被新的社会生活代替。各个时期的记忆载体不在, 城市面临失忆, 失忆的城市里也不再存在经久而稳定的城市基因, 真正表达城市精神的城市基因, 都早已找寻不到过往的踪迹, 这就是城市基因的一种非稳态现象。人们越来越被孤立, 与我们的城市分割开来, 整个城市处在一种分离的状态。我们与城市的互动越来越少, 在城市中的体验也越来越单薄。城市也愈发被现代化所引领, 并且在消费主义的引导下, 遗忘那份与城市的特殊情感。兴教寺拆迁中, 保留的“古建筑”并不能代表兴教寺所有的记忆与价值, 僧舍、兴慈楼它们都在是各个时期的兴教寺发展的“印证”, 抽离某段时期的建筑或记忆, 会磨灭其多重价值。

所以“兴教寺拆迁”无疑又是一次试验, “我们是在寻找曾经有过的辉煌时代的证据或者能找到的任何传统痕迹, 还是在判断和评估历史, 从中选出较珍贵的部分, 来留住我们眼中最好的东西?……因为他们最典型地代表了那个时代?因为它们作为群体象征的重要性?因为它们在当代具有的内在价值?因为它们作为历史信息资源而具有的特殊用途?或者我们应该 (正如我们通常所做的那样) 让随机性来代替我们的挑选, 只保留那些能够在上个世纪幸存的东西?[3]”。

6小结

一拆一建下, 城市差异性不断消失, 每座城市都与其他城市相像。城市失忆, 城市基因不再传承, 呈现出非稳态现象。就在这种城市基因的非稳态现象之下, 给城市所带来的结果就是所有的东西都冠上了“符号”来代表城市精神。就如方程式一般——固定有解。每座城市都在努力追求新的节奏中忘记了“记忆既不是短暂易散的云 雾, 也不是干爽的透明, 而是烧焦的生灵在城市表面结成的痂, 是浸透了不再流动的生命液体的海绵, 是过去、现在与未来混合而成的果酱, 把运动中的存在给钙化封存了起来[4]”。

参考文献

[1] .Matthew Carmona, Tim Health, Taner Oc, Steven Tiesdell编著, 冯江, 袁粤, 万谦, 傅娟, 张红虎译.城市设计的维度.南京:江苏科学技术出版社, 2005.

[2] .伊塔洛·卡尔维诺.看不见的城市.南京:译林出版社, 2006.

[3] .阿尔多·罗西.城市建筑学.北京:中国建筑工业出版社.2006.

[4] .南·艾琳, 张冠增译.后现代城市主义.上海:同济大学出版社, 2007.

计算机数值模拟稳态磁场重联 篇5

计算机数值模拟方法在科研工作中有着广泛而重要的应用。磁场重联是等离子体能量转化和耦合过程的一种重要的机制,是空间等离子体中非常重要的物理现象和研究内容之一。关于磁场重联的数值模拟方法一般有三种:磁流体动力学(MHD)模拟方法、全粒子模拟方法、混合模拟方法。MHD研究模型和结果在重现磁层大体形状和动力学问题方面取得了成功,但是这种模拟中没有考虑到粒子运动的影响;既包括离子动力学也包括电子动力学的全粒子模拟强调了X线处的物理特性,但是对于研究出流区的大尺度(离子尺度)结构,这种模拟方法是具有数值局限性的。[1]于是一维和二维的混合模拟方法就发展起来了。混合模拟方法将离子视为粒子,运动满足牛顿方程,而作为无质量流体的电子满足动量方程。由于混合模拟主要研究与离子运动有关的低频电磁波,所以满足准中性条件:离子数密度等于电子数密度。

Petschek[2]模型是稳态磁重联的早期理论模型之一。在这个模型中,重联结构包括三部分:一个入流区,两个出流区及一个很小的中央扩散区。这个模型描绘的是电流片两侧具有相同的等离子体密度及等大反向磁场强度的对称情况。该模型认为X线两边分别出现一对慢激波,通过这些慢激波磁能快速的被转化为等离子体的热能和动能。这个模型可以应用到地球远磁尾的磁场重联研究中。

本文采用二维三分量(2.5维)[1,3]的混合模拟方法研究Petschek稳态磁重联模型下磁场重联的结构,研究发现磁场重联后形成了突出的隆起状磁场位形,在X中性点附近的一个很小的区域内,磁场分量By出现了四极分布特征。

二、数值模型

1. 基本方程

离子的运动方程为

电子的动量方程为

电子流的宏观速度由安培定律计算

式中,为电荷耦合常数。

最后,我们用法拉第电磁感应定律推进磁场

2. 初始条件

等离子体片的中心区是电流片,初始时刻位于z=0处,其法线方向沿着z轴方向。电流片两侧分布着反向平行的磁场,反平行的磁场分量沿x方向、导向场沿y方向。

磁场沿x方向分量的初始位形为

是初始电流片的半宽度。初始时电流片两侧具有相同的导向场,假定模拟区域的导向场初值处处为零,By0(28)0。为了快速的磁场重联,我们将z方向磁场分量B z设为零值。

离子数密度初始位形为

在模拟区域的中心(7)x,z(8)(28)(7)0,0(8)处放入非均匀常数电阻率用以触发和控制重联,其空间位形为

3. 边界条件

入流边界初始时无入流扰动;出流边界取自由边界条件;磁场以无旋、无散条件外推;粒子的边界条件通过在入流边界放置的缓冲区控制。

三、模拟结果

3.1重联结构

磁场重联过程中不同时刻的磁场位形如图1所示,图中的是重联层以外区域中离子的回旋频率,本模拟中时间用无量纲化;是重联层外侧区域中离子的惯性长度,本模拟中的所有长度都用其无量纲化。从图1中我们可以看出磁场位形随时间演化的情况。时,中心电流片两侧具有反向平行的磁场分量,在计算开始的初期,磁场位形随时间的演化比较缓慢。磁场重联开始发生后,在模拟区域中心部分的磁场拓扑结构发生了明显地改变,少量磁力线变得高度弯曲。随着模拟时间的增长,磁场位形改变地十分迅速,发生重联的磁力线的范围进一步沿电流片扩张,电流片附近磁场高度变形,重联区的大小也随时间的增长而变大。在计算的后期,磁场位形变化就呈现出无明显变化的现象,说明磁场重联的发展最后会进入稳态饱和阶段,磁场位形结构相对稳定。

如图2所示给出的是重联过程中不同时刻离子数相对密度的位形图。从图2中可以看出离子数密度分布随时间演化的情况。重联开始前,随着离开电流片距离的增加,磁场强度逐渐增强,等离子体密度逐渐下降。磁场重联过程中,X中性点附近的离子数相对密度与初始时的分布相比逐渐降低,但是磁场重联使得出流区中的等离子体密度变大了。这是因为等离子体在X中性点被磁张力加速,射向出流区,所以在X中性点出现了密度分布的极小值,出流区出现了密度分布的极大值。随着重联过程的进行,由磁重联而形成的新的磁力线逐渐往出流区外移动,这个过程中等离子体会随磁力线一起被带出重联区,因此在重联的演化过程中,高密度等离子体区的范围会逐渐向出流区边界推移。

3.2 By的四极分布

如图3所示的是X中性点附近,磁场的y分量By在xz平面内的等值线图,图中的实线表示的是正值,虚线表示的是负值。我们选取的区域是X中性点附近的一个很小的区域,该区域的尺度与离子的惯性长度尺度可比,但仍远大于电子的惯性长度尺度。这就导致在该尺度范围内电子仍可视为流体,但离子的粒子性特征就不能再被忽视,离子的运动会退耦脱离磁场线,在电磁场中离子和电子的运动情况不同,因而会产生Hall电流效应。考虑Hall效应后,即使初态导向场By0=0,t>0后仍会产生四极形的By结构,这是无碰撞磁场重联的一个基本特征[4]。因此从图3中我们看到了X线周围的磁场分量By出现了明显的四极形分布结构。这一模拟结果与全粒子模拟[5]结果基本一致。

四、结论

本文采用2.5维混合模拟的计算机数值模拟方法研究了Petschek稳态磁重联模型下X中性点周围关于电流片对称的矩形区域内磁场重联的过程,主要结论如下:

1. 磁场重联使得X中性点周围的磁场位形发生了明显地改变,磁力线高度弯曲,重联过程中形成了一种突出的隆起状磁场位形。等离子体在X中性点被磁张力加速,射向出流区,X中性点附近离子数相对密度减小,而出流区离子数相对密度增加。随着磁场重联的不断进行,高密度区逐渐向计算边界推移。

2. 在X中性点附近和离子惯性长度可比的一个很小的区域内,由于Hall电流效应的作用,即使初态导向场By0=0,重联发生后X线附近仍会产生四极形的By结构。

混合模拟这种计算机数值模拟方法虽然无法了解电子的动力学行为,但它可以部分地获得离子动力学的信息及其对重联发生处的等离子体的影响,并且与完全粒子模拟相比具有计算量要小得多的优势,因此在磁重联相关问题特别是离子动力学起主导作用的问题研究中具有广阔的前景。

参考文献

[1]王旭丹.磁尾磁场重联的二维混合模拟研究.大连理工大学硕士论文[D].2006.

[2]H.E.Petschek.Magnetic field annihilation[A].National Aeronautics and Space Administration,NASA Spec.Publ.SP-50[C],Washington:1964:425-439.

[3]Lin Y,Swift D W.A two-dimensional hybrid simulation of the magnetotail reconnection layer.Journal of Geophysical Research[J].1996,101(A9):19859-19870.

[4]周国成,等.地球近磁尾无碰撞磁重联事件.空间科学学报[J].2003,23(1):25-33.

线接触非稳态弹流问题的数值计算 篇6

非稳态弹流润滑状态是工程问题中的常见问题, 如凸轮与挺柱副、齿轮与齿、轮柴油机的活塞与汽缸套等线接触问题。随着人们对非稳态载荷认识的加深, 研究非稳态润滑已成为一种趋势[1-4]。Wijnant等[5]通过实验和理论研究揭示了结构振动对点接触弹流润滑膜厚度的影响;王静等[6]用多重网格技术研究了脉冲载荷对弹流润滑的影响;严珩志等[7]利用多重网格技术探讨了谐波动载频率和幅值对线接触弹流润滑的影响规律;文献[8-10]探讨了非稳态润滑的数值求解方法和建模。由于非稳态润滑方程的求解较为困难, 常用的多重网格技术编程难度较大, 因此本文研究利用间接迭代法对余弦载荷作用下的弹流润滑进行求解, 分析其油膜及压力分布规律。这种方法随着当前计算机性能的提高显得简单可行, 易于研究者上手应用。

1数学模型及求解

1.1基本方程

设ρ为润滑油的密度;h为油膜厚度;p为油膜压力;η为润滑油的黏度;v为卷吸速度;t为时间参数;x为空间坐标;x1为入口边界坐标;x0为出口边界坐标。h0 (t) 为待定油膜厚度常数, R (t) 为曲率半径;E′ 为综合弹性模量;s为积分变量。a为润滑油的黏压系数;η0为大气压下的黏度。ρ0为大气压下的密度。

为了便于计算分析, 引入下列量纲一参数:速度参数Dsp=v0η0/ (E′R0) ;载荷参数Dlp=w0/ (E′R0) ;卷吸参数Sk=3 Ksp (t) Dsp (π/Dlp) 2;挤压参数Sv=3Dsp[ (π/Dlp) 3/2]0.5;膜厚参数h=hR0/b2;压力参数p=p/p0;黏度参数η=η/η0;密度参数ρ=ρ/ρ0;时间参数T=v0t/R0;坐标参数x=x/b;入口坐标参数xa=x1/b;出口坐标参数xb=x0/b。其中, b为赫兹接触半宽;pH为最大赫兹压力;v0为初始速度;w0为单位长度上的载荷;R0为初始的曲率半径;Ksp (t) =v (t) /v0为速度系数;Kr (t) =R (t) /R0为曲率半径系数;

将控制方程量纲一化, 得到量纲一Reynolds方程:

离散的量纲一膜厚方程:

式中, D (i, j) 为量纲一弹性变形离散系数。密度和压力关系:

黏度和压力关系:

式中, z为Roelands积分常数。

载荷方程:

1.2数值求解

采用间接迭代法求解非稳态润滑问题。在Fortran Powerstation 6.0平台上进行编程求解。 Reynolds方程中时间项的离散采用向后的差分格式, 沿时间轴将每个周期分解为240个瞬时点, 每计算完一个周期, 将结果与前一周期的值比较, 直到两者一致。 迭代初值为Hertz压力, 选用Gauss-Seidel迭代格式和双岐子迭代格式。

2计算结果分析

为了全面考虑线接触条件下多种振动冲击工况对润滑状态的影响, 对周期载荷中的总载荷的大小、变载占总载的比例及周期频率分别进行了系列的计算与分析, 得到其规律。计算的基本参数如下:E′=2283MPa;η0=0.08Pa·s;圆粒半径R =0.02m;黏压系数α=2.19×108m2/N。载荷采用余弦曲线W =W0+Wd (cos2πft) , 4个瞬时观察点A1~ A4如图1所示。

2.1载荷的影响

采用上述参数, 频率为1500Hz, 保持动载幅值与总载荷的比例为0.5不变, 逐步增大总载荷, 当量纲一总载荷分别为3.0×10-5、5.0×10-5和8.0×10-5时, 获得的计算结果如图2~ 图4所示。

比较点A1点和A3点的膜厚曲线可知, 曲线在入口区和中心区的走向呈相反趋势, 这同载荷在这两点的走向相反的特点相一致。A2点和A4点的曲线亦呈现此特点。瞬时点A1的载荷最大, 为 (W0+Wd) , A3点最小, A2点、A4点载荷相同, 这与压力分布图上的压力宽度相符。而出口处的颈缩现象随负载的加大而减小。

2.2动载比例的影响

采用上述参数, 频率为1500Hz, 保持量纲一总载荷为8.0×10-5不变, 逐步增大动载比例, 当动载幅值分别为1.0×10-5和3.0×10-5时, 获得的计算结果如图5、图6所示。

经比较可以看出, Wd幅值小时, 对各点压力和膜厚分布影响不大。Wd幅值大时, 各点压力曲线可能出现双峰, 出口处的颈缩现象减弱或消失。

2.3频率的影响

采用上述参数, 保持量纲一总载荷为8.0× 10-5、动载为2.0×10-5不变, 逐步增大频率, 当频率分别为500Hz和2500Hz时, 获得的计算结果如图7、图8所示。

经比较可以看出, 频率增大动态效应明显。 动态效应主要发生在入口区和中心区, 出口区 (颈缩区) 和最小油膜厚度变化不大;随着频率的增大, 压力分布曲线由抛物线变为双峰, 压力最大值增加。

3结论

(1) 随载荷的增大, 出口处颈缩现象减弱, 挤压效应明显。

(2) 周期性载荷的频率和幅值对油膜压力和膜厚分布有较大影响, 频率越高, 幅值越大, 挤压效应越明显。

(3) 动载幅值和频率的表现特征有所不同。 动载幅值只改变油膜的中心区, 随着动载幅值的增大, 材料表面在中心区发生局部凹陷的程度增强;频率对入口区和中心区的影响较大, 频率达到一定程度, 油膜形状呈波浪状。

摘要:非稳态弹流润滑是工程中常见状态, 为了研究非稳态润滑时设备的工作性能, 采用间接迭代法求解了周期载荷线接触弹流问题。通过算例研究工程上常见的余弦载荷作用下弹性流体动力润滑特性, 探讨载荷、动载幅值和频率变化对油膜的影响。计算结果表明:随着载荷的增大, 出口的颈缩现象减弱;随着频率和动载幅值的增大, 压力曲线的形状发生改变, 挤压效应明显。

关键词:线接触,弹流润滑,非稳态,数值计算

参考文献

[1]任志强, 张彬彬, 王静.冲击载荷条件下线接触往复运动的弹流润滑特性[J].润滑与密封, 2013, 38 (1) :28-34Ren Zhiqiang, Zhang Binbin, Wang Jing.Elasto Hydrodynamic Lubrication in Line Contact Reciprocating Motion under Impulse Load[J].Lubrication Engineering, 2013, 38 (1) :28-34.

[2]刘德良, 孙昂, 徐久军.载荷对线接触热弹流温度场的影响[J].机床与液压, 2012, 40 (1) :134-136.Liu Deliang, Sun Ang, Xu Jiujun.The Effects of Loads on Temperature Field of the Line Contact Thermal EHL[J].Machine Tool and Hydraulics, 2012, 40 (1) :134-136.

[3]梅甫良, 龙连春.线接触非稳态弹流问题的数值解[J].石油大学学报, 1996, 20 (5) :65-68.Mei Fuliang, Long Lianchun.Numerical Solution for the Isothermal Elasto-Hydrodynamic Problems under Line Contacts[J].Journal of the University of Petroleum, 1996, 20 (5) :65-68.

[4]艾晓岚, 俞海清.线接触非稳态弹流润滑完全数值解及其实际应用[C]//第四届全国摩擦学学术会议论文集.北京, 1987:21-27.

[5]Wijnant Y H, Venner C H, Larsson R, et al.Effect of Structure Vibrations on the Film Thickness in an EHL Circular Contact[J].ASME Journal of Tribology, 1999, 121 (4) :259-264.

[6]王静, 杨沛然.时变等温线接触弹流润滑问题球解的多重网格技术[J].润滑与密封, 2000 (3) :9-11.Wang Jing, Yang Peiran.Multileve Technique for the Solution to Time-dependet Isothermsl EHL Line Contact Problem[J].Lubrication Engineering, 2000 (3) :9-11.

[7]严珩志, 王红志, 钟掘.谐波振动下线接触弹流润滑的仿真及分析[J].摩擦学学报, 2002, 22 (1) :62-65.Yan Hengzhi, Wang Hongzhi, Zhong Jue.Simulation and Analysis of Elastohydrodynamic Ubrication Line Contact Problem under Sine Wave Vibration[J].Tribology, 2002, 22 (1) :62-65.

[8]梅雪松, 谢友柏.一种快速求解非稳态弹流问题的直接迭代算法) [J].应用力学学报, 1991, 8 (3) :20-28.Mei X S, Xie Y B.Direct Overlapping Method for Quick Solution of Unstable Hydrodynamic Problem[J].Journal of Applied Mechanics, 1991, 8 (3) :20-28.

[9]Ilya I Kudish.On Formulation of a Non-steady Lubrication Problem for a Non-conformal Contact[J].Tribology Transanctions, 1999, 42 (1) :53-57.

页岩储层非稳态渗透率测试方法综述 篇7

页岩气是指储集在富含有机质泥页岩中的非常规天然气,一部分以游离或溶解态赋存于孔隙和裂缝中,一部分吸附在有机质和黏土矿物的内表面,可以是生物成因、热解成因或混合成因,在一定地质条件下聚集成藏并达到经济开采价值[1,2,3,4,5]。研究表明[6,7],全球页岩气资源量约为456×1012 m3。美国、加拿大已在多个盆地进行了页岩气的商业性开采。我国主要盆地和地区的页岩气资源量非常丰富[8,9,10,11,12],约为(15 ~30)×1012 m3,勘探开发前景广阔。页岩储层物性差[13,14],孔隙喉道已达到纳米级,孔隙度约为4%~6%,基质渗透率小于0.001 mD。目前国内在页岩气成藏机理、资源潜力等方面取得了巨大的进步[15,16,17,18,19,20,21,22,23,24,25],但在页岩储层渗透率测试等方面研究相对较少。渗透率是页岩气藏储层评价、产能计算及开发方案制定所必需的关键参数之一,合理准确的渗透率测试方法对于油田开发方案的编制、经济分析具有至关重要的指导作用。

室内测试储层渗透率的方法可分为稳态法和非稳态法,稳态渗透率测试方法在常规油气储层中的应用已有50余年的历史。然而,对于煤层和页岩等致密储层,传统稳态渗透率测试方法效率低、实验过程易受环境温度影响、流速计量误差偏大[26,27,28]。基于稳态法测试致密储层渗透率的上述缺点,非稳态渗透率测试方法在页岩等致密储层中得到了较好的应用。为了研究页岩储层渗透率测试方法,对国内外致密储层非稳态渗透率测试方法相关文献进行了广泛深入地调研,详细阐述了岩芯柱脉冲衰减法、岩屑脉冲衰减法、脱气法三种脉冲衰减渗透率测试方法的实验原理、数学处理方法、影响因素、适用条件和存在的问题。针对压汞法,给出了基于不同理论的预测模型,为页岩储层渗透率测试的深入研究提供参考。

1 脉冲衰减法

脉冲衰减法[29,30,31]是基于一维非稳态渗流理论,通过测试岩样一维非稳态渗流过程中孔隙压力随时间的衰减数据并结合相应的数学模型和测试仪器所限定的初始条件和边界条件,对渗流方程进行精确解答和合适的误差控制简化,从而获取储层的渗透率参数。与常规稳态渗透率测试方法相比,脉冲渗透率测试方法不需要记录岩样的出口流速,具有较高的测试效率和精度,其在页岩等致密储层中的应用越来越广泛。脉冲衰减渗透率测试方法主要划分为岩芯柱脉冲衰减法、岩屑脉冲衰减法和脱气法[32,33,34,35,36,37,38]。

1.1 岩芯柱脉冲衰减法

1968年,Brace[30]首先提出了岩芯柱脉冲衰减渗透率测试方法,该方法在核废料存储、致密砂岩气藏储层评价、稠油提高采收率中应用较为广泛。随着页岩气藏的商业性开发,岩芯柱脉冲衰减渗透率测试方法在页岩储层评价中也得到了广泛的应用,许多学者对该方法进行了探讨[39,40,41,42,43,44,45,46,47,48]。

岩芯柱脉冲衰减渗透率测试实验中,对样品施加一定围压并在样品两端连接容器,使得实验装置中的上游腔室、下游腔室和岩样孔隙压力相等且达到平衡状态。对上游容器内气体增压待系统稳定后,打开阀门并连续监测样品两端容器内压力数据。最后,利用解析、数值等数学方法结合相关参数对压力数据进行计算分析,获取样品的渗透率。图1给出了脉冲衰减渗透率测试装置示意图。

根据质量守恒原理及达西定律,在不考虑重力作用时,气体在柱塞岩样中流动的一维非稳定渗流数学模型为:

Ρ(x,t)x=cμϕΚΡ(x,t)t(0<x<L,t>0)(1)

初始条件:P(x,0)=Pu(0)=Pd(0)

(0<x<L) (2)

边界条件:P(0,t)=Pu(t) (t≥0) (3)

P(L,t)=Pd(t) (t≥0) (4)

dΡudt=kcμϕLVsVuΡx|x=0(t0)(5)

dΡddt=kcμϕLVsVdΡx|x=L(t0)(6)

式(1)~式(6)中,P(x, t)为岩样孔隙压力,Pa;c为岩样孔隙流体压缩系数,Pa-1;μ为流体黏度,Pa·s;ϕ为岩样孔隙度,无量纲单位;k为岩样渗透率,m2;L无岩样长度,m;t为脉冲衰减时间,s;Vu和Vd为上游和下游腔室体积,m3;Pu和Pd为上游和下游腔室的气体压力,Pa。

数学处理过程中,假设岩样为均质,结合初始条件和边界条件,可利用解析方法和数值方法对方程(1)进行求解,进而获得岩样的渗透率。由岩芯柱脉冲衰减渗透率测试方法的数学模型和实验装置示意图可知:气体在岩样内部的压缩存储效应、实验装置的密封性、气体在岩样中的吸附存储及实验过程中压力监测精度均影响实验结果。在实验影响因素分析中,Trimmer,Lin [49,50]使用不同实验配置分析了气体在样品孔隙体积中的压缩存储效应对实验结果的影响。实验结果表明,在某些情况下,气体在样品孔隙体积中的压缩存储效应对实验结果影响明显。Bourbie,Chen,Amaefule,Maini,Kwan等人[51,52,53,54,55]对不同上下游容器体积组合进行了分析。Hsieh[56]基于无因次实验变量和标准曲线,提出了考虑气体在样品孔隙体积中压缩存储效应的解析求解方法。该方法能够对不同上下游容器体积组合实验结果进行求解,解除了对实验中容器尺寸的设计限制。Ning[57,58]在非稳态压力脉冲衰减实验中,分析了气体滤失对实验结果的影响,并给出了考虑气体滤失的实验结果处理方法。Luffel分析了在岩样中吸附作用对实验结果的影响,给出了考虑气体吸附作用的数学处理方法。

岩芯柱脉冲衰减渗透率测试过程中,如果页岩样品存在天然裂缝,则实验测试的渗透率为样品基质和天然裂缝渗透率的合成,难以区分基质和天然裂缝的性质。1990年Kamath[59]指出了压力脉冲衰减实验中均质岩芯和含有微裂缝岩芯的压力变化特征的差别,建立了相应的有限差分计算模型。并利用模型计算结果与实际实验测试结果进行拟合来获得测试样品基质和裂缝的物性参数。Ning将含微裂缝岩芯简化为水平平行六面体,中间存在一条裂缝,给出了压力脉冲衰减方法测试渗透率的解析解。利用该方法测试基质和裂缝渗透率。后来又给出该方法岩芯中含有多条天然裂缝和一条天然裂缝在岩芯一侧时对应的数学模型,将数值计算结果和解析方法求得结果进行对比分析。

除此之外,Holder[60]利用脉冲衰减方法测试渗透率时,假设岩芯为非均质,给出了相应的处理方法。Kamath等人利用非稳态压力脉冲衰减法测试岩芯的孔隙度和渗透率,同时应用实验过程中的压力变化数据定量求解岩芯的非均质性。基于室内岩芯柱脉冲衰减渗透率测试方法,Jochen[61]提出了一种矿场测试页岩储层单层性质(渗透率,表皮因子,储层压力)的方法。利用封隔器封隔目的层(套管或裸眼完井均适用),向井筒注入氮气并保持压力低于裂缝开启压力。关井监测压力下降过程数据,通过分析压力数据获取储层渗透率等属性参数。该测试方法可视为氮气柱脉冲衰减测试方法。

室内岩芯柱脉冲衰减渗透率测试方法具备较高的测试效率和精度,该方法能够测试给定围压条件下的页岩样品的渗透率,为储层原始应力条件下的渗透率测试和页岩储层应力敏感性研究提供了有效途径。岩芯柱脉冲衰减渗透率测试方法对岩样制备、实验装置密封性和压力监测仪器精度要求较高,目前数学处理上难以区分天然裂缝和基质的作用。因此,岩芯柱脉冲衰减渗透率测试方法适用于实验室内对页岩储层进行高精度的渗透率测试。

1.2 岩屑脉冲衰减法

页岩储层脆性较强,标准岩芯柱制备过程存在较多问题。因此,Luffel提出了岩屑脉冲衰减渗透率法测试岩芯切片或钻井岩屑原始含水饱和度条件下的渗透率。岩屑脉冲衰减渗透率测试方法将岩屑容器和装有实验气体(通常为氦气)的标准容器对接,通过标准容器内气体形成压力脉冲,连续记录岩屑容器和标准容器的压力变化数据。最后,利用相应的数学处理方法获得岩样的渗透率。

图2给出了岩屑脉冲衰减渗透率测试装置示意图,主要包括两个容器,一系列控制阀门和高精度压力传感器。实验过程中,首先打开阀门3对岩屑样品容器抽真空,待系统稳定后测试压力并关闭阀门3;打开阀门1向标准容器冲入高压气体,稳定后测试压力并关闭阀门1;打开阀门2将岩屑样品容器与标准容器连通,利用高精度压力传感器连续测试两容器的压力变化数据直至系统稳定;系统稳定后测试整个系统的平衡压力。该方法的数学处理方程如下:

ϕ=[Vr(Ρr0zr0-Ρeze)+(Vs-Vb)(Ρeze-Ρ0z0)](Ρeze-Ρ0z0)Vb(7)

k=Ra2(ϕ+(1-ϕ)Κa)μcgs1α12(8)

Κa=qρ(9)

FR=1-(Κc+1)(ρc0-ρ)ρc0-ρ0(10)

Κc=ρbVcΜ(ϕ+(1-ϕ)Κa)(11)

Vc=Vr+Vs-Vb (12)

ρc0=ρr0Vr+ρ0(Vs-Vb)Vc(13)

tanα=3α3+Κcα2(14)

式(7)—式(14)中,Vr为标准容器容积,m3;Vs为岩屑样品容器容积,m3;Vb为样品总体积(包含孔隙体积),m3;P0为岩屑样品容器抽真空后系统平衡压力,Pa;z0为岩屑样品容器抽真空后气体压缩因子,无量纲单位;Pr0为标准容器充气平衡后的压力,Pa;zr0为标准容器充气后气体压缩因子,无量纲单位;Pe为整个系统平衡后压力,Pa;ze为整个系统平衡后气体压缩因子,无量纲单位;Ra为岩屑或岩芯切片的半径,m;q为吸附气密度,m3·g-1;ρ为气体密度,g·m-3;cg为气体压缩系数,Pa-1;s1为ln(FR)随时间变化曲线中后期直线段的斜率,无量纲单位;FR剩余气体与样品排开气体比值,无量纲单位;Kc为样品容器剩余体积与样品孔隙体积的比值,无量纲单位;ρc0为标准容器和样品容器内气体平均密度,g·m-3;ρ0为标准容器初始气体密度,g·m-3;ρ为气体密度,g·m-3;ρb为样品体积密度,g·m-3;Vc为标准容器和样品容器总孔隙空间,m3;M为样品总重量,g;α1为超越方程式(14)的根。

与岩芯柱脉冲衰减渗透率测试方法相比,岩屑脉冲衰减渗透率测试方法测试效率高,测试条件不受限制,可在现场进行测试;能够测试不同含水饱和度条件下的气体渗透率,为页岩储层气水相对渗透率研究提供了一种有效测试途径;样品外形不受限制,能够测试岩芯切片、钻井岩屑等不规则样品;缩小了样品测试尺寸,理论上能够更容易去除天然裂缝的影响,为页岩储层基质渗透率的测试研究提供参考。岩屑脉冲衰减渗透率测试方法的缺点是样品不能施加围压,测试精度相对较低,测试样品形状对测试结果有一定的影响。目前在实验结果处理上将测试样品简化为规则圆柱状外形,对实验结果有一定的影响。岩屑脉冲衰减渗透率测试方法适用于现场快速岩样渗透率的测试,也可用于室内精度要求不高的渗透率测试。

1.3 脱气测试法

脱气渗透率测试方法是对岩芯柱脱气过程进行监测,通过压力和产量数据求取样品的渗透率[30,62,63,64]。室内试验测试过程中,在一定压力条件下向装有样品的容器中注入气体(通常为甲烷或氦气)。样品一端与给定容积的真空容器相连。系统平衡后打开岩样一端,气体从岩样打开端溢出,并连续测量压力变化。通过压力和气体产量数据计算样品渗透率,数学方程如下:

G=Μ2AΔΡt12(15)k=zRΤμG2ΡS(16)

式(15)—式(16)中,G为脱气量,g·cm-1·Pa-1·s-0.5;M为产气量,g;A为脱气面积,m2;P为平均压力,Pa;t为时间,s;k为渗透率,m2;z为气体压缩因子,无量纲单位;R为气体常数,J·g-1·K-1;T为绝对温度,K;μ为气体黏度,Pa·s;S为斜率,g·m-3·Pa-1。

脱气渗透率测试方法主要适用于现场密闭取芯的测试,在测量页岩储层含气量等参数的同时定量给出渗透率的大小。与上述渗透率测试方法相比,脱气渗透率测试方法测试精度较低。除上述数学处理方法外,Cui[38]还对脱气渗透率测试方法的数学处理方法进行了研究,并给出了不同样品尺寸对应的渗透率求解方法。

2 压汞法

压汞法测试渗透率是指利用压汞毛管压力数据建立预测模型对渗透率进行预测。1921年,Washburn[65]提出了压汞毛管力方程,描述了毛管压力,界面张力和接触角的关系:

R=2(0.417)σΗg-Aircos(θΗg-Air)1Ρc(17)

式(17)中,R为多孔介质孔喉半径,μm;σHg-Air为汞和空气的界面张力(通常为480 dynes·cm-1),dynes·cm-1;θHg-Air为汞和空气的接触角(通常为140度),度;Pc为毛管压力,psia。

压汞毛管力方程是压汞渗透率预测模型的基础方程。利用压汞毛管压力曲线预测渗透率主要有两种基础理论[66,67,68],渗流(特征尺度)理论和Poiseuille理论。渗流(特征尺度)理论认为流体在多孔介质中的流动过程受一个或几个特征尺度控制。Poiseuille理论将多孔介质处理为具有不同孔隙半径的微观模型。两种理论的区别在于Poiseuille理论认为流体在多孔介质中的流动路径是确定性分布的,而渗流(特征尺度)理论认为流体在多孔介质中的流动路径是随机的。

2.1 基于Poiseuille理论压汞渗透率预测模型

2.1.1 Purcell模型

Purcell,Rose等人[69,70]以Poiseuille理论为基础,通过引入岩性因子描述多孔介质的迂曲度,提出了压汞渗透率预测模型:

kΡurcell=14254FΡurcellϕ0100d(SΗg)Ρc2(18)

Purcell利用压汞渗透率预测模型对27块Upper Wilcox和Paluxy砂岩储层岩样进行了分析,预测了无围压条件下样品的空气渗透率。

2.1.2 Huet-Blasingame模型

Huet[71]以Purcell,Burdin研究成果为基础,结合Brooks[72]和Nakornthap[73]函数拟合参数,用指数函数拟合毛管压力曲线,对89块砂岩样品进行了分析,样品渗透率分布范围为(0.004 1-834 0)×10-3μm2,孔隙度分布范围0.003—0.340。Huet压汞渗透率预测模型为:

Ρc=Ρd[Sw-Swi100-Swi]-1λ(19)

kΗ-B=81718.86691(Ρd)1.7846[λλ+2]1.6575(100-Swi)0.5475ϕ1.6498(20)

式(18)—式(20)中:kPurcell为Purcell渗透率,mD;kH-B为Huet-Blasingame模型渗透率,mD;FPurcell—Purcell岩性因子(与岩芯迂曲度有关,岩性因子随岩石类型变化而变化),无量纲单位;Pd为毛管压力,psia;Swi为模型参数,无量纲单位;λ为Brooks-Corey 孔隙分布指数,无量纲单位。

2.2基于渗流(特征尺度)理论压汞渗透率预测模型

2.2.1 Winland模型

Winland[74]以渗流(特征尺度)理论为基础对56块砂岩和26块碳酸盐岩样品进行了分析,建立了进汞饱和度35%对应的孔隙半径,空气渗透率,孔隙度之间的预测模型:

lg(R35)=0.996+0.588lg(kWinland)-0.864lg(ϕ) (21)

kWinland=49.4R351.7ϕ1.47 (22)

2.2.2 Swanson模型

Swanson[75]认为,在压汞毛管压力曲线结果中,绘制进汞总体积和毛管压力的比值曲线,且总体积和毛管压力比值的顶点代表了所有主要控制渗透率的孔隙空间均已经被汞充满。Swanson对319块砂岩和碳酸盐岩样品进行了分析,建立了空气渗透率与进汞体积和毛管压力的关系式:

kSwanson-air=399[SbΡc]A1.691(23)

同时,他还研究了58块样品在1 000 psi围压条件下盐水测试渗透率与进汞体积和毛管压力的关系,提出了盐水渗透率模型:

kSwanson-brine=431[SbΡc]A2.109(24)

2.2.3 Walls-Amaefule模型

Walls[76]以Swanson公式为基础对渗透率低于0.01 mD致密砂岩气储层进行了研究。尽管使用了不同的方法确定(Sb/Pc)A,但预测结果与Swanson公式预测结果一致。

kWalls-Amaefule=30.5[SbΡc]A,Ηg-air1.56(25)

2.2.4 Katz-Thompson模型

Katz,Thompson[77,78,79]指出流体在多孔介质中的流动过程主要受三个特征尺度控制:特征长度,压汞试验中,汞最先进入的孔隙对应的直径,对应的压力称为“门限压力”,是压汞毛管压力曲线的初始点;最大水力长度,最高水力电导率对应的有效孔隙直径;最大导电率长度,最大离子电导率对应的有效孔隙直径。同时,还给出了压汞毛管压力曲线预测渗透率的Katz-Thompson模型和不同导电率比值和不同特征长度的对应数值:

kLc=[1013226]Lc2[σoσw](26)

[σoσw]=ϕm(27)

[σσo]=[LEmaxLc]ϕSLEmax(28)

kLE=[1013226]Lc2[LEmaxLc]ϕSLEmax(29)

kLΗ=[101389]LΗmax2[LΗmaxLc]ϕSLΗmax(30)

2.2.5 Kamath模型

Kamath[80]对毛管压力数据预测渗透率模型进行了回顾,并使用进汞饱和度与毛管压力比值的最大值确定渗透率(致密气藏储层渗透率小于1.0 mD):

kΚamath=413[SbΡc]A,Ηg-air1.85(31)

2.2.6 Pittman模型

Pitrman[81]对202块来自14个不同储层的样品进行了分析,渗透率分布范围(0.05—100)×10-3 μm2,孔隙度分布范围0.03—0.28。在Winland模型基础上,利用门限压力对应了孔喉半径取代进汞饱和度35%时对应的孔喉半径,建立了压汞渗透率预测模型:

kΡittman=32.3Rapex1.185ϕ1.627(32)

2.2.7 OU模型

Dastidar[82]对150块砂岩样品进行了测试分析,样品渗透率分布范围(0.000 1-10 000)×10-3 μm2,孔隙度分布范围0.005—0.300。根据孔隙半径分布几何平均值提出特征尺寸,建立压汞渗透率预测模型:

Rwgm=exp[i=1nwiln(Ri)i=1nwi](33)wi=αiαΤ(34)kΟU=4073Rwgm1.64ϕ3.06(35)

式(21)—式(35)中,kwinland为Winland模型预测渗透率,mD;kSwanson-air为Swanson模型预测气测渗透率,mD;kSwanson-brine为Swanson模型预测盐水渗透率,mD;kWalls-Amaefule为Walls-Amaefule模型预测渗透率,mD;kKamath为Kamath模型渗透率,mD;kPittman为Pittman模型渗透率,mD;R35为进汞饱和度35%对应的孔喉半径,μm;[Sb/Pc]A为进汞体积与毛管压力比值曲线的顶点;Lc—特征长度,μm;m—胶结指数,无量纲单位;σo—饱和水岩芯传导率,S·m-1;σw—饱和岩芯所用水的传导率,S·m-1; Ri—第i个毛管压力对应的孔隙半径;wi—第i个毛管压力值对应的孔隙半径Ri的权重;n—毛管压力测试总数据点个数。

压汞法预测渗透率主要受样品类型影响,不同样品类型渗透率预测结果差异较大。然而,压汞法可以作为室内渗透率测试的辅助工具预测渗透率趋势。

3 结论

(1)对国内外致密储层非稳态脉冲衰减渗透率测试方法和压汞法进行了广泛深入地调研,详细阐述了岩芯柱脉冲衰减法、岩屑脉冲衰减法、脱气法三种脉冲衰减渗透率测试方法的实验原理、数学处理方法、影响因素、适用条件和存在的问题。针对压汞渗透率测试方法,给出了基于不同理论的预测模型,研究结果为页岩储层非稳态渗透率测试方法的深入研究提供参考。

(2)岩芯柱脉冲衰减渗透率测试方法具备较高的测试效率和精度,能够测试不同围压条件下的储层渗透率,为页岩储层原始条件下的渗透率测试和应力敏感性研究提供了有效的实验手段。然而,岩芯柱脉冲衰减渗透率测试方法对实验仪器和样品制备要求严格,适用于室内对页岩等致密储层进行精密的测试。数学处理上还需要进一步的研究,以区分基质和天然裂缝的作用。

(3)与岩芯柱脉冲衰减法相比,岩屑脉冲衰减法测试条件不受限制,可在现场进行测试;能够测试不同含水饱和度条件下的气体渗透率,为页岩储层气水相对渗透率研究提供了一种有效测试手段;样品外形不受限制,能够测试岩芯切片、钻井岩屑等不规则样品;缩小了样品测试尺寸,理论上能够更容易排出天然裂缝的影响,为页岩储层基质渗透率的测试研究提供参考。岩屑脉冲衰减方法适用于现场快速岩样渗透率的测试,也可用于室内精度要求不高的渗透率测试。目前在实验结果处理上将测试样品简化为规则圆柱状外形,对实验结果有一定的影响。

(4)脱气渗透率测试方法主要适用于现场密壁取芯的测试,在测量页岩储层含气量等参数的同时定量给出渗透率的大小,可作为渗透率测试的参考数据。

(5)压汞法测试致密储层渗透率受样品类型影响较大,不同类型样品渗透率预测差异较大,测试结果同样可作为渗透率测试的参考数据。

摘要:为了研究页岩储层非稳态渗透率测试方法,对国内外脉冲衰减法和压汞法相关文献进行了广泛深入地调研。详细阐述了岩芯柱脉冲衰减法、岩屑脉冲衰减法、脱气法三种脉冲衰减渗透率测试方法的实验原理、数学处理方法、影响因素、适用条件和存在的问题。针对压汞法,给出了基于不同理论的预测模型。研究结果表明:岩芯柱脉冲衰减法具备较高的测试效率和精度,能够测试不同围压条件下的储层渗透率,为页岩储层原始条件下的渗透率测试和应力敏感性研究提供了有效手段。岩芯柱脉冲衰减法对实验仪器和样品制备要求严格,适用于室内对页岩等致密储层进行精密的渗透率测试,且数学处理上还需要进一步的研究,以区分基质和天然裂缝的作用。岩屑脉冲衰减法不受样品形状限制,缩小了样品测试尺寸,理论上能够避免天然裂缝的影响,为页岩储层基质渗透率测试提供有效手段。脱气法和压汞法精度相对较低,测试结果可作为参考数据。

非稳态模拟 篇8

气动装置如气动摩擦离合器(PFC)和气动摩擦制动器(PFB)等,广泛应用于工业领域的生产实际当中,采用压缩空气作为动力源驱动零部件的运动。气动装置在工作时需经常进行进排气操作,特别是排气动作时,气缸内的气体在短时间内通过气动阀排出到大气中,会产生冲击性很强的高达125d B(A)的非稳态排气噪声,远高于国家标准规定的85d B(A)及其以下的要求,对周围环境造成严重污染,同时会严重危害车间工人的身心健康[1,2,3]。

气动装置非稳态排气噪声与常见的稳定性排气噪声有所不同,具有明显的脉冲冲击特性。该排气噪声主要由非稳态质量流引起的单极子声源和湍流四极子声源组成,在频域上呈现马鞍状,即有较高的低频和高频段的噪声声压级[4,5,6]。然而,如果仅是单纯地从时域或者频域的角度进行分析,并不能完全表征该噪声信号的特征,比如无法反映出各频率信号随时间变化的规律等[7]。因此有必要通过时频分析手段,即采用时间和频率的联合函数来表征该气动装置排气所辐射出的噪声信号,对该噪声的产生机理和辐射特性进行分析与研究,从而为工程实际中对该噪声实施有效控制奠定基础。

1气动装置排气系统结构

1.1气动排气系统结构

典型的气动系统包括气源、气缸、气动阀、不同口径的管道以及排气消声器装置等,如图1所示。气源(如空气压缩机)提供的高压气体,通过气动阀及管道流入气缸内,控制气动装置动作;随后气动阀换向或打开排气阀,使得气缸内的高压气体经由阀和管道排出。本文研究的是气动系统的排气过程,因此主要对排气通路进行研究。

1.2气动排气系统的工作原理

图1所示气缸为一个单作用气缸,这种单作用气缸被广泛应用于气动执行装置中,如机械压力机的PFC和PFB气缸等,通过高压气体的通入与排放实现摩擦盘的接合与脱离。当高压气体通入时,气缸内的活塞被高压气体推动,摩擦盘接合;而当排气时,气缸内压力降低,活塞通过复位弹簧复位,摩擦盘实现脱离动作。气动阀与排气过程的空气动力学特性密切相关,通常可以采用二位三通换向阀来实现气缸的进气和排气的切换,有时也会设置单独的排气阀或者快速排气阀来完成气缸的排气动作。图1中的气动阀为一个电气比例换向阀(如日本SMC公司生产的VY1700型气动阀)。不同于普通的二位三通换向阀,电气比例换向阀除了通过换向操作实现气缸进排气的切换外,还可以设定调节压力来控制气缸端的压力值。

气动阀VY1700通过调节先导气腔的压力值来调节进气端(P口到A口)和排气端(A口到R口)的提升阀芯的开启与闭合。当气缸端的压力值与给定先导气腔压力相等时,调节活塞保持在平衡位置,进气端和排气端的阀芯均处于关闭不导通的状态;当气缸端的压力值低于给定的先导气腔压力时,排气端阀芯闭合,调节活塞下压带动进气端的阀芯开启,导通P口和A口,使得气源提供的高压气体通入气缸;当气缸端的压力值高于给定的先导气腔压力时,进气端阀芯闭合,调节活塞上升带动排气端的阀芯开启,将A口和R口导通,使得气缸内的气体经过气动阀排出。

在排气口处安装消声器可有效降低排气噪声,本文为了分析气动装置非稳态排气所辐射噪声的特性,主要研究未安装消声器装置时的直接排气过程所辐射的噪声。气动装置排气噪声的试验测试平台基于JH23-63型机械压力机上的PFC排气系统,如图2所示。采用丹麦B&K公司的Type4189型传声器和Type2270型声级计对排气过程辐射的瞬时噪声进行测量,传声器和声级计设置在距离气动阀排气口1m,与排气口轴线夹角45°的位置。噪声信号通过研华的PCI-1712型高速多功能数据采集卡记录到计算 机中进行 处理和分 析 , 采样频率 为200k Hz。

2排气噪声辐射特性分析

2.1短时傅里叶变换

短时傅里叶变换作为最早提出的一种时频分析方法,概念清晰、物理意义明确,成为研究非平稳信号的主要方法,广泛应用于很多领域。为了达到时域的局部化,短时傅里叶变换采用加窗的方式,随时间窗的移动对信号进行傅里叶变换,得到一组局部频谱,通过局部频谱的差异来考察信号的时变特性。给定一个时间宽度很短的窗函数w(t),随着时间窗的滑动,信号x(t)的短时傅里叶变换(STFT)定义为[8]:

在实际应用中,由于采集的信号往往是离散的,因此对STFT(t,f)进行离散化处理,在等间隔的时间和频率网格点(mΔt,nΔf)处采样,以Δt和Δf分别作为时间变量和频率变量的采样间隔,则短时傅里叶变换的离散形式为:

2.2 气动装置排气噪声的时频分析

2.2.1排气噪声的时频域声压级结果

当PFC/B-63型机械压力机上的PFC气缸初始气压0.48MPa时,对通过VY1700型电气比例换向阀直接排气所辐射的噪声信号进行短时傅里叶变换,采用海明窗作为窗函数,时间间隔1.25ms,频谱分辨率100Hz,并计算时频域上的声压级,如图3所示,其中排气过程从0.045s处开始。

从这一时频域上的声压级谱图可以看出,该噪声在各个频率上均表现出随时间变化的脉冲冲击特性;在有排气流量喷出的一段时间间隔内,均表现出喷注噪声的宽频特性。

2.2.2排气噪声不同时刻的功率谱密度

为了反映出气动装置排气噪声的频谱随时间的变化规律,图4给出了在不同时间段的噪声信号功率谱图。选取的几个时间段中心点为0.04s、0.048s、0.06s、0.075s、0.1s和0.18s,分别代表了排气开始前、排气开始后、排气初期、排气流量最大、排气后期以及排气结束后的典型时刻。

从不同排气阶段的频谱声压级可以看出:在排气开始前(图4a),各频带的声压级均低于背景噪声60d B(包含数据采集中的干扰噪声),信号主要以背景的低频段噪声为主;当气动阀打开后(图4b),气流从排气口开始排出并产生排气噪声,此时主要是低频段噪声增大,高频段的增加量较小,显示在排气刚刚开始时的辐射噪声主要是以排气流量突然变化产生的单极子声源辐射为主;随着排气的进行(图4c),高频段噪声的声压级明显增大,说明在排气口处的喷注射流开始产生剧烈的湍流,并辐射出四极子的高频噪声;当排气流量达到最大值时(图4d),注意到和排气初期相比低频段的声压级基本保持不变,而高频段噪声继续增大,使得噪声的声压幅值达到瞬态峰值;其后,随着气缸内压力降低、排气流量下降(图4e),湍流噪声开始减弱,而由于质量流变化引起的单极子声源使得低频段噪声继续维持在75d B左右,随后各频段的声压级逐渐减弱直到排气过程结束;在排气结束后(图4f),高频段噪声迅速衰减,低频段声压级也逐渐恢复到背景噪声以下。

2.2.3排气噪声各频率功率谱密度随时间的变化

图5是气动装置排气噪声信号在不同频率的功率谱密度随时间变化的结果,频率值分别取为100Hz、200Hz、500Hz、1k Hz、2k Hz、4k Hz、8k Hz和16k Hz。

从图5可以看出,各频率噪声的功率谱随时间的变化有明显的起伏,其幅值脉动基本和总的噪声时域声压幅值波动接近,但各频率的噪声也具有不同的特点。和低频段及高频段的噪声相比,500Hz和1k Hz的噪声能量较低,这也体现出了气动装置排气噪声频谱的马鞍状特点;在低频段频谱中,0.1s附近出现了比较明显的峰值,分析原因为PFC气缸内活塞复位时产生的机械冲击噪声,这从单纯的时域或频域分析是很难看出的;在8k Hz处的噪声功率谱密度很高,说明通过PFC气缸通过气动阀VY1700直接排气时产生的气动装置排气噪声的能量主要集中在这一频率段,这与频谱分析的结果也是一致的;另外可发现,在排气开始到0.05s前的这一段时间内,低频率噪声(100Hz和200Hz)的功率谱密度值达到0.02Pa2s左右,而中高频率噪声(1k Hz~16k Hz)则非常小,说明在这一时间段内的噪声主要由低频段的单极子噪声组成,而高频的湍流噪声还未完全产生,这与前述分析结果相似。

总体看来,通过对气动装置排气噪声进行时频分析,可以反映出噪声信号在时间域与频域的变化及分布情况,从而更清楚地分析出排气过程与辐射噪声间的内在联系,了解噪声的产生原因和声源构成等情况。

3结论

本文来自 360文秘网(www.360wenmi.com),转载请保留网址和出处

【非稳态模拟】相关文章:

稳态系统05-14

稳态控制06-25

稳态优化07-08

稳态电能质量06-08

稳态分析论文范文05-12

稳态分析论文提纲11-15

生物稳态与环境题07-16

高三生物稳态与环境08-12

高中生物稳态与环境08-13

《稳态与环境》教学反思04-08

上一篇:专项心理能力下一篇:财会人才教育