1、1减压阀动态特性的数值模拟及故障分析刘君，董海波，徐春光 (大连理工大学 航空航天学院 辽宁 大连 116024)摘 要：采用 ANSYS 有限元软件计算得到减压阀膜片的反力时程曲线和运动部件的动态特性，在此基础上将结构简化为两个单自由度质量弹簧阻尼动力学模型，采用 Newmark 算法计算。流体控制方程为三维非定常积分形式的 ALE(Arbitrary Lagrangian Eulerian)方程，采用基于弹簧近似的动网格有限体积格式求解，应用了新的离散几何守恒律和高精度界面算法；同时采用虚拟网格通气技术实现阀门部件运动过程所引起的拓扑变化。计算表明，在较宽的上游增压速率范围内减压阀出口压力
2、存在振荡，均值接近按照静态性能设计的理论值。通过对计算流场进行分析，确定了造成两种开启故障的主要机理，修改模型参数可以排除异常。关键词：减压阀；流固耦合；非结构动网格；压力振荡；动态特性中图分类号：O354; V435 文献标识码：A 文章编号： 引言减压阀通过其内部构件调节管路系统的流量和压力，使上游的高压气体在下游出口处维持相对稳定的低压。传统设计流程中主要考虑静态性能，目前已有比较成熟的计算方法能够预测减压阀稳定工作状态下的出口压力和流量特性 1-2。由于分析手段有限，通常情况下只能依靠样品实验方法进行动态特性研究，由于阀门内部结构复杂、体积小，很难进行数据测量工作，动态试验往往只能给出
3、总体性能是否满足设计指标的定性结论 3。近年来有许多学者采用数值仿真手段开展阀门研究，根据前期调研情况看，对液体阀门管路系统中产生的水锤、空化、气蚀等现象进行仿真研究相对较多；而气体管路系统中也会出现所谓的“气柱振荡”现象，在上世纪 70 年代以后逐渐引起学术界关注，利用这种共振产生激波进行能量传导可用于制冷 4-6。减压阀的减压功能主要依靠调整流道最小截面产生壅塞来实现，由于气体可压缩性强、声速较低，在阀门内部复杂的结构中流动时常会形成激波等超声速流动现象。文献7 、8 通过数值模拟得到的结论可以看出，由于阀门的限流作用，内部流场中形成的激波导致了气动载荷剧烈变化，这是美国 Stennes
4、航天中心(SSC) 在火箭发动机试车过程中稳压阀出现故障的主要原因。文献9、10 进行固体火箭发动机调压阀的数值模拟时也发现存在马赫数高达 4 的局部区域。由于流场内存在激波，常用的商业软件中 SIMPLE 或 PISO 等基于不可压缩流动建立的算法不再适用，需要在弱解理论指导下发展能够捕捉激波的新方法。除此之外，阀门内部复杂结构网格生成、启闭引起的空间拓扑变化、描述运动部件的网格技术、流固耦合的界面算法、计算结果的验证等问题也是数值模拟应用于阀门动态性能分析遇到的挑战。图 1 是一种膜片式减压阀的结构示意图。低压腔与出口管路相连接，其下方密封膜片盒内的主弹簧 K1 使膜片 K2 向上变形，推
5、动阀芯 M02 向上运动，在阀芯和阀口(下限位 K5)之间形成流通面，高压腔和低压腔之间的压差在该处形成壅塞实现降压功能。阀芯和阀座顶盖之间的副弹簧 K3 和上限位 K6 控制阀芯运动规律，因此流通面大小受很多因素所决定。如果出口压力超过预定值，作用在膜片的气动力破坏弹簧 K1、K2 和 K3 之间受力平衡，推动膜片向下运动，带动阀芯使缝隙变窄，减少进入低压腔的气体流量，致使压力降低；反之，如果出口压力小于预定值，气动力推动阀芯使缝隙变宽，更多气体进入低压腔导致压力升高，这是实现其稳压的原理。为了避免高压冲击载荷引起的下游波动，采用阀芯与膜片中心落座可以分离的结构设计，如果膜片向下运动位移大于
6、下限位 K5，膜片脱离阀芯，副弹簧 K3 使得阀芯和阀口完全贴合，阀门完全关闭，该流通面的流量为 0。同时为了避免膜片回弹过量，膜片盒内设置了下限位 K4。阀芯和上限位 K6 之间的空腔称为卸荷腔，阀芯和阀座之间缝隙使得高压气体进入卸荷腔进而影响阀芯受到的气动载荷，为此阀芯中间开孔与低压腔相通，即使2阀门完全关闭也存在一定流量。图 1 减压器结构示意图Fig.1 Schematic of pressure relief valve基于“虚拟现实” 的优势，本文通过数值模拟的方法预测这种减压阀开启过程的动态特性，分析设计参数的影响机理，进而为工艺流程控制和产品检验提供依据。1 结构有限元模拟尽管
7、在结构分析中常用 ANSYS、FLUENT 等商用软件，这些软件中也包含流体单元的计算模型，但这些模型都是针对不可压流动方程建立，不适用于气体减压阀。由于 ANSYS 软件和流场求解器之间存在耦合问题，本文把减压阀内部运动简化为如下两个单自由度的质量弹簧模型。如图 1 示意，在阀芯与膜片的分离面处将阀门结构分为上下两部分。膜片系统的位移变量采用表示，方向向上为正，运动方程如下：1x(1)01141(/3)()()mMxKLFx式中，膜片质量 ，主弹簧刚度 、质量 、0K1M预压缩量 。膜片运动中，调整块与下限位 碰1 4L撞过程采用“接触压缩回弹脱离”非线性弹簧模型，弹簧刚度 表示为：4(2)
8、843.0(/),NmxKL在阀芯与上限位 K6 接触的全开状态作为初始位置，运动到 位置时下限位 K4 发挥10x1.8x作用。 是阀芯和膜片相互间作用力； 是作用F1F在膜片上的气动力，可以根据低压腔内流场压力积分得到，在与上部接触时去除阀芯对应的面积；是膜片变形产生的反力，由于膜片厚度很薄，1()mx当位移相对其厚度较大时，需要考虑位移大变形效应，膜片内部应力很大，可能超出其屈服限，因此，膜片变形提供的力和位移之间的关系是非线性的。对于图 1 中刚度 进行预处理建模，在流固耦合2K模拟前先用 ANSYS 软件对主阀膜片进行有限元分析，得到位移变形与反力的关系 供本文计算1()mFx时使用
9、。阀芯系统的位移变量采用 表示，运动方程2如下：(3)0232562(/)()MxCKLKxfsignF式中，阀芯质量 ，副弹簧刚度 、质量 、0133M预压缩量 。阀芯碰撞过程用非线性弹簧 和3 5近似，下限位行程为 ，初始位置6 5.8741Lm是全开状态，上限位行程为 ； 是阀20x30f芯和阀座之间涨圈的运动摩擦力。如果阀芯脱离膜片中心落座块的情况下，相互间作用力 ；如F果阀芯接触膜片一起运动，以上 2 个动力学模型合并成 1 个方程，这时 为内力，相互抵消； 是F2作用在阀芯上的气动力，受低压腔、高压腔和卸荷腔内压力影响。由于模型中受力均是非保守力，构成非保守系统，采用 Newmar
10、k 算法计算。在本文数值模拟过程中，由于三维流动的复杂性，每个迭代步内流场计算耗时远大于固体位移场计算耗时。2 流体数值方法气体的粘性比液体小，在气体阀门中，局部空间会形成超声速流动，由激波和膨胀波引起的载荷变化远远大于粘性的贡献，表现为 NS 方程中时间导数项比粘性项的影响要大。为了提高计算效率，从三维 Euler 方程出发进行数值模拟，本文的算例也验证了这种简化方式是可行的。为了描述网格变形，必须把流体控制方程变换为 ALE 形式：(4)0QdVFndSt式中， 是控制体， 是控制体边界； 为控制体边界外法向单位向量 ；守恒变量(,)xyz和对流项 为：()， uQvwe0xyztnuFn
11、Upvwae3其中 为流体相对于网格的速度矢量，UVx为流体的速度矢量， 为网(,)uvw (,)ttxyz格运动速度矢量， 为网格运动txtyan的法向速度。采用守恒变量存贮在单元中心的格心格式进行空间离散，内部复杂的空间结构均采用四面体非结构网格，在进出口和低压腔布置测点如图 2 所示。图 2 模型网格Fig.2 The grid of model假设物理量在单元内线性分布，可以构造出空间二阶精度的格式。在跨声速和超声速绕流计算中为了捕捉激波，依据 Riemann 问题提出了多种守恒型的通量计算格式。本文采用 Van Leer 矢通量分裂格式计算通量，时间离散采用二阶时间精度的四步 Run
12、ge-Kutta 方法，采用弹簧模型处理网格变形，在流场边界运动以后带动内部网格节点。采用这种动网格技术除了需要发展可靠性高的网格变形算法，还要注意采用离散几何守恒律克服网格变形给流场计算引入的误差。关于时空二阶的有限体积法推导过程、抑制非物理振荡的限制器、弹簧模型网格变形方法、虚拟挡板通气技术和新的离散几何守恒律的具体细节参考文献11。虚拟挡板通气技术的基本原理 11-12如图 3 所示，生成初始网格时在接触/分离位置或小尺度缝隙处设计一层网格，一般情况下这层网格是全场最小尺度的网格，根据运动物体在网格内部设置虚拟档板，即图 3 中阴影部分作壁面处理，其他部分作通气处理。当运动物体的位移等于
13、虚拟档板的初始网格以后，采用正常的变形网格拉伸来描述宽度增加，当拉伸距离较大导致网格质量变差，采用局部重构进行网格更新。阀门尺度为 100mm 量级，最小网格尺度大约在 0.5mm，运动行程有限，一般进行 2到 4 次重构可以满足要求。这种新技术可以计算包括接近零的无限小距离，模拟从接触到逐步远离的真实物理过程，同时也不会因网格尺度太小而影响计算时间推进步长。本文在 3 个地方设置了虚拟挡板：(1)阀芯和下限位 K5 之间缝隙设置模拟开度小于0.4041mm，直到完全关闭状态的开度变化；(2)阀芯和阀座的涨圈结构存在缝隙使得高压气体进入阀芯上方的卸荷腔，影响阀芯受到的气动载荷，卸荷腔与低压腔相
14、通，也影响出口压力，因此，在缝隙进口位置设置网格宽度为 0.45mm 的网格来模拟缝隙宽度为 0.009mm 的微小流量；(3)高压气流在进入高压腔之前经过了一个多孔网状结构的过滤器，根据前期实验数据有明显压降，也需要设置虚拟挡板。图 3 虚拟档板示意图Fig.3 Sketch of virtual baffle faces3 流固耦合界面算法和程序验证流固耦合系统的求解通常有两种方法：全耦合方法(monolithic method)和分解方法(partitioned method)。全耦合算法是指流体控制方程和结构方程作为一个系统方程同时求解，公式推导复杂，建立高精度算法困难，求解过程涉及到
15、矩阵运算，难以应用于实际工程中的多自由度结构系统，目前使用很少。分解算法是把耦合问题分解为流场和结构两部分，分别进行求解，流场和结构之间的耦合通过流固耦合界面上的信息交换来实现，又称为交错迭代算法(the conventional serial staggered CSS)。相对于全耦合算法，分解算法具有以下优点：(1)可以充分利用已有的 CFD 和 CSD 计算方法和程序；(2)可以减少程序开发的难度；(3)可以保持程序的模块化。所以，自从这一算法提出以来得到就受到4重视。为提高精度和计算效率，在流体和结构方程交错求解中又引入内迭代，内迭代达到收敛后全系统统一时间。这样一来，根据是否存在内迭
16、代，进一步把 CSS 分为紧耦合算法(tightly-coupled method)和松耦合算法(loosely-coupled method)。实际应用表明，紧耦合算法在有些情况下达到收敛标准的内迭代计算量很大，计算效率比时间小步长推进的松耦合算法还低。在文献13中结合新的几何守恒律算法，构建了高精度的松耦合算法。如下图 4 示意，高度 、长度为 的1H10L二维空间均匀分布无量纲密度 、压力 、0p声速 的静止气体，右侧连接质量 的活01am塞，活塞和固定点之间有刚度系数 的弹簧。k根据非定常一维均熵流理论可以得到活塞运动速度与压力的关系式：(5)2100()rs smukpua图 4 验
17、证模型及其结果比较Fig.4 Comparison of displacement for piston problem采用 Runge-Kutta 方法求解以上常微分方程可得到活塞运动特性，用来验证前面介绍的流固耦合计算方法的有效性。采用 2844 个网格点、5260 个单元的三角形网格离散初始空间，动网格模拟流场，在 时两种方法计算得到的活塞速度和位移比3t较如图 4，计算结果非常一致，表明本文建立的流固耦合计算程序具有良好的计算精度，为下面的应用研究打下很好的基础。4 数值计算与结果分析减压阀上游通过管路连接高压气源，下游管路有孔板限流器。工作介质为常温条件下的空气。为了避免减压阀内部压
18、力波对进出口边界影响，在三维计算域的入口和出口处连接足够长的一维等截面管路模型，一维模型与三维模型之间通过延拓两排虚拟网格点来传递流场物理量，一维模型的自由端分别提无反射入口和出口边界条件。入口长度为10m，出口长度为 2.5m。按照静态性能设计模型，主弹簧的预压力为1050N，出口压力 1.7MPa 为理论值。标准工况的上游开启总压为 21MPa，初始流场为标准大气参数 0.1MPa。在一维模型的入口处总压在给定时间内从 0.1MPa 线性增加到 21MPa，开启时间大约在10ms 到 100ms 之间。计算得到的压力随着时间的变化规律见图 5。尽管低压腔测点位置不同，曲线差异很小，表明在低
19、压腔内主要进行亚声速流动。算例基本围绕比理论值略高的均值波动，考虑到粘性影响，计算值高于理论值也是合理的。图 5 开启压力为 21MPa 时不同监测点压力变化时程Fig.5 Pressure curves of observation points at the case of 21MPa upstream pressure计算过程中在过滤器虚拟挡板两侧形成压力梯度，涨圈结构也有明显的压力差存在，说明本文采用虚拟挡板成功模拟滤网结构的阻挡作用和高压气体进入卸荷腔的流动过程。如果主弹簧预压力678N，出口压力调整到 1.2MPa 为理论值，数值模拟的低压腔测点压力规律与图 5 类似，也是围绕理论
20、值波动。以上这些论据验证本文计算是合理有效的。减压器在正常工作过程中，由于流固耦合相互干扰，膜片脱离阀芯后向下运动，与限位器碰撞后反向运动，逐渐在限位器和原点之间发生高频振荡；阀芯振幅均值较小，除了接触膜片一起向上运动，在脱离膜片后也有高频小幅振动。同时，在数值模拟过程中出现了减压器失效的故障状态，对多种工况计算结果进行总结分析，主要的故障状态可以分为两种，第一种如图 6(a)所示，5在工作一段时间以后，膜片的运动位移逐渐增大直至发散，这样可能导致膜片损坏；第二种如图 6(b)所示，减压阀可以稳定工作，但是这种工况下的出口压力低于理论值。 (a) 开启压力为 15MPa(a) The set
21、pressure P=15MPa(b) 开启压力为 2MPa(b) The set pressure P=2MPa图 6 故障状态下阀芯与膜片位移历程Fig.6 Displacement curves of valve core and diaphragm 图 6(a)为第一种故障状态的计算结果，当开启压力为 15MPa 时，减压阀阀芯和膜片以小幅振动的形式工作了一段时间后便振动发散，低压腔压力均值从 1.7MPa 变化到 7.5MPa，如图 7 所示。由这个故障工况的结果推测，采用原来设计的低压腔外形，当近距离布置出口孔板时，低压腔连通孔板前管路的容积相对较小，而高压腔压力很高，低压腔压力对
22、阀芯开度过于敏感，每一次波动都会导致低压腔压力出现一次跃升，如果压力的振动与结构的运动产生共振就会导致颤振，最终必然导致结构损坏。图 7 压力变化时程曲线Fig.7 Pressure curves of low-pressure cavity分析图 6(b)的第二种故障原因，由于阀芯有下限位 K5，因此在图中坐标原点处以下是膜片脱离阀芯后的单独运动状态，原点以上可以是阀芯单独运动也可以包含膜片的重合运动。从图中可以看出，膜片与限位器 K5 碰撞后继续向上运动，并与阀芯顶杆碰撞，碰撞时两者的动量方向相反但大小接近，因此碰撞后两者并保持接触，直到关闭状态并保持速度不变。低压腔压力最终稳定在一定值，
23、但不是设计值。由于出口孔板一直在向外界排气，低压腔压力之所以保持稳定，是因为涨圈结构的漏气进入卸荷腔再经由阀芯中心孔进入低压腔与出口排气达到了平衡。此时弹性元件的受力也达到平衡，图 8为稳定工作时的马赫数云图。这种非正常工作状态是由于建模时设定的涨圈漏气量过大造成的，说明涨圈漏气量过大会导致减压阀失效。图 8 开启压力为 2MPa 时阀门结构的马赫数云图Fig.8 Mach contours in failure state at the case of 2MPa upstream pressure上文给出的模拟结果都是近距离布置出口孔板的情况，有些工况出现故障状态。增大低压腔容积虽有助于改善
24、弹性元件的工作状态，但效果并不明显。对原始模型进行改进，模拟远距离布置出口孔板的工况，设定的孔板放置位置距离低压腔 2 米，但对于 2 米长的出口管路真实建模将会导致计算域不同方向的尺度相差太大，因此本文将 2 米长的出口管路等效成一个体积相同的容腔，置于减压阀出口附近，其模型及网格如图 9 所示.6图 9 新模型网格Fig.9 The new grid of model采用这种新模型重新进行计算，除了计算外形不同外，计算模型中还重新设定了膜片硬芯与下限位的碰撞能量损失，由原来的无损失改为损失 20%，另外卸荷腔虚拟挡板的通气量也根据设计值重新进行了设定。图 10 给出了开启压力为 15MPa
25、 时采用新模型进行计算得到的稳定工作状态下的马赫数云图。图 10 开启压力为 15MPa 新模型稳定工作时的马赫数云图Fig.10 Mach contours within the new model at the case of 2MPa upstream pressure分别针对两种故障工况，采用新模型重新进行模拟，计算结果得到明显改善，说明修改模型以后，可以抑制导致故障状态的振动发生。图 11 压力变化时程曲线Fig.11 Pressure curves of low-pressure cavity从图 10 可以看出，更改卸荷腔虚拟挡板通气量后，顶杆中心孔内的流速非常小，高压气体几乎不
26、通过涨圈结构进入卸荷腔；从图 11 可以看出，开启压力为 15MPa 和 8MPa 的工况在稳定工作时低压腔压力没有振动而是接近理论值 1.7MPa 左右，开启压力为 21MPa 的工况低压腔压力呈现出有限幅值的波动。低压腔压力趋于稳定的工作状态，将不会再出现颤振现象，有效的防止了减压器的结构出现损坏。5 结 论以上研究表明：(1)在两种主弹簧预压力的大部分算例中，计算得到的出口压力的均值接近按照静态性能设计出的理论值，基于非结构动网格的流固耦合计算方法有效开展减压阀的动态特性研究，采用虚拟网格技术实现了减压阀从完全关闭到开启的真实工作过程。(2)在进口气体的冲击载荷作用下，这款减压阀正常工作
27、的出口压力表现为周期振荡，而且对于进口气体的不同增压速率有很好的适应性。(3)减压阀会出现气动力与动力学系统相互干扰形成的颤振现象，定性分析这种流固耦合现象是导致颤振的主要机理，可采取相应措施可以消除颤振。(4)由于涨圈结构的设计问题，导致涨圈漏气量过大造成减压阀失效，修改模型参数可以消除这种故障，在工艺控制时需要关注缝隙宽度的检测。参考文献：1 谭建国,江燕平,王振国.减压器增压过程中的振荡 抑制措施J.推进技术,2013,34(1):145-149.TAN Jian-guo, JIANG Yan-ping, WANG Zhen-guo. Suppression approach on os
28、cillation of pressure relief valve during pressurization processJ. Journal of Propulsion Technology, 2013,34(1):145-149.2 叶奇昉,严诗杰,陈江平,等.气动先导式电磁阀的自激振动J.机械工程学报,2010,46(1):115-121.YE Qi-fang, YAN Shi-jie, CHEN Jiang-ping, et al. Self-excited vibration in a pneumatic pilot-operated solenoid valveJ. Jour
29、nal of Mechanical Engineering, 2010,46(1):115-121.3 戴梧叶,马彬,张国舟,等.减压器特性实验研究J.北京航空航天大学学报,1999,25(6):711-713.DAI Wu-ye, MA Bin, ZHANG Guo-zhou, et al. Experimental studies on pressure regulator characteristicJ. Journal of Beijing University of Aeronautics and Astronautics, 1999,25(6):711-713.4 孙建红,纪益明,
30、明晓.开口圆管气体的非线性振荡研究及雾化应用J. 空气动力学学报,1997,15(2):177-184.SUN Jian-hong, JI Yi-ming, MING Xiao. Investigation of 7nonlinear gas oscillations in open pipes and application in atomizationJ. Acta Aerodynamica Sinica, 1997,15(2):177-184.5 李学来,郭荣伟.振荡管内气柱谐振的研究J.航空学报, 1999,20(2):97-99.LI Xue-lai, GUO Rong-wei. O
31、n resonance of gas column in an oscillatory tubeJ. Acta Aeronautica et Astronautica Sinica, 1999,20(2):97-99.6 刘君,郭正,瞿章华.闭口管内气体非线性振荡研究J. 推进技术,2000,21(4):32-35.LIU Jun, GUO Zheng, ZHAI Zhang-hua. Numerical simulation of nonlinear gas oscillations in closed pipeJ. Journal of Propulsion Technology, 200
32、0,21(4):32-35.7 Vineet A, Ashvin H, Peter C, et al. Simulations of Instabilities in Complex Valve and Feed SystemsR. AIAA-2006-4758.8 Jeremy S, Ashvin H, Vineet A. Unsteady Analyses of Valve Systems in Rocket Engine Testing EnvironmentsR. AIAA 2004-3663.9 李哲,魏志军,张平.调压阀内流场数值模拟及动态特性分析J北京理工大学学报,2007,27
33、(5):390-394.LI Zhe, WEI Zhi-jun, ZHANG Ping. Internal flow field numerical calculation and dynamic characteristic study of pressure-regulating valveJ. Transactions of Beijing Institute of Technology, 2007,27(5):390-394.10 李哲,魏志军,张平.燃气调压阀内流场三维数值模拟J.固体火箭技术, 2007,30(3):210-213.LI Zhe, WEI Zhi-jun, ZHAN
34、G Ping. 3D numerical simulation on internal flow field for pressure regulating valve of gas generatorJ. Journal of Solid Rocket Technology, 2007,30(3):210-213.11 刘君,白晓征,郭正.非结构动网格计算方法及其在包含运动界面的流场模拟中的应用M.长沙:国防科技大学出版社,2009.12 郭正,刘君,陈广南.安溢活门动态特性的流固耦合数值分析J.机械工程学报,2008,44(8):61-66.GUO Zheng, LIU Jun, CHEN
35、 Guang-nan. Numerical analysis of dynamic properties of safety valve by using method of fluid and structure interactionJ. Journal of Mechanical Engineering,2008,44(8):61-66.13 郭攀,刘君,武文华.爆炸冲击载荷作用下流固耦合数值模拟J.力学学报,2013,45(2):283-287.GUO Pan, LIU Jun, WU Wen-hua. Numerical modeling for fluid-structure in
36、teraction under blast and impact loading response J. Chinese Journal of Theoretical and Applied Mechanics, 2013,45(2):283-287.Numerical simulation of dynamic properties of pressure relief valve and fault diagnosisLIU Jun, DONG Hai-bo, XU Chun-guang(School of Aeronautics and Astronautics, Dalian Univ
37、ersity of Technology, Dalian, 116024)Abstract：Using ANSYS to calculate the force response curves with time history of pressure relief valve and dynamic characteristics; Solid structure was simplified to two single degree of freedom mass spring damper dynamic models which could be calculated by Newma
38、rk algorithm. The three-dimensional unsteady integral forms of ALE equation were solved by finite volume scheme based on spring analogy method and dynamic grid. A new discrete geometric conservation law and high-precision fluid and solid coupling algorithm were used to calculate the governing equati
39、ons; the virtual mesh ventilation method was adopted to resolve the problem when multi bodies moved from contact to separation and the mesh topology changed. The results show that the outlet pressure of valve oscillated within wide range upstream pressurization rate, the average pressure is close to
40、 the theoretical value of static performance design. Through the flow field analysis, it can determine the main mechanism of two kinds of open fault and modify the model parameters to exclude anomalies.Key words：pressure relief valve; fluid and solid coupling; dynamic unstructured grid; pressure oscillation; dynamic properties作者简介：刘君(1965)，男，教授。电话：(0411)84707176；E-mail：
