数值仿真模型（精选9篇）
数值仿真模型 篇1
海洋工程领域中的浮式结构物在海域展开应用时, 会不可避免地遭遇各种恶劣海况。浮式结构物受到海浪产生的绕流力、惯性力、冲击力、浮托力等作用, 这些作用直接关系到结构物的生存安全。因此在浮式结构物设计时, 必须要充分考虑海浪的影响, 精确计算海浪载荷, 校验和优化结构, 研究浮式结构物的耐波性和稳性, 保证其满足设计要求的稳定性和安全性。
目前评估海浪载荷对浮式结构物的影响主要采用物理模型试验的方法, 然而这种方法费时费力, 而数值模拟具有参数设置灵活、计算结果精确等优点, 正逐步成为设计领域研究海浪载荷的重要手段。由于受到计算机硬件发展水平和波浪理论不成熟等的制约, 早期的海浪数值仿真主要以二维为主, 但是对于需要研究海浪和浮式结构物相互作用过程中产生的波浪折绕射、漩涡等现象时, 显然不能满足工程的需要。因此研究三维海浪数值模型, 实现对海浪现象更加真实准确的描述, 是海上复式结构物设计领域中研究海浪和结构物相互作用的必然发展方向。
本文利用谱分析的方法, 在MATLAB环境下对三维随机海浪模型进行了数值模拟仿真, 并给出了三维随机海浪波面图, 为浮式结构物设计中计算海浪载荷提供了参考。
1 Longuet-Higgins长峰波海浪模型
海浪现象十分复杂, 具有随机性和非线性的特征, 因此要建立精确的海浪模型是十分困难的。根据海浪谱和随机海浪理论, 可将实际海浪看成为不同频率、不同传播方向、不同波高和不同初始相位的正弦波叠加的结果[1]。Longuet-Higgins模型就是把长峰波海浪的波动用无数个随机余弦波的叠加来描述, 其波幅的表达式为:
其中ξai、ωi和εi分别表示第i个余弦波的波幅、频率和相位。
在实际应用中, Longuet-Higgins模型常用海浪的频谱来表达, 即:
ζ (t) =i=∑∞12SζΔωcos (ωi+εi)
式中Sζ (ωi) 为海浪的频谱。目前应用较多的海洋频谱有Pierson-Moscowitz谱 (PM谱) 、Neumann谱 (N谱) 、ITTC双参数波谱等, 其中PM谱应用最为广泛, 能很好地代表实际充分发育的随机海浪, 故本文采用PM谱。PM谱的表达式为:
其中v为离海面19.5m处的平均风速。
平均风速为10m/s、12m/s、15m/s的PM谱的仿真曲线如图1所示。从图1中可以看出, 海浪的PM谱为窄带谱, 能量集中于某一频段内, 随着风速的增加, 即高海情情况, PM谱峰值增大, 能量更为集中, 如当风速达到15m/s时 (六级海情) , 谱峰值为3.6m2∙s, 能量集中于0.3-1.7Hz频带内。因此对海浪进行数值仿真时, 可以选取能量集中的频段内的有限个谐波成分进行叠加, 可减少计算量, 从而加快仿真速度。
2 三维不规则短峰波随机海浪仿真
2.1 基于谱分析的三维不规则短峰波随机海浪模型
Longuet-Higgins模型假设海浪只沿一个固定方向传播, 且假设其波峰和波谷线相互平行且垂直于前进方向。实际上, 海面上的波浪不仅波高不同, 频率不同, 还会向各个方向传播。这些谐波除产生在主风向上的主浪向以外, 在主浪向两侧的±π2角度范围内都会有谐波的扩散, 这样的海浪称为三维不规则短峰波海浪。
与Longuet-Higgins模型不同, 三维不规则短峰波海浪是由振幅不等、频率不等、初始相位不同、传播方向不同的谐波多重求和得到。双叠加模型[2]为:
其中ζaij、ωi、ki、μj、εij分别为组成谐波的波幅、角频率、波数、方向角和随机初相位, εij为[0, π2]内的随机数, (ξ, η) 为波面上某点的坐标。
根据谐波的波幅和频谱的关系, 三维不规则短峰波海浪的波面方程可进一步写为:
其中Sζ (ωi, μj) 为短峰波海浪的方向波谱。
通常认为波能的方向和频率是无关的, 则短峰波海浪方向波谱可用两个独立函数的乘积来表示, 即:
其中Sζ (ω) 为海浪的PM谱, μ为海浪的方向角, ϕ (μ) 为波能的扩散函数。ITTC提出的建议性形式为:
波数ki表示谐波的传播方向上2π距离内波的个数, 于是有ki=2πλi, 其中λi为第i个谐波的波长。因为波长不易获得, 由线性波动理论可知, 在不考虑表面张力的情况下, 数值上的波数和角频率有如下关系, 即:
式中g为重力加速度, D为水深。
2.2 三维不规则短峰波随机海浪仿真
某浮式结构物的作业和生存海域的海洋条件为:风向30°, 风速13.8m/s, 五级浪, 为校验该结构物所能承受的波浪载荷, 需要对该海情下的海浪进行模拟。
选择100m×100m大小的海域, 利用上述理论进行三维随机波浪仿真。各参数设置为:仿真频段[0.3Hz, 1.7Hz], 海浪的方向角μ=30°, 风速v=13.8m/s, 时间t=10s。图2为该海域的模拟海浪图。
该仿真海浪中, 最大波高Hmax=3.32m, 根据标准浪级波高的参考值[3], 五级浪对应的波高范围为[) 2.5, 4.0 m, 最大波高Hmax位于允许的波高范围内, 说明利用海浪谱来模拟三维随机海浪能够得到比较精确的海浪波面图和波高值。进一步根据流体的势流理论就可以分析计算出该结构物受到的海浪载荷, 为校验结构物的结构强度提供了必要的基础。
3 结论
海上浮式结构物结构强度校验需要计算分析海浪载荷, 该文利用海浪谱分析的方法, 实现了在开阔海域主要由风力引起的海浪的模拟, 该仿真海浪的波面图和波高符合标准浪级波高的参考值。进一步利用流体的势流理论就可以分析计算出结构物受到的海浪载荷, 为进一步的结构强度校核提供基础。
参考文献
[1]俞聿修.随机波浪及其工程应用[J].大连理工大学出版社, 1992.
[2]Khatri S K.In the search of a coastal ocean wave model.OCEANS 97, MTS/IEEE Conference Proceedings, 1997, 213-218.
[3]Kinsman B.Wind waves.Prentice Hall In., 1965.
数值仿真模型 篇2
航空母舰尾流数值仿真研究
为研究舰载飞机进舰着舰过程中航空母舰尾流的影响,选取了一个尾流模型在微机的MATLAB环境中模拟其速度的空间分布.根据成因将紊流分为4个组成部分:自由紊流分量的.功率谱用成形滤波器法模拟,稳态分量通过线性插值得到,周期性分量可直接计算得到,而航母诱导的随机紊流分量则通过MATLAB的SIMULINK工具包模拟.模拟表明周期性分量占主导地位.此模型需要实验数据和船上实测数据进一步确认.
作 者：彭兢 金长江 PENG Jing JIN Chang-jiang 作者单位：北京航空航天大学,飞行器设计与应用力学系刊 名：北京航空航天大学学报 ISTIC EI PKU英文刊名：JOURNAL OF BEIJING UNIVERSITY OF AERONAUTICS AND ASTRONAUTICS年，卷(期)：26(3)分类号：V211.78关键词：航空母舰 舰载飞机 尾流 紊流 功率谱
数值仿真模型 篇3
假定竹材为各向同性材料, 利用MIDAS/Gen建立数值模型。通过调整数值模型竹材的弹性模量可使数值模型与实际模型的自振频率相等;通过调整数值模型阻尼比可使数值模型与实际模型的加速度时程数据最大值相等。由此方案建立的数值模型为进一步的结构模型动力学研究提供了参考。
1振动台试验
1.1试验参数
1.1.1试验模型
模型用0.35 mm厚本色侧压双层复压竹皮和502胶水制作, 胶水层厚度为0.1 mm, 考虑胶水层的竹材 (以下简称竹材) 容重测试结果为9.3×10-6N/mm3。竹质框架模型共三层, 底层层高240 mm, 二层和三层层高均为480 mm, 平面柱网尺寸为240mm×200 mm, 柱与梁均为箱型截面 (长×宽×厚=10 mm×10 mm×0.8 mm) , 楼面板和屋面板均为0.8 mm厚, 模型立体网线尺寸及楼层定义见图1。用热熔胶将竹质框架与330 mm×330 mm×8 mm竹质底板连接, 竹质底板与振动台之间用螺栓连接。振动台为北京波谱世纪科技发展有限公司生产的8通道WS—Z30—50小型精密振动台, 试验分别测试台面、模型三层楼面和顶面的加速度时程数据。试验装置照片见图2。
1.1.2荷载和地震作用
各楼层铁块荷载为:台面和二层楼面空载, 三层楼面4.661 kg, 顶面2.238 kg。
地震作用采用单向加载, 加载方向为X轴方向 (见图1) 。输入的地震波为汶川波, 取自2008年汶川地震中什邡八角站记录的南北方向加速度时程数据, 原始记录数据点时间间隔为1/200 s,
全部波形时长为205 s, 峰值加速度5.81 m/s2。模型试验截取原始记录中10~42 s区间内的数据, 通过调整数据点时间间隔和加速度幅值得到三种不同振动强度的汶川地震波:第一级汶川波数据点时间间隔为1/200 s, 峰值加速度为0.57 m/s2;第二级汶川波数据点时间间隔为1/250 s, 峰值加速度为0.86 m/s2;第三级汶川波数据点时间间隔为1/300 s, 峰值加速度为1.00 m/s2。其中, 调整后的第三级汶川波加速度时程曲线见图3。
1.2试验加载工况
在进行汶川波加载前, 先用白噪声激振法测试框架结构模型 (含铁块荷载, 简称实际模型) 的自振频率[2], 白噪声激振频率为0~20 Hz, 激励通道为台面通道。试验进行了3种工况的测试。
工况1:模型自重+铁块荷载+第一级汶川波;
工况2:模型自重+铁块荷载+第二级汶川波;
工况3:模型自重+铁块荷载+第三级汶川波。
1.3试验结果
振动台输入白噪声测试得到实际模型的基本频率为1.59 Hz, 第二频率为4.64 Hz (见图4) 。
为了便于与数值仿真结果的加速度时程曲线走势比较, 给出测试的各工况顶面加速度时程曲线。测试的各工况三层楼面和顶面峰值加速度值及其发生时刻见表1。测试的工况1、工况2和工况3实际模型顶面加速度时程曲线分别见图5~图7。
2数值仿真模拟
采用MIDAS/Gen软件建立数值模型, 模型几何尺寸、竹材容重和铁块荷载取值均与实际模型相同, 竹材的弹性模量初始值和所有振型的阻尼比初始值分别取10 000 N/mm2和0.012, 底端四个柱脚边界条件按固结考虑。在小型振动台试验中, 由于上部结构的作用, 台面加速度时程数据实际测试值与理论输入值不相等。为了使数值仿真模拟结果更可靠, 各工况仿真分析所取的地面加速度分别采用振动台试验测试的台面加速度时程数据。建模后, 首先作特征值分析, 然后输入地震波进行线性时程分析。
由于各工况测试的加速度均为绝对加速度, 为便于直接对比, 三层框架模型仿真分析的振动方程用绝对坐标表示为
忽略模型制作误差, 在认为数值模型与实际模型的几何尺寸和荷载相同的条件下, 根据各参数对自振频率计算结果的影响[3], 若要求各振型数值模型与实际模型的自振频率对应相等, 则仿真分析的可调参数仅限于竹材的弹性模量。另外, 为了使数值模型与实际模型对应部位的最大值加速度值相等, 则可以通过调整数值模型所有振型的阻尼比来实现。仿真分析流程见图8。图中, f0为计算的数值模型自振基本频率, f为测试的实际模型自振基本频率, amax0为计算的数值模型顶面最大绝对加速度, amax为测试的实际模型顶面最大绝对加速度。
利用MIDAS/Gen按图8操作, 竹材的弹性模量最终调整结果为8.805×103N/mm2, 此时, 计算的数值模型在X方向的自振基本频率等于测试的实际模型自振基本频率值1.59 Hz, 计算的数值模型在X方向的第二振型自振频率为4.64 Hz, 也与测试结果吻合;工况1的最终调整结果为ζ=0.027 5, 工况2的最终调整结果为ζ=0.028 8, 工况3的最终调整结果为ζ=0.040 6。计算的工况1、工况2和工况3数值模型顶面加速度时程曲线分别见图9~图11。
各工况数值模型与实际模型的三层楼面和顶面峰值加速度及其发生时刻对比见表1。由数值模型与实际模型顶面加速度时程曲线的对比和表1可知, 数值仿真模拟的精确度较高。
在数值仿真模拟结果比较精确的情况下, 各工况阻尼比并不相同, 其值随地震波振动强度的增大而增大, 文献[4]的实验数据和理论分析同样证实了这一结论, 阻尼比取值与振动强度的关系曲线见图12。
3结论
数值仿真模拟方案建立的数值模型为进一步的结构模型动力学研究提供了参考。通过三层竹质框架模型振动台试验的数值仿真模拟得出以下结论:
(1) 在弹性工作条件下, 分别通过调整数值模型竹材的弹性模量和阻尼比, 使数值模型与实际模型的自振频率和加速度时程数据最大值相等可以得到较高精度的动力学仿真效果。
(2) 数值仿真模拟工作间接确定了竹材弹性模量和结构阻尼比, 确保了弹性模量和阻尼比取值的有效性。
(3) 在数值仿真模拟结果比较精确的情况下, 证实了结构的阻尼比取值随地震波振动强度的增大而增大的结论。
摘要：为了确保进一步研究设置耗能减振装置结构仿真分析的可靠性, 首先借助于小型振动台, 对无耗能减振装置的三层竹质框架模型进行了振动台试验, 测试了实际模型的自振频率和不同测点加速度时程数据。其次利用有限元计算软件MIDAS/Gen建立竹质框架数值模型, 分别通过调整数值模型竹材的弹性模量和阻尼比, 使数值模型与实际模型的自振频率和加速度时程数据最大值相等的方式进行数值仿真研究。仿真研究得到的数值模型与实际模型的加速度时程曲线吻合良好, 证实了结构的阻尼比取值随地震波振动强度的增大而增大的结论。
关键词：竹质框架,振动台试验,加速度,仿真模拟
参考文献
[1] 蒋玉川, 傅昶彬, 阎慧群.MIDAS在结构计算中的应用.北京:化学工业出版社, 2011:177—220
[2] 庄表中.白噪声及其作为激振信号的应用.力学与实践, 1987; (5) :32—35
[3] 王社良.抗震结构设计.抗震结构设计 (第4版) .武汉:武汉理工大学出版社, 2012:37—47
数值仿真模型 篇4
多轴振动试验系统传递函数估计的数值仿真
文章应用频响矩阵建立了多轴振动系统输入输出之间的`数学模型.通过数值仿真的方法,研究了多轴振动试验中传递函数估计的方法,针对不同激励信号组合、不同估计算法和使用多次平均技术情况下进行了数值模拟,并对结果进行了分析.可以为多轴振动试验控制系统的研究提供参考.
作 者：魏军 冯咬齐 樊世超 邱汉平Wei Jun Feng Yaoqi Fan Shichao Qiu Hanping 作者单位：北京卫星环境工程研究所,北京,100094刊 名：航天器环境工程 ISTIC英文刊名：SPACECRAFT ENVIRONMENT ENGINEERING年，卷(期)：200926(1)分类号：V416.2 V416.8关键词：多轴振动试验 传递函数估计 仿真
数值仿真模型 篇5
研究人体上呼吸道气流流场特性，对探索气溶胶在人体上呼吸道内的扩散及沉积模式和分析大气环境状况等因素对人体健康的影响有非常重要的意义。近年来，国内外研究人员采用数值仿真的方法对人体上呼吸道气流流场特性进行了大量有意义的研究工作。徐新喜等[1,2,3]通过建立三维的口喉模型，数值模拟了循环呼吸条件下气流组织特性，进行了初步的实验研究;Wang等[4]采用数值方法，模拟了人体上呼吸道内的气流组织形式;Nithiarasu等[5]也对人体上呼吸道中的层流和湍流等形式进行了仿真研究;Wlofgang等[6]对考虑壁面弹性的人体前四级支气管内气流运动进行了数值模拟研究;X.G.Cui等[7]对理想的人体口喉模型内的气流组织进行了精细的大涡模拟，指出喉部射流将和二次涡流在射流边界处产生相互作用，且喉部射流尾部将发生严重变形。但是目前研究也存在不足，一是采用的模型多数为简化模型，与真实人体上呼吸道结构差异较大;二是大多数数值模拟模型采用低Reynolds数k-ω模型[8,9]和标准k-ε模型[10]来描述气体流动情况，缺乏更高精度的仿真计算。同时，在以往的研究中，通常假设呼吸道壁面是刚性的。流固耦合力学的发展和成熟能够克服这一假设存在的缺陷，使数值仿真的理论分析与呼吸道内的实际气流情况更为相符，从而开拓研究呼吸道内气流运动的思路和方法。基于此，本文运用流固耦合力学的方法，对循环呼吸模式下人体口喉模型内气流和呼吸道壁的相互耦合作用进行了数值模拟，系统地分析了流固耦合作用下口喉模型内的气流运动特点。
2 数值仿真模型和方法
2.1 仿真模型
参考ARLA(aerosol research laboratory of alberta)模型[11]和Stapleton K.W模型[12]的尺寸，文中建立和采用如图1所示的口喉仿真模型。
2.2控制方程
2.2.1流体控制方程
文中采用流体湍流模型为Menter’s SST两方程模型，即耦合了LangtryMenter的转捩模型。该模型基于间歇因子和转捩动量厚度雷诺数2个输运方程，前者与SST湍流模型耦合，后者用来捕捉非当地湍流强度的影响，其张量表达式[13,14]为:
间歇因子输运方程:
式中生成Pγ项的表达式为：
其中，S为应变率张量，Flength用以调整转捩区长度，Fonset用以启动Pγ，包含了“当地转捩”判据，当该比值大于1时，Pγ启动即γ开始增长。
转捩动量厚度雷诺数输运方程：
式中的源项定义如下：
2.2.2固体（壁面）结构控制方程
本文避免了传统的刚性壁面假设，而是将口喉模型壁面材料假设为线弹性材料，这样与实际情况更为接近。参照文献[15]，设定壁面的杨氏弹性模量为9 MPa，泊松比为0.4，为均匀各项同性材料，即服从胡克定律的薄壁弹性模型。
气管壁面方程可表达为：
式中：σs是应力张量，ρs是密度，as是加速度。
在流固耦合交界面上，应当满足以下条件：
式中，d为位移，n为边界法向，下标s、f分别表示固体和流体。
2.3 边界条件和计算方法
假设口腔入口气流流量随时间按正弦函数形式周期性地变化，且气流在入口边界是匀速的，方向为入口平面的法向方向，呼吸流量为30 L/min，呼吸频率为15循环/min。支气管出口采用压力出口边界条件，压力为0 N，其余量梯度为0。在一个循环周期内，选取呼吸吸气和呼气状态下的峰值T3和T7的仿真结果分别进行分析，其中T3=1 s,T7=3 s，如图2所示。
3 结果与讨论
3.1 吸气时刻仿真结果分析
如图3所示，循环吸气过程中呼吸流量达到最大值T=1 s时，口喉模型中不同切面内气流轴向速度分布，其中A-A’为冠状面，其余为水平面。图3(a)显示的是口喉模型的正中矢状面(将口喉模型分为对称的左、右两部分的纵切面)的气流流场分布结果。观察发现，气流在口喉模型内的咽部和喉部，形成了2个速度增长点。在喉部，气流速度达到全场最大值12.33 m/s，而在咽部气流发生分离现象，在咽部外壁形成分离区，如图3(d)所示。咽部外壁分离区的产生，将使气溶胶在此部位回旋、滞留，导致其在咽外壁的沉积几率增加。图3(e)、(f)显示，湍流喷射现象出现在了声门位置，接着气流在声门下方的气管内出现流动分离，这主要是声门的形态结构所致。图3(f)显示声门下方近气管内壁气流速度较高，而图3(g)结果表明，与声门距离增加时气管内壁气流速度仍高于外壁，但气流速度差，并呈现出逐渐减小的趋势。气流速度高的部位，其内部气溶胶颗粒的惯性碰撞作用多，因此将造成多数气溶胶颗粒在气管内壁沉积。
3.2 呼气时刻仿真结果分析
如图4所示，为循环呼气过程中口喉模型内呼吸流量达到最大值T=1 s时刻不同切面内的轴向速度分布图，其中图4(a)显示的是口喉模型的正中矢状面(将口喉模型分为对称的左、右两部分的纵切面)的气流流场分布结果。从图中可以看出，在y=0截面内，在咽部和喉部由于受到截面收缩和气流运动方向改变的影响，分别在咽部前壁、喉部后壁和声门后部形成高速区，气流速度的最大值为13 m/s，如图4(c)～(e)所示。气流速度的增加，容易造成气溶胶颗粒在咽部和喉部的沉积。图4(b)显示，与吸气阶段不同，在呼气阶段在口腔内部发生流动分离现象，形成分离区，分离区的产生使得气溶胶颗粒在口腔中的沉积几率增加。在呼气阶段，由于呼气气流与吸气阶段相比气流流向发生改变，在咽部和喉部未出现气流分离现象，而且气管内气流流速分布比较均匀，在气管不同截面内速度梯度不明显，如图4(f)、(g)所示。如图5所示，为吸气阶段和呼气阶段各个时刻流场内的最大速度曲线。从图中可以看出，呼气阶段气流最大速度略高于吸气阶段，对气道的冲击力也将更强一些。分析原因，一方面是由于气流流向不同，另一方面则与气道的膨胀与收缩有关。
4 结论
盾构的数值仿真分析 篇6
1.1 ANSYS简介
ANSYS是目前广泛应用于机械工程领域的数值分析软件,通过在计算机上模拟现场的各种条件,通过节点间的连接模拟构件的实际受力情况,通过加载受力,在计算机上提前预测各种不利效应,从而判断该结构是否满足安全、使用要求。
1.2 ANSYS软件主要功能
ANSYS软件能模拟大部分物理力学作用,可广泛的用于地下工程、航空航天、机械摩擦和制造、物理碰撞、桥梁施工、温度效应、大型水利工程、等诸多专业学科。能对结构进行自适应单元划分,从而进行非线性分析,通过对参数的调整,确定各影响因素的敏感性,对各种设计方案进行优化,从而在实际应用中,通过各种物理化为作用,对敏感度高的影响因素重点改进,对耗费经济少,功效较高的影响因素重点优化,确定最优设计方案。
在盾构施工中,由于土层的受力比较复杂,根据经验难以断定地层的损失量,采用过于保守的设计会加大工程造价,通过ANSYS模拟盾构隧道施工全过程,可以有效的控制工程造价,保证地下工程施工安全。
2 盾构隧道分析
采用四结点单元划分网格如图2所示:
首先我们研究盾构开挖时所产生的沉降变形,由于衬砌单元具有很强的刚度,开挖过程中并没有衬砌单元的存在,于是先采取单元生死技能杀死衬砌单元便可以得到,要激活“单元死”的效果,ANSYS程序并不是将“杀死”的单元从模型中删除,而是将其刚度矩阵乘以一个很小的因子[ESTIF]。从而不对载荷向量生效。
从云图4可以看出,地表沉降在隧道盾构正上方处沉降量最大,离正上方水平距离越远,沉降越来越小,成正态分布,下面采用peck理论公式计算进行对比。
Peck公式为:Sx=Smaxexp(-x2/(2i2))
式中:Sx——x处的沉降量。
x——离中线的水平距离;
Smax——中心处的最大沉降值;
i——宽度系数。
最大沉降量采用下式估算:
Vs——盾构引起的地层损失;
i——沉降槽的宽度系数。
计算结果显示,沉降呈正态分布,根据数据作图如5所示。
由云图6可以看出,在盾构周围有较大的应力,为了保证盾构隧道的强度,及时给盾构后的四周支护衬砌。激活衬砌单元,使衬砌的刚度发生作用,目前的地铁施工中,采用拼装管片的方法使圆形管片形成封闭的衬砌环,可以得到后续第二步变形如图7所示。
从云图7可以看出,衬砌的激活对沉降并没有什么影响,衬砌的目的是用了加强隧道壁的刚度的。
按一层楼1.8KN/m2的恒重加载,分别在隧道上方加5层楼和3层楼的荷载,每层楼按100m2折算。
由图9可以看书,加载后,地面沉降有所加强,地铁通常在城市施工,建筑物较为密集,周边建筑物产生的荷载要大的多,为了减小盾构隧道对地面沉降的影响,减小地面建(构)筑物及地下管线受到的破坏,我们对土体进行加固处理,使土层的承载力增大,粘聚力变大,抗剪切强度得到提高,摩擦角增大,实际施工中有采用高压旋喷桩、水泥搅拌桩、高压注浆等方法。
根据图10、11对比未加固土层的云图可以看到,土层经过加固后,变形有了改善,只是在隧道正上方可以看到比较明显的变形,符合peck公式的理论分析结果。
摘要：随着城市经济的发展,地铁等地下工程的施工越来越多,盾构施工使周边建筑物及管线的基础应力得到释放,或由于附加应力的影响引起地基下沉或隆起,这种变形对地面建筑物影响程度与有效间距、施工范围与建筑物间的土质情况等有关。
关键词：地下工程,地基沉降,盾构施工
参考文献
[1]朱伟(译).隧道标准规范(盾构篇)及解说.北京:中国建筑工业出版社2004.196～172.
[2]覃仁辉,隧道工程.重庆:重次大学出版社.2001,1～6.
[3]沈培良,张海波,殷宗泽.上海地区地铁隧道盾构施工地面沉降分析[J].河海大学学报(自然科学版),2003,31.556～559.
[4]张海波,殷宗泽,朱俊高.地铁隧道盾构法施工过程中地层变位的三维有限元分析.岩石力学与工程学报,2005,24期755～760.
火炮射击稳定性数值仿真 篇7
随着近年来轮式火炮、车载火炮和超轻型火炮的迅猛发展, 火炮轻量化趋势越来越明显。轻量化火炮的核心特性就是射击稳定性, 这通常需要设计高效反后坐装置和炮口制退器来实现[1,2,3]。火炮射击过程中翻转方向要发生变化, 因而稳定条件也要发生变化, 另外翻转的驱动力 (后坐阻力和复进合力) 变化很复杂, 这导致射击稳定性变化也十分复杂。本文针对某型地面火炮车轮着地射击的工况 (稳定性问题最突出) , 基于后坐部分运动微分方程、射击过程稳定条件和稳定参量变化公式基础上的稳定性分析数学模型, 对射击过程稳定性时变特性进行了数值仿真, 重点考虑了后坐行程和火炮射角对稳定性的影响。
1 稳定性分析数学模型
该数学模型涉及内容较多, 对火炮设计理论教材中已有内容只作简要介绍[4,5,6], 重点阐述教材中讨论较少或者没有涉及的问题, 如后坐复进驱动力一体化描述、稳定参量变化公式等。
1.1 后坐部分运动微分方程
其中, V为后坐速度, R为后坐阻力。
其中, U为复进速度, r为复进合力。
在求解运动微分方程时, 可将后坐和复进运动微分方程合并为
其中, V为后坐部分的速度, X为后坐部分行程。通过求解式 (6) , 可以得到后坐部分在后坐和复进过程中的V和X, 继而求得相关各力 (如R、r等) 的大小。这些力决定了火炮射击过程中的稳定性。
1.2 射击过程稳定条件
后坐时, 全炮受力分析如图1所示, 其中Qz为火炮战斗全重, NA和NB分别为车轮和驻锄处对火炮的垂直支撑力, TA和TB分别为车轮和驻锄处产生的水平阻力, 为后坐部分惯性力。
将Ppt向惯性力作用线简化后, 全炮受力情况如图2所示。后坐时火炮有绕驻锄后翻的趋势, 其稳定条件为
其中, D为驻锄支点B到全炮重心的水平距离, h为R对驻锄支点B的力臂, e为后坐部分质心到炮膛轴线的距离, 规定质心在炮膛轴线上方e为负, 图中刚好是e为负的状态。
复进分为加速复进和减速复进两种情况, 加速复进时期, 全炮受力如图3所示。该受力情况与图2基本一致, 只是Ppt·e=0, R成为了r, 其稳定条件仍然可以采取式 (8) 的形式。
复进减速时期, 全炮受力如图4所示。此时火炮有绕A点前翻的趋势, 故稳定条件为Qz (L-D) ≥rhf (9) 其中, L为车轮与地面接触点A到驻锄支点B的水平距离, hf为力r到A点的力臂。
总之, 火炮射击过程中, 稳定条件有两个, 式 (8) 和式 (9) 分别对应火炮后翻和前翻两种情况, 后翻时驱动力方向向后, 前翻时驱动力方向向前, 可统一用驱动力Qr来表征, 后坐时驱动力为R, 复进时驱动力为r。
1.3 稳定参量变化公式
射角变化和后坐行程变化时, 稳定条件中各参数也要发生变化, 这是研究稳定时变特性必须考虑的。
1.3.1 全炮重心随射角的变化
如图5所示, O为摇架耳轴中心, Lr为耳轴到驻锄支点的水平距离, C为火炮平射时火炮起落部分质心, C′为射角为φ时火炮起落部分质心, OC=OC′=l1;平射时, OC连线与水平线夹角为δ, 射击前全炮重心到驻锄支点水平距离为D00;射角为φ时, 射击前全炮重心到驻锄支点水平距离为D0, 起落部分重量为Qq, 则由图可知
实际上, 由于l1不大, 射角变化后, 全炮重心的变化很小。
1.3.2 全炮重心随后坐行程的变化
如图6所示, G为发射前后坐部分质心位置, G′为对应后坐行程X的后坐部分质心位置, 则后坐行程为X时, 全炮重力对B点的后翻稳定力矩为
后坐行程为X时, 全炮重心对车轮与地面接触点的前翻稳定力矩为
1.3.3 力臂h随射角的变化
如图7所示, Hz为耳轴中心到地面垂直距离, △H为驻锄支点B到地面垂直距离, 由图可知,
1.3.4 力臂hf随射角的变化
如图8所示, L为车轮与地面接触点A到驻锄支点B的水平距离, 由图可知, hf=Hzcosφ+ (L-Lr) sinφ+e (15)
2 某型火炮射击稳定性数值仿真分析
基于建立的稳定性分析数学模型, 采用MATLAB语言编制程序[7], 可以得到某型火炮车轮着地射击时的稳定性时变特性 (表现为稳定力矩和翻转力矩随时间的变化曲线) 。
2.1 射角对稳定性的影响
射角20°射击时全炮驱动力Qr (后坐时为R, 复进时为r) 的变化如图9所示, 由于复进时r较小, 将其进一步放大。Qr为正时, 全炮受力如图2所示, 火炮有绕驻锄向后翻转的趋势;Qr为负时, 全炮受力如图4所示, 火炮有绕车轮向前翻转的趋势。图中椭圆虚线框内力的突然下降是由于后坐与复进的交替瞬间, 阻力突然换向所致 (参见式 (2) 和式 (4) ) 。复进过程中, 有两个时间段Qr为正, 两个时间段Qr为负, 这说明复进过程中, 火炮前翻后翻呈现一定的交替。
射角20°射击时, 稳定和翻转力矩如图10所示。力矩为正时, 表示全炮有后翻趋势, 力矩为负时, 表示全炮有前翻趋势。整个射击过程, 火炮2个阶段有后翻趋势, 2个阶段有前翻趋势, 与图9是对应的。由图10还可以看出, 不同阶段稳定力矩在数值上 (负数取绝对值) 均大于翻转力矩, 火炮能够保持稳定。但是不同阶段稳定程度变化很大, 譬如, 后翻趋势第二阶段的稳定性明显好于第一阶段。
射角0°射击时, 稳定和翻转力矩如图11所示, 可以看出, 后翻趋势阶段变为3个, 且在后翻趋势第一阶段, 翻转力矩大于稳定力矩, 火炮不能保持稳定。该炮的稳定极限角为13°, 故在射角0°射击, 必然发生车轮跳起的现象。
射角45°射击时, 稳定和翻转力矩如图12所示, 此时, 后翻趋势阶段有2个, 但是在这两个阶段, 翻转力矩均变为稳定力矩, 火炮不存在绕驻锄翻转的问题, 这实际上是由于图2中翻转驱动力的作用线从驻锄右侧移到了驻锄左侧所致。
2.2 火线高、驻锄长度及驻锄埋深变化对稳定性的影响
火线高、驻锄长度及驻锄埋深变化分别反映在参量Hz、Lr和△H上。以射角20°射击为参照 (图10) , 火线高降低20%、驻锄长度增加20%、驻锄埋深增加20%后的稳定力矩和翻转力矩分别如图13、图14、图15所示。可知, 火线高降低和驻锄长度增加能显著提高射击稳定性, 而驻锄埋深增加将导致射击稳定性有所下降。这些结论通常我们只是定性地了解, 这里通过仿真可以给出量化的描述。
2.3 后坐过程中Ppt·e对稳定性的影响
如式 (8) 所示, 后坐过程中, Ppt·e (不妨称之为附加矩) 对火炮稳定性有一定影响。Ppt数值很大, 但是e值很小, 附加矩Ppt·e对稳定性的影响大小却很难了解, 通过仿真获得后坐过程中的翻转力矩和去掉附加矩的翻转力矩如图16所示。由于该型火炮后坐部分质心在炮膛轴线上方, e为负值, 故附加矩将使翻转力矩变小。同时我们也发现, 附加矩对翻转力矩影响并不是很大, 或者说对火炮射击稳定性影响不大, 这主要是由于Ppt作用时间很短, 同时, 附加矩在翻转力矩中所占比例很小。图中考虑附加矩的翻转力矩在曲线顶部的短期增大是因为这个阶段为后效期, 该炮炮口制退器作用显著, 弹丸出膛后, Ppt方向改变。
3 结论
火炮射击稳定性是火炮轻型化必须考虑的核心特性, 通常我们对火炮射击稳定性的认识大都停留在定性的层面上, 本文通过仿真获得了某型火炮射击全过程稳定性的时变曲线, 为深入分析火炮射击稳定性提供了重要依据。
火炮射击过程中的稳定性表现出十分复杂的特征, 主要原因是:驱动力变化复杂、翻转方向存在交替, 稳定条件不断变化, 稳定力矩和翻转力矩随后坐行程、射角、火炮结构参数发生变化。这些原因导致其数学模型十分复杂。
对某型地面火炮车轮着地射击时稳定性仿真结果表明:火炮射击过程中稳定性不断变化 (表现为稳定力矩和翻转力矩的复杂变化) ;射角对稳定性有较大的影响, 火线高、驻锄长度、驻锄埋深对稳定性也有影响;Ppt·e对火炮稳定性影响不是很大。
摘要：为描述火炮射击全过程稳定性的时变特性, 基于稳定性分析数学模型, 采用数值仿真方法, 给出了某型火炮射击过程中稳定力矩和翻转力矩的时变曲线。仿真表明, 射击过程中全炮的稳定性呈现复杂的变化, 尤其在复进阶段, 火炮前翻和后翻趋势呈现交替变化的特征, 另外, 火炮射角对射击稳定性影响极大。
关键词：射击稳定性,数值仿真,数学模型
参考文献
[1]谈乐斌, 候保林, 陈卫民.降低火炮后坐力技术概述[J].火炮发射与控制学报, 2006 (4) :69-72.
[2]余子平, 钱林方, 徐亚东.基于动力学仿真的车载火炮射击稳定性优化[J].火炮发射与控制学报, 2006 (3) :36-40.
[3]宗士增, 钱林方, 何永等.轻型车载火炮二维后坐动力学研究[J].弹道学报, 2006, 23 (2) :76-78.
[4]高树滋, 陈运生, 张月林, 等.火炮反后坐装置设计[M].北京:国防工业出版社, 1995.
[5]张相炎, 郑建国, 扬军荣.火炮设计理论[M].北京:北京理工大学出版社, 2005.
[6]高跃飞.火炮反后坐装置设计[M].北京:国防工业出版社, 2010.
基于数值仿真的电价灵敏度分析 篇8
在过去的20多年中, 世界上许多国家开展了以“厂网分开、竞价上网”为主要特征的电力工业市场化改革, 电力企业从垄断型经营模式向竞争性经营模式转变, 电力市场成为国家经济格局中的重要组成部分。其中, 电价作为市场的核心问题, 直接关系到整个市场的发展和稳定。在电力市场环境中, 电价具有较高的波动性和随机性, 其正常波动是电力市场赖以生存和发展的基本条件, 对市场的健康发展具有积极的正面效应;反之, 电价的异常波动将加剧市场成员的运营风险, 并影响电力市场的稳定发展。尤其是2000年美国加州电价的暴涨给美国加州的经济和生活带来了巨大损失, 也引起人们对电价风险问题的普遍关注。然而, 由于电力市场的特殊性, 如何借鉴金融市场风险管理理论进行电价风险分析还有待研究。目前, 已有一些学者对电价问题展开了较为深入的研究, 并取得了丰硕的研究成果, 其研究主要集中在电价预测[1]、波动性分析[2,3,4]、概率分布[5,6,7,8]等方面。
总体而言, 灵敏度分析是常用的风险分析技术, 能够直观地表征自变量与因变量变化之间的内在联系, 这对于电价风险研究而言, 具有重要的理论意义和实用价值。究其原因, 电力市场因素的不确定性具有其客观必然性, 所以, 进行市场因素变化与电价变化之间的定性和定量分析, 以挖掘市场因素不确定性与电价不确定性之间的潜在规律, 实现对电价风险的深层次分析, 具有极其重要的意义。
目前的灵敏度分析方法已在金融资产市场风险分析中得到较为深入的研究和应用[9], 其研究方法大多是直接对目标变量与自变量之间的数学模型进行求导计算灵敏度系数。在电力系统分析领域, 灵敏度分析被广泛应用于系统仿真[10]、保护控制[11,12]、系统稳定性分析[13]等方面, 这些研究大多基于明确的数学模型或系统, 通过灵敏度分析实现对模型变量及参数的动态分析。然而在电价研究方面, 由于价格本身的复杂性和电力市场的特殊性, 对于电价及其影响因素的显示建模存在较大困难, 从而制约了解析法在电价灵敏度分析方面的应用。为解决该问题, 文献[14]提出了电量电价灵敏度矩阵的计算方法, 此方法所求的电价灵敏度是一种静态量度, 对于电价的动态灵敏度问题不适用。
本文从数值仿真的角度, 提出了一种新的电价灵敏度分析模型。该模型首先根据统计学习理论, 建立基于支持向量机的电价仿真模型, 并对电价影响因素进行小扰动模拟, 然后根据智能仿真结果, 计算电价的灵敏度。该模型将支持向量机的黑箱建模技术与灵敏度分析有机结合, 所建模型原理简单、适用性较强。
1模型原理
首先对电价灵敏度模型进行基本定义。设y (t) 为目标交易时段t的电价, x1 (t) , x2 (t) , …, xn (t) 为电价解释变量, 即模型计及的影响因素, 例如历史电价、负荷、网间交换功率、系统可调发电容量等, x (t) =[x1 (t) , x2 (t) , …, xk (t) ]T∈Rk, 且为避免影响因素数据级和单位差异所引起的灵敏度数量级上的差别, x已进行归一化处理。这里定义电价影响因素与电价之间的映射关系为:
简记为y=f (x) 。则电价y的全增量Δy可近似表达为:
式中:βi定义为相对于负荷、历史电价等不同影响因素的电价灵敏度, 且由于x为归一化数据, 所以灵敏度单位与电价单位一致。
进一步将Δy分解为:
于是, 构建电价灵敏度模型的步骤如下:首先, 建立基于支持向量机的电价仿真模型;然后, 利用该模型, 进行电价的数值仿真;最后, 根据仿真样本, 求解灵敏度系数。
2模型实现
针对动态电价灵敏度研究中存在的有限样本和价格本身的复杂性等问题, 本文提出基于支持向量机的黑箱建模技术, 通过对初始状态附近的各市场因素波动及电价响应的数值仿真, 研究影响因素变化与电价变化之间的动态关系, 实现电价灵敏度分析。
2.1建立基于支持向量机的电价回归模型
电力市场中电价与市场因素之间的回归建模属于复杂的高维非线性回归建模问题。在诸多的回归算法中, 支持向量机作为一项新的数据挖掘技术, 借助于最优化方法解决机器学习问题, 有效地克服了传统算法的“维数灾”和“过学习”等问题, 成为解决高维非线性回归问题的新工具, 已在时间序列预测[15,16]方面得到较广泛的应用。因此, 为保证后续灵敏度模型的准确性, 本文提出基于支持向量机的电价仿真模型。设对于任意交易时段t, 给定电价及其影响数据的训练样本集为:
式中:x∈Rk;y∈R。
根据统计学习理论, 模型采用非线性映射Φ将数据映射到高维特征空间, 然后在该空间内进行线性回归。设电价仿真模型为:
采用ε-不敏感函数时, 得到优化问题:
通过引入拉格朗日算子αi, α*i, γi≥0, 得到优化问题的对偶形式[17]:
记K (xi, xj) =ΦT (xi) Φ (xj) 。由于上述优化过程只考虑高维特征空间中的内积运算, 因此可采用核函数, 实现电价非线性回归建模, 建立电价仿真模型如下:
在此模型的基础上, 进一步利用测试集样本进行电价仿真模型的效验。具体地, 基于支持向量机的电价仿真模型实现流程如图1所示, 通过枚举迭代进行寻优, 在参数调整中考察训练误差和测试误差 (见式 (8) ) , 并选取测试精度最高的模型为电价仿真模型。根据该流程, 对各交易时段分别建立模型。
2.2数值仿真和电价灵敏度系数求解
该模块利用前述的电价仿真模型, 实现市场因素变化与电价变化之间的数值仿真。
设市场因素与电价的初始状态为[x (0) , y (0) ], 其中, x (0) =[x1 (0) , x2 (0) , …, xk (0) ]T∈Rk。对于任意的电价影响因素xi (i=1, 2, …, k) , 其第s次迭代新值为xi (s) , 且x
式中:j=1, 2, …, i-1, i+1, …, k。
根据上述样本集绘制电价波动与影响因素xi发生小扰动时的数据曲线, 其切线的斜率即为式 (3) 中电价对市场因素xi的灵敏度系数。以次类推, 建立最终的电价灵敏度模型。
3算例分析
本算例采用2002年10月—12月美国新英格兰现货电能量市场的负荷、网间交换功率和电价实际数据, 构建基于智能仿真的电价灵敏度模型, 并进行有效性验证。其中, 10月和11月的相关数据用于电价仿真模型的训练, 建立12月1日的电价灵敏度模型。
设y为t交易时段的系统边际价格, 电价影响因素x1和x2分别为同时段的负荷和网间交换功率, 则电价灵敏度模型形如Δy=β1Δx1+β2Δx2。首先对负荷和网间交换功率数据进行归一化处理, 根据本文的建模方法, 建立基于支持向量机的电价仿真方程, 其中, 核函数取线性核函数。在此基础上, 对负荷、网间交换功率及电价变化进行智能仿真, 仿真过程中负荷和网间交换功率均采用±3%基值以内之间的均匀波动, 输出为电价波动的仿真结果, 从而建立电价灵敏度模型。所建电价灵敏度模型的相关数据如图2和图3所示。图2为2002年12月1日6:00的电价灵敏度曲线。由图知, 电价与负荷、网间交换功率变化之间存在明显的线性相关性。于是, 求出该交易时段电价对负荷和网间交换功率的灵敏度系数分别为3.972 6美元/ (MW·h) 和-3.312美元/ (MW·h) 。由此发现负荷和网间交换功率对电价变化的不同贡献, 总体而言, 该交易时段电价对负荷和网间交换功率的灵敏度数值大致相同, 但不同因素引起的电价变化却刚好反向。
图3为2002年12月1日24个交易时段的电价灵敏度系数曲线。易见, 各交易时段电价对负荷和网间交换功率的灵敏度系数并不相同。但总体而言, 电价对负荷的灵敏度系数为正数, 说明电价变化与负荷变化方向相同, 即负荷增大时电价有增大的趋势;而在多数交易时段内, 电价对网间交换功率的灵敏度系数却为负数, 且绝对值略小, 说明多数时段内电价变化与网间交换功率变化方向相反, 即网间交换功率增大时电价有减小的趋势。
此外, 算例还将电价对负荷与网间交换功率的灵敏度系数进行了数值比较, 所得的灵敏度系数比率曲线如图4所示。可见, 在绝大多数交易时段, 电价对负荷的灵敏度系数与电价对网间交换功率灵敏度系数的比率都大于1。尤其是在9:00和12:00两个交易时段, 甚至超过10多倍。这说明在这些时段, 负荷对电价的影响占主要地位, 网间交换功率的影响占次要地位。例如12:00的电价灵敏度模型为:
当负荷和网间交换功率都变化3%时, 电价变化0.281 7美元/ (MW·h) , 其中负荷引起的电价变化为0.288 6美元/ (MW·h) , 网间交换功率引起的电价变化为-0.006 9美元/ (MW·h) , 易见2个影响因素对于电价变化贡献的较大差异, 这种波动量分解有助于更好地分析和理解电价波动。另一方面, 当比率远小于1时, 则说明此交易时段网间交换功率占主要地位, 负荷影响占次要地位, 例如在第4个交易时段。
从以上分析可见, 本文灵敏度分析模型原理简单, 易于理解和操作, 且建模过程可解释性强。
4结语
本文基于数值仿真的思想, 提出了一种新的电价灵敏度模型。该模型充分利用支持向量机的黑箱建模技术, 通过对市场因素波动的模拟, 考查不同市场情景下的电价动态响应。再根据激励与响应之间的对应关系, 成功求解出电价灵敏度系数。而且, 本文模型适用于电力市场下的小样本条件, 模型泛化能力较强。但是, 由于本灵敏度分析以回归模型为基础, 电价模型的回归精度会对灵敏度准确性有一定的影响。
数值仿真模型 篇9
进入21世纪后, 随着仪器设备和计算机技术的发展, 细观尺度下沥青混合料力学行为的相关研究逐渐成为热点。在构建细观力学数值模型过程中, 有限单元法将沥青混合料视为粗集料和砂浆两相[1], 离散单元法则划分为三类接触:集料内部、砂浆内部、集料与砂浆界面[2]。无论何种数值方法, 首先需获得沥青砂浆的宏观力学参数。通常静态参数如杨氏模量、抗压或弯曲强度等较易直接测得, 而由于试验操作难度、材料本身性质等原因, 砂浆动态参数通常采取对混合料加载过程进行反演获得。
本文借助基于离散单元法的三维颗粒流软件 (Particle Flow Code in Three Dimensional, PFC3D) , 利用室内试验与数值模拟相结合的方法, 对典型级配的沥青砂浆进行动态模量试验, 获取沥青砂浆的动态黏弹性参数, 研究材料高温黏弹性本构和动态力学行为, 同时为沥青混合料细观力学研究提供材料参数。室内试验与数值仿真的相互验证, 还能为后续使用虚拟试验部分代替真实试验提供借鉴。
1 室内试验与数值模型
1.1 沥青砂浆试件
为体现研究成果的普适性, 选取4种矿料级配 (AC-13、AC-20、SMA-10、SMA-13) 作为成型沥青砂浆的依据。文中采用玄武岩和SBS改性沥青作为原材料, 将粒径小于2.36 mm的细集料与沥青结合料的混合物视为砂浆, 根据沥青混合料级配换算出砂浆的各筛孔质量通过率, 即沥青砂浆级配, 如表1所示。
参考Superpave体积设计法2和《公路沥青路面施工技术规范》 (JTG F40—2004) 推荐公式, 假定矿粉与细集料表面均匀裹附着沥青, 颗粒表面积与其所裹附沥青质量为正比例关系, 按照比表面积法[5]获得沥青砂浆的油石比, 如表1所示。
沥青玛蹄脂混合料结构类型属于骨架密实型, 相对密级配混合料而言, 细料比例较少, 但油石比和矿粉含量较高, 因此通过换算得到砂浆的沥青用量偏高。可以认为, SMA系列成型的沥青砂浆的黏弹性能更为显著, 此外, 砂浆较高的沥青用量使得常规试验方法较难直接获取材料的高温黏弹性参数, 需根据时温等效原理予以确定。
文中动态模量试验试件尺寸为高150 mm、直径100 mm的圆柱体, 部分试件采用旋转压实 (SCG) 成型;部分试件在高温下呈现流动状态, 可直接倒入模具, 借助静压法成型, 待冷却后切割为目标尺寸。针对每类级配制备3个平行试件, 共4组12个试件进行测试, 各试件周边粘贴3组位移传感器, 如图1 (a) 所示。根据《公路工程沥青及沥青混合料试验规程》 (JTG E20—2011) 的T0738—2011章节, 通过Superpave简单性能试验机 (SPT) 开展不同温度 (15℃、20℃、25℃及30℃) 与频率 (0.1 Hz、0.5 Hz、1 Hz、5 Hz、10 Hz、25 Hz) 下沥青砂浆的单轴压缩动态模量试验, 如图1 (b) 所示, 对试件施加偏移正弦波轴向压应力, 控制应变模式 (85~110με) 。
1.2 数值模型
在PFC3D内编写程序生成规则排列的单元, 即每单元四周与6个单元相接触, 设定单元半径1 mm, 最终构建如图2 (a) 所示由148 200个离散单元组成的高150 mm、半径100 mm的圆柱体数字试件。文中将砂浆看作均质体, 因此所建立的数值模型由同样大小的单元组成, 并对每一单元赋予相同的微观特性, 相邻单元之间接触本构行为可以看作砂浆单元与砂浆单元之间的接触。
PFC3D软件提供的用户自定义Burger's接触模型如图2 (b) 所示, 由麦克斯韦尔 (Maxwell) 和开尔文 (Kelvin) 两部分模型串联组成, 包含切向和法向作用。Maxwell模型和Kelvin模型被广泛应用在沥青混合料黏弹性分析中, 分别侧重于描述材料的蠕变和应力松弛行为。而沥青结合料/砂浆/混合料的黏弹性能往往呈现出蠕变和松弛的叠合效应, 因此, 综合了Maxwell和Kelvin模型特点的Burger's模型, 结构简单明了、参数数目适宜, 是研究人员最为常用的经典黏弹性模型。
离散单元法通过微观力学响应表征材料宏观力学特征, 室内试验测得的材料参数实际是细观力学行为的宏观反映, 不能直接作为离散元模型的输入参数, 有关Burger's模型中宏观力学属性与微观力学元件之间联系, Liu[6]已在其论文中给出详细推导过程, 文中直接引用, 图2 (b) 所示宏观本构行为与法、切向微观接触行为之间的关系为
式中:E1、E2分别为Maxwell和Kelvin模型内弹簧劲度;η1、η2分别为Maxwell和Kelvin模型内黏壶黏度;L为相邻单元球心距 (文中即为单元直径) , v为沥青砂浆泊松比, μ为摩擦系数;分别为Burger's接触模型法向微观特性;分别为Burger's接触模型切向微观特性。
综上所述, 通过沥青砂浆室内单轴蠕变试验获取宏观Burger's四参数, 借助公式 (1) 和 (2) 即可得到离散元模型的微观特性。
2 结果与分析
2.1 动态模量
PFC3D软件提供了“WALL”命令, 通过控制“WALL”单元的速度实现对模型的加载。在如图2 (a) 所示模型顶底面各生成直径100 mm的圆形“WALL”单元, 固定底面自由度, 不断调整顶面“WALL”单元的移动速度实现正弦波加载。在数值模拟过程中, 实时采集模型顶底面“WALL”单元所受的不平衡反力, 取均值后与作用面积相除即可获得轴向应力;通过实时监控“WALL”单元位移, 与模型原始高度相除即可获得轴向应变。
以AC-13型砂浆为例, 15℃条件下, 0.1 Hz、1 Hz、5 Hz和10 Hz下其应力应变如图3所示。
离散元模拟与室内试验获得的应力应变曲线发展趋势基本一致, 由模拟结果可以看出, 荷载应力按照正弦波形式施加于数值模型;随着频率的增大, 轴向应力幅值逐渐提高。由图3可见轴向应变幅值随着荷载作用时间的延长而呈现出增大趋势 (0.1 Hz低频作用下尚未体现) , 表明与沥青混合料相比, 沥青砂浆具有更为显著的黏弹特性, 在重复荷载作用下的蠕变变形累积效应更为明显。此外, 相比AC型砂浆, SMA型砂浆在相同频率作用下反映出的轴向正弦应力偏小, 此现象也与之前提及SMA型砂浆具有更明显的黏弹特性这一认识相符。
沥青砂浆的动态模量E*和相位角Φ可以通过动态模量室内试验或离散元模拟得到的应力应变曲线并结合式 (3) 求得:
式中:σmax、σmin、εmax和εmin分别为加载应力、应变响应的峰值和谷值;Δt为相邻峰值应力和峰值应变的时间差;T为加载周期, 加载频率的倒数。
基于图3所示应力应变曲线, 按照式 (3) 计算得到离散元模拟的沥青砂浆动态模量和相位角。对15℃、25℃条件下5个加载频率的数据进行比较, 结果如图4所示。
由图4可以看出相同温度与加载频率下, 无论数值模拟还是室内试验, 最大粒径较大的AC-20型砂浆动态模量均略高于AC-13型, 而最大粒径较小的SMA-10型砂浆动态模量却高于SMA-13型;由表1数据可知, SMA-10型砂浆中含有的0.3~0.6 mm、0.6~1.18 mm和1.18~2.36 mm这3档细料分别高于SMA-13型6.5%、5.5%和2.5%, 同样AC-20型砂浆在0.3~2.36 mm区间内的细料含量也略高于AC-13型;此外, AC-20型和SMA-10型砂浆的沥青用量也分别高于AC-13型和SMA-13型。分析认为, 细集料含量的增加仍然提高了沥青砂浆的抵抗变形能力, 而过高的沥青用量增大了砂浆的蠕变变形。
由图4还可以看出, 各温度、不同频率以及不同类型沥青砂浆条件下, 离散元模拟得到的动态模量均高于室内试验结果;分析认为, 虽然砂浆是由粒径2.36mm以下细集料与沥青胶结料组成, 但仍是非均质混合物, 强度存在空间分布的随机性, 而文中数值模型将砂浆作为均质结构考虑, 每一离散单元均赋予相同的微观特性, 是一种理想状态, 因而导致了模拟结果与室内试验的误差, 最大13%, 平均8%。此外, 不同温度条件下, SMA型砂浆模拟与室内结果的匹配程度 (最大误差8%, 最小4.8%, 平均6.2%) 要优于AC型砂浆 (最大误差13%, 最小7.4%, 平均9.5%) , 说明文中引入的Burger's接触本构模型适用于高油石比的间断级配类型沥青砂浆时, 能够更好地发挥本构模型的黏弹性特征。
在未开展材料微观特性试算校核的前提下, 文中离散元模拟结果与室内试验结果匹配度较好, 说明文中建立的数值模型、引入的Burger's接触本构模型以及离散元方法, 能够较好地实现具有黏弹特性的沥青砂浆动态模量数值仿真, 通过不断调整砂浆的微观特性, 可以取得与室内结果更优的匹配度。
2.2主曲线建立
参考AASHTO Design-PP62所述方法, 选用包含于力学经验法的S型函数, 选用AASHTO规范中的二次多项式移位因子。将移位因子方程带入S型函数 (即动态模量主曲线通用方程) , 得到最终形式如下:
式中:α, β, δ, γ, a1, a2均为拟合参数;f为试验温度下加载频率;TR和T分别为参考温度和试验温度。
限于篇幅, 文中选取AC-13和SMA-13型砂浆的室内试验数据, 以20℃为主曲线参考温度, 利用Origin8.5软件根据式 (4) 进行拟合, 相关参数如表2所示。2种不同级配类型的沥青砂浆在各温度条件下的动态模量曲线, 按照移位因子变换后如图5所示。需要说明的是, 在室内试验过程中, 由于沥青用量较高, SMA类型沥青砂浆在温度超出25℃即呈现软化现象, 因此在研究时不考虑30℃。
由图5可以看出, 由于荷载频率范围较小, 沥青砂浆动态模量主曲线没有呈现出与沥青混合料动态模量主曲线形式相似的S形分布, 而是在一定程度上表现出均匀的变化状态。可以认为, 沥青砂浆在结构上确是类似均质特征, 所体现出的宏观力学行为较为稳定, 也为文中在数值模拟时假定砂浆为均质材料提供了支持;而沥青混合料则为典型的非均质结构, 即使处于相同荷载与温度条件下同样呈现出力学行为差异。
3结论
(1) 离散元模拟实现了正弦波加载, 砂浆轴向应变幅值随荷载作用呈现增大趋势, 表明沥青砂浆与混合料相比具有更显著的黏弹特性。
(2) 细集料含量较多的砂浆动态模量较大;离散元模拟的动态模量略高于室内结果, 在不调整材料微观特性前提下, 总体匹配较好;油石比较大的SMA类型砂浆模拟结果与室内试验误差更小。
(3) 沥青砂浆动态模量主曲线一定程度上呈现出均匀变化状态, 而非沥青混合料通常的S形分布, 体现了砂浆的类似均质材料特性;主曲线的建立, 可以突破试验条件限制, 得到沥青砂浆在全频范围内动态模量数值和变化规律。
(4) 引入Burger's接触本构的离散元方法, 能够较好地对重复荷载作用下沥青砂浆的动态模量进行仿真, 并可以进一步推广应用在疲劳试验仿真中。
参考文献
[1]Wang Jiangyang, Qian Zhendong.Indirect tension test of epoxy asphalt mixture using microstructural finite-element model[J].Journal of Southeast University (English Edition) , 2011, 27 (1) :65-69.
[2]王江洋, 钱振东, Wang Lin-bing.环氧沥青混凝土裂纹发展的三维可视化离散元模拟[J].湖南大学学报, 2014, 41 (6) :112-119.
[3]谭凤吉.CRTS II型水泥乳化沥青砂浆灌注质量通病及防治措施[J].铁道建筑, 2012 (4) :139-142.
[4]徐浩, 王平, 徐井芒, 等.反复荷载下水泥乳化沥青砂浆力学性能试验[J].建筑材料学报, 2014, 17 (4) :715-720.
[5]汪海年.沥青混凝土微细观结构及其数值仿真研究[D].西安:长安大学, 2007.