Cosserat理论的有限元实现及微梁弯曲的尺寸效应模拟
两种强度准则在Cosserat扩展模型弯曲变形中的应用_
2 二维Cosserat扩展模型
2.1 Cosserat 理论基本方程 Cosserat 理论是 Cosserat 兄弟[19]于 1909 年首先 提出，也称偶应力理论。其研究对象是具有一定特 征尺度、含有微结构的微元体；它除受常规应力作 用外，还考虑了偶应力的影响。 现取 x1-x2-x3 直角坐标系，如图 1 所示。其中， x3 轴正向为垂直于平面向外。 相较于经典弹性理论， Cosserat 介质理论在平面问题中由于多考虑了偶应 力作用，故将原先的应力、应变分量扩展为广义的 应力、应变，定义如下：
第 37 卷第 1 期 2016 年 1 月
DOI：10.16285/j.rsm.2016.01.033
岩 土 力 学 Rock and Soil Mechanics
Vol.37 No.1 Jan. 2016
两种强度准则在 Cosserat 扩展模型 弯曲变形中的应用
张建成 1，贾金青 1，涂兵雄 2，何慧荣 3
形屈服面是不光滑且有尖角，这些尖角可能会导致其应用于 Cosserat 弹塑性分析时的计算困难，Drucker-Prager 强度准则正 可被看成 Mohr-Coulomb 强度准则为避免这些困难而做的光滑近似， 因此很有必要对基于 Drucker-Prager 强度准则的 Cosserat 理论进行研究。 利用 MATLAB 分别编写了基于 Mohr-Coulomb 和 Drucker-Prager 强度准则的 Cosserat 扩展模型的有限元程序， 并对在互层岩体中开挖的洞室进行了变形分析。结果表明，两种强度准则均可用于 Cosserat 扩展等效模型，但基于 Drucker-Prager 强度准则的有限元程序收敛速度更快，稳定性更好，能得到更为理想的结果。 关 键 词：Cosserat 扩展模型；弯曲变形；Mohr-Coulomb 准则；Drucker-Prager 准则；MATLAB；弹塑性有限元 文献识别码：A 文章编号：1000－7598 (2016)01－0279－08 中图分类号：TU 473.1
平纹织物三点梁弯曲有限元模拟与分析
平纹织物三点梁弯曲有限元模拟与分析靳欢欢;杜赵群【摘要】基于纱线与织物性能原位综合测试系统(CHES-FY),采用有限元软件ABAQUS构造平纹织物的三点梁弯曲模型,并实施弯曲模拟分析.通过模拟计算得到了织物中纱线应力分布和压针上抗弯力-位移曲线,提取弯曲特征指标(抗弯力、抗弯功和抗弯斜率),与试验测得的抗弯力-位移曲线对比分析,并讨论了纱线模量、托针间距和纱线摩擦因数对弯曲性能特征指标的影响.结果表明:织物弯曲时,经纱方向受力较大,纬纱方向受力较小;压针上的抗弯力-位移模拟曲线和实测曲线趋势相同,呈先上升到最高点、后下降的趋势;随着纱线模量的增加,托针间距增加,纱线摩擦因数增加,压针上抗弯力、抗弯功和抗弯斜率也随之增加.模拟分析与试验结果一致,从而证明了有限元方法模拟织物弯曲的可行性.【期刊名称】《东华大学学报（自然科学版）》【年(卷),期】2016(042)003【总页数】6页(P344-349)【关键词】弯曲;平纹织物;有限元模拟;抗弯力-位移曲线【作者】靳欢欢;杜赵群【作者单位】东华大学纺织面料技术教育部重点实验室,上海201620;东华大学纺织面料技术教育部重点实验室,上海201620【正文语种】中文【中图分类】TS116织物弯曲性能为抵抗弯曲形变的能力，既决定着外套用布料的硬挺和保形性能，又决定着内衣用布料的柔软度和悬垂特性，故研究快速客观测量织物弯曲性能的方法极为必要.目前，最为经典的测试方法为斜面法、纯弯曲法和三点梁弯曲法.斜面法是将织物试样由水平面推出至触及与水平面成规定角度的斜面时的推出长度，获得织物弯曲长度和弯曲刚度来评价织物的弯曲性能[1]，该方法简易、快速，适宜纺织技术开发人员快速评价织物的弯曲硬挺性.澳大利亚联邦科学院的SiroFAST系统采用此法测量织物的弯曲性能[2].纯弯曲法是将织物试样竖直放置，以垂直于竖直方向的弯矩往复弯曲试样获得弯曲刚度表征织物弯曲性能，可避免斜面法中重力的影响.日本计量与标准化委员会研发的KES-FB2弯曲仪采用纯弯曲法测试低应力下织物的弯曲性能[3].三点梁弯曲法则是将织物试样水平放置在一对水平托针上，由位于一对水平托针中间的压针弯曲试样获得织物弯曲力，据此获得织物的弯曲刚度来表征织物的弯曲性能.东华大学研制的纱线与织物性能原位综合测试系统(CHES-FY)就是基于三点梁弯曲原理实施织物的弯曲性能测试，其可模拟织物手感风格评价中手指与织物的交互作用[4-5].相对而言，斜面法和纯弯曲法测试边界难以模拟手指的触摸行为，三点梁弯曲法则可以近似获得织物在手指间的摩擦交互作用. 为了明晰三点梁弯曲法测试过程中交互作用行为，精准测量织物的弯曲性能，需要进行系统参数的理论研究.基于有限元方法的纺织材料结构性能研究愈加普遍，文献[6-8]采用不同的方法建立织物模型；文献[9]模拟不同层叠结构芳纶纤维织物增强复合材料的弯曲试验，得到了试样的弯曲应力分布和弯曲模量及弯曲刚度；文献[10]计算三维多重纬经角连锁玻璃纤维织物准静态三点弯曲破坏过程，比较试验和有限元计算得到的载荷-位移曲线和破坏形态；文献[11]模拟机织物的顶破过程，分析顶破过程中机织物的应力-应变分布、能量吸收机制以及摩擦作用的影响；文献[12]分析芳纶纱线和织物的拉伸断裂性能等.与试验测试和力学分析相比，有限元方法可更加直观地展示弯曲过程中的受力变形状态，研究织物弯曲过程中应力分布，尤其是仪器系统参数对弯曲测试指标的影响.因此，本文基于ABAQUS有限元分析软件实施平纹织物三点梁弯曲模拟，分析三点梁和织物的结构性能参数对织物弯曲性能测试结果的影响.1.1 织物试样规格本文采用平纹棉织物，其规格参数如表1所示，将织物裁剪成30 mm×20 m m 的长方形.抽取织物中纱线对其进行拉伸性能测试，其力学性能如表2所示.1.2 织物三点梁弯曲试验所用试样试验前均在温度(20±2)℃和相对湿度(65±3)%条件下平衡24h，放置备用.采用CHES-FY中的三点梁机构进行测试，测试原理如图1所示.测试时织物放在两托针上，两端不固定，驱动机构控制压针向下运动，测试速度为0.2mm/s.与压针相连的力学传感器记录所受的力和压针移动的位移.通过扫描电子显微镜(SEM)获得平纹织物试样的横截面图(如图2所示)，据此建立织物中经纱和纬纱的部件模型.建立模型前做如下假设：(1)忽略纱线是由纤维组成的，不考虑纱线内部空隙，纱线是实体单元;(2)纱线屈曲轨迹为正弦曲线，经、纬纱线截面为椭圆形;(3)纱线为各向同性材料;(4)织物经纱和纬纱的屈曲相同，截面形状也相同.由扫描电子显微镜图得到纱线厚度和宽度分别为0.221 和 1.314mm，纱线屈曲波长(即单胞尺寸) 为1.728mm，所得的平纹织物的有限元模型如图3所示.平纹织物采用八节点线性积分缩减积分单元(C3D8R)进行网格划分，其中，纱线厚度和宽度方向分别划分为2个和6个单元，单胞长度方向划分为50个单元.根据织物中纱线的屈曲状态，假设经、纬纱弯曲形变过程中整体发生变形，将纱线看成实体；弯曲形变过程中的经纱和纬纱的摩擦因数、泊松比和密度变化等均忽略不计，其中，纱线弹性模量为279 MPa，密度为0.85 g/cm3，泊松比为0.3，摩擦因数为0.3.3.1 织物弯曲过程中形态及应力分布设定平纹织物的基本尺寸为30.032 mm×19.664 mm，经纱与x轴方向平行，纬纱与z轴方向平行.由23根经纱和35根纬纱构成织物平面为xOy面，法线为y轴，托针和压针设为解析刚体，设定托针和压针直径均为2mm，长度均为30mm，压针沿y轴方向向下运动弯曲试样.建立如图4所示的织物三点弯曲模型，对建立的有限元模型进行约束和加载，织物两端为自由端，托针1和2固定，压针只有y 轴方向上的平移自由度，织物弯曲模拟过程如图5所示.由图5可知，经纱受力较大，纬纱受力较小.这是由于经纱垂直于压针运动方向，纬纱平行于压针运动方向，纱线的交织和摩擦相互作用使经纱方向受力较大，主要在于经纱抵抗弯曲变形.压针与织物接触，经向纱线先拉伸受力，随着压针继续弯曲织物，由于纱线的交织和摩擦相互作用，应力以压针为轴对称中心向两侧呈梯度递减传递.为了对比模拟与试验测试的差异，提取压针上的抗弯力，获得弯曲过程中的平纹织物抗弯力-位移曲线，其模拟与试验测试曲线如图6所示.由图6可知，有限元模拟和试验得到的织物抗弯力-位移曲线趋势相同，随着弯曲程度的增加，抗弯力先上升达到最大抗弯力点，后下降.初始阶段织物抗弯力呈线性上升，因为压针与织物接触，织物受压产生反作用力，作用于针上的力随压针位移的增大而增大，之后织物开始在托针表面滑移.模拟和试验曲线在达到峰值前都比较光滑，且近似线性增加.有限元得到的模拟结果与试验结果存在差异，这是由于有限元模拟将纱线看成均匀实体，忽略了纱线弯曲过程中纤维之间的交互作用力.为了比较试验测试结果与有限元模拟结果的一致性，从图6中提取最大抗弯力、抗弯功(曲线上力值从0到峰值与位移围成的面积)和抗弯斜率特征参数，如表3所示.由表3可知，当托针间距为8 mm时，有限元模拟曲线和试验曲线提取的3项特征参数中，其中两项的相对误差比较小，由此可以得出有限元模拟结果和试验结果具有一定的一致性.3.2 托针间距对抗弯力-位移曲线的影响设定纱线模量为279 MPa、摩擦因数和泊松比均为0.3、压针直径为2 mm时，有限元模拟托针间距分别为6，8和10 mm时平纹织物抗弯力-位移曲线图，如图7所示，不同托针间距下的特征参数结果如表4所示.由图7可知，不同托针间距时抗弯力-位移曲线趋势相同，都是先上升后下降，并随着托针间距的增大，曲线由陡峭变得平缓.由表4可知，随着托针间距的增加，最大抗弯力、抗弯功和抗弯斜率均减小.压针直径不变时，托针间距大，同等外力作用于织物时，织物与支撑两端的弯矩大，织物相对容易弯曲，表现为抗弯力小；托针间距减小时，抗弯斜率增加，故压针下压时织物受到的抗弯力增加，最大抗弯力增加，抵抗弯曲变形的抗弯功增加.3.3 纱线模量对抗弯力-位移曲线的影响设定托针间距为10 mm、纱线泊松比和摩擦因数均为0.3时，采用有限元模拟纱线模量(E)分别为44.74，447.4和4 474 MPa时对应的平纹织物抗弯力-位移曲线图，如图8所示，提取对应的最大抗弯力、抗弯功和抗弯斜率结果如表5所示.由图8可知，纱线模量不同时对应的抗弯力-位移曲线变化趋势一致，均随着弯曲程度的增加，抗弯力先增加后减小.由表5可知，纱线模量增加10倍，最大抗弯力也增加10倍左右，且抗弯功和抗弯斜率与纱线模量成比例增加，趋势相同.这表明纱线模量越大，即抗弯刚度越大，织物弯曲抗弯力越大，这与理论相一致.3.4 摩擦因数对抗弯力-位移曲线的影响设定纱线模量为279 MPa、泊松比为0.3、托针间距为10 mm时，摩擦因数分别为0，0.1和0.3的平纹织物抗弯力-位移曲线图，如图9所示，提取对应的最大抗弯力、抗弯功和抗弯斜率结果如表6所示.由图9可知，当纱线的摩擦因数增加时，纱线间的相互作用增加，织物弯曲过程中所受到的抗弯阻力增加.由表6可知，摩擦因数从0增加到0.1时，最大抗弯力、抗弯功和抗弯斜率明显增加.当摩擦因数从0.1增加到0.3时，最大抗弯力和抗弯功变化不明显，但抗弯斜率增加.从模拟结果可知，织物试样与三点梁结构之间的摩擦交互作用，也影响织物的抗弯曲变形，即抗弯力-位移曲线形态.总体而言，随着摩擦因数的增加，抗弯斜率增加，最大抗弯力和抗弯功变化相对平稳，有限元模拟结果与理论基本一致，这从而证明了有限元方法模拟织物弯曲的可行性.本文采用有限元软件ABAQUS模拟平纹织物的三点梁弯曲，由织物弯曲性能有限元分析和试验测试的抗弯力-位移曲线对比可知，抗弯斜率、最大抗弯力和抗弯功相对误差较小，有限元模拟和试验结果具有较好的一致性.另外，讨论了纱线模量、托针间距和纱线摩擦因数对机织物三点梁弯曲性能的影响.纱线模量对织物弯曲性能影响较大，随着纱线模量的增加，有限元模拟的抗弯力成比例增加，纱线模量越大，织物抗弯力值越大；随着托针间距的减小，织物抗弯力增加较快；随着纱线摩擦因数的增加，织物抗弯斜率增加，最大抗弯力和抗弯功变化较小.【相关文献】[1] PEIRCE F T. The handle of cloth as a measurable quantity [J]. Journal of the Textile Institute, 1930, 21(9): 377-416.[2] POSTLE R. Fabric objective measurement, part VI: Product development and implementation [J]. Textile Asia, 1989, 20(10): 59-68.[3] KAWABATA S. The standardization and analysis of hand evaluation [M]. 2nd ed. Osaka: The Textile Machinery Society of Japan, 1980.[4] DU Z Q, YU W D. Physical interpretation of pulling-out curve based on a new apparatus [J]. Journal of the Textile Institute, 2008, 99(5): 399-406.[5] DU Z Q, SHEN H, ZHOU T, et al. Comparison of properties characterization between CHES-FY, KES-F and FAST[J]. Industrial Textila, 2011, 62(3): 123-128.[6] HUA L N, MIKE J C, ANDREW C L. A finite element approach to the modelling of fabric mechanics and its application to virtual fabric design and testing[J]. Journal of the Textile Institute,2012,103(10): 1063-1076.[7] MASOUMEH V, STEPAN L. Finite element simulation of a yarn pullout test for plainwoven fabrics[J]. Textile Research Journal,2010, 80(10),892-903.[8] HASSEN H, ADEL G, HEDI B. Dynamic fabric modelling and simulation using deformable models[J]. Journal of the Textile Institute,2011,102(8):647-667.[9] 寇兴才，李旭东，张斌，等.芳纶纤维复合材料弯曲性能试验仿真[J].塑料，2013,12(4)：78-81.[10] 周亮，唐予远，聂建斌.多胞织物结构复合材料静态弯曲性能与有限元分析[J].天津工业大学学报，2008,27(5):4-7.[11] 张天阳，丁辛.机织物顶破过程的有限元模拟[J].东华大学学报(自然科学版)，2012,38(12)：688-694.[12] 程建芳，肖露，柴晓明，等.有限元分析法研究Kevlar纱线及织物的拉伸性能[J].浙江理工大学学报，2013，30(5)：649-653.。
正畸不锈钢弓丝弯曲成型与回弹的有限元数值模拟研究
正畸不锈钢弓丝弯曲成型与回弹的有限元数值模拟研究东南大学硕士学位论文正畸不锈钢弓丝弯曲成型与回弹的有限元数值模拟研究姓名孔灿申请学位级别硕士专业机械制造及其自动化指导教师汤文成20090301摘要正畸不锈钢弓丝弯曲成型与回弹的有限元数值模拟研究研究生姓名孔灿导师汤文成教授学校名称东南大学摘要弯曲是成型工序中最为普遍的成型方法之一其应用相当广泛。
在口腔正畸临床中也需要对提供正畸力的不锈钢丝弯曲成型。
传统的弯制过程在操作上缺乏定量指标弯曲出来的钢丝与指定形状匹配度较低不仅在弯制过程中耗费了大量的时间而且由于反复的弯曲容易引起正畸丝的疲劳断裂现象。
本文运用有限元中的优化方法针对不锈钢正畸丝弯曲成型和回弹问题进行了研究得出了恰当的弯曲形状使正畸丝在弯曲回弹成型后恰好达剑临床应用的要求。
通过理论解析方法验证了有限元计算结果的可靠性以后设计出了一种正畸不锈钢丝弯制器。
使用该弯制器可以方便、快速的将这种正畸不锈钢丝弯曲成型实现了临床中正畸丝弯制过程中操作的标准化。
从而可以减少正畸丝弯制操作的人为性提高正畸丝弯制过程的效率并最人程度上减少患者的痛苦避免医疗事故的发生。
在研究过程中得出了以下结论。
实际不锈钢丝弯曲加工中变形速度比较慢弯曲的速度比较恒定在模拟过程中可以采用静态分析。
使用静态分析模拟的精度完全可以满足临床医学中的要求。
使用无模法进行仿真可以将复杂的接触等非线性转换为力或位移的加载问题从而在保持精确度的同时大大节省计算的时间。
使用参数化语言进行仿真可以通过改变其中的参数实现同一类别的模拟。
本文所使用的语言经过简单修改后即可计算其它型号止畸丝的弯曲同弹问题。
本文的研究展示了在口腔正畸生物力学研究领域中采用有限元方法和数值模拟方法计算正畸不锈钢丝的弯曲成型问题可以为临床正畸学提供较人的方便。
而且将本研究扩展开来可以适用于不同类型的类形工件的弯曲同弹问题为其它复杂零件的弯曲成型加工时的回弹问题提供分析和解决的方法。
基于COSMOSWorks的电弧炉钢包车梁分析
The Analysis of the Ladle Car Beam for EAF Based on COSMOSWorks Software
LEI Pei，WANG Xu-peng
工业炉作为主要的炼钢设备在冶金行业得到广泛的 应用。钢包车梁作为其重要的部分由于结构复杂、零件 数量较多，特别是其涉及强度的校核，给设计带来诸多 不便 。然而，目前工业炉产品的设计主要以二维 CAD 软件为主， 其强度计算只能通过力学的方法等效为简支 梁进行人工计算，整个计算过程比较繁琐，且计算结果 准确性无法保证。本文以 COSMOSWorks 有限元分析 软件对工业炉钢包车梁进行分析， 不仅大大简化了分析 过程，同时通过各类图表、梯度表、曲线图直观地表达 了应力、应变、设计安全系数等结果，为结构设计优化 提供了依据。
软件简介
COSMOSWorks 是一种基于有限元分析技术 （即 FEA 数值技术） 的设计分析软件， 是众多 CAE 软件中的一种。 凭借其先进的快速有限元技术（FFE） , 使工程师能非 常迅速地运用它实现对大规模复杂设计的分析和验证, 获 得修正和优化设计所需的必要信息；同时与 SolidWorks 无缝集成，可将分析结果直观地列入图形化的树形分析 管理器中。此外，COSMOSWorks 有不同的包和应用程 序，以适合不同用户的需求。 本文所使用的 COSMOSWorks Professional 是 COSMOSWorks 模块中的一种，能够进行零件或装配件的静 态、热传导、扭曲、频率、掉落测试、优化及疲劳分析， 且操作简便, 结果准确可靠 。
ZnO （s） ＋ Coal ZnFe 2 O（s） Coke Breeze 4 挥发 碳热还原 Zn l） （ △ △ 精炼 Zn （粗锌） （s）
横观各向同性岩土材料应变局部化现象的有限元模拟
横观各向同性岩土材料应变局部化现象的有限元模拟常江芳;徐远杰;楚锡华【摘要】The anisotropic properties of geomaterials are significantly related to their inherent microstructures.In this paper,a modified Drucker-Prager yield criterion for transversely isotropic materials is developed by evaluating the internal friction angle with the stress state,the microstructure tensor and the material principal direction.Then,based on the Cosserat continuum theory,a consistent return mapping algorithm for the modified criterion is formulated,and a consistent tangent modulus matrix is achieved.Moreover,the codes are implemented through the user defined element subroutine (UEL) in the finite element software Abaqus,and the correctness of the program is verified by comparing the theoretical results with the relationships of the material strength to the principal direction and anisotropic degree of the material in the integration points.Finally,the influences of the principal direction and anisotropic degree of the material on the bearing capacity and the failure mode of the structure are emphatically analyzed by numerical examples,which are then compared with the results based on the classical continuum theory.It is found that the above-mentioned method is effective in simulating the strain localization of transversely isotropic geomaterials.%岩土材料的各向异性性质与其微结构紧密相关,文中考虑到内摩擦角是应力状态、组构张量及材料主方向的函数,发展了适用于横观各向同性岩土材料的修正的Drucker-Prager屈服准则.基于Cosserat连续体理论,推导了该准则的一致性映射返回算法,形成了一致性切线模量矩阵,并利用有限元软件Abaqus的用户单元子程序(UEL)进行了数值实现.通过将积分点材料强度随材料主方向及各向异性程度的变化关系与理论结果进行比较,验证了程序开发的正确性.数值算例重点分析了材料主方向和各向异性参数对结构极限承载力及破坏模式的影响,且与经典连续体的结果进行了比较.结果表明,文中方法能够较好地模拟具有横观各向同性的岩土材料的应变局部化现象.【期刊名称】《华南理工大学学报（自然科学版）》【年(卷),期】2016(044)010【总页数】11页(P70-80)【关键词】岩土材料;岩土力学;横观各向同性;Cosserat连续体;屈服准则;应变局部化【作者】常江芳;徐远杰;楚锡华【作者单位】武汉大学工程力学系,湖北武汉430072;石家庄铁道大学工二程力学系,河北石家庄050043;武汉大学工程力学系,湖北武汉430072;武汉大学工程力学系,湖北武汉430072【正文语种】中文【中图分类】O34自然界的岩土材料一般均呈现很强的各向异性.各向异性又分为固有各向异性和诱导各向异性.由于其沉积过程、孔隙、裂纹等固有属性而使其在原位状态下表现出的各向异性称为固有各向异性；岩土材料在非比例加载或塑性变形的影响下，由于主应力方向不断变化，导致其细观结构发生演化而表现出的各向异性称为诱导各向异性.关于此两种各向异性，国内外学者做了相关的理论及实验研究[1- 3]，提出了多种各向异性屈服准则[4- 7].姚仰平等[4]以横观各向同性土为例，阐述了同时考虑材料外部应力分布与材料内部强度分布下各向异性土的破坏机制.贾乃文[5]以Hill正交异性屈服准则为基础，讨论了横观各向同性条件下，岩土中出现塑性屈服的Drucker准则修正表达式.Duveau等[6] 对一些强各向异性材料的屈服准则做了一个综合评价.Gao等[7]在屈服函数的摩擦系数中引入了一个各向异性变量，考虑它为应力不变量及微结构张量的联合不变量的函数，提出了一个适用于横观各向同性材料的屈服准则.然而基于唯象方法提出的各向异性屈服准则一般包含较多的材料参数，且这些参数与材料微观结构之间的关系不明确.Pietruszczak等[8- 9]提出了一个各向异性屈服准则，基于物质点给出了各向异性程度及加载方向对材料强度影响的理论分析，其材料参数是应力和微结构张量联合不变量的显式函数，使得实验验证更易实现.余天堂等[10- 11]基于Pietruszczak等[8- 9]提出的各向异性屈服准则，考虑各向异性参数和单轴抗压强度是一个由微结构张量和加载方向表示的分布函数，给出了沉积岩的一种各向异性模型.Zhong等[12]对Pietruszczak-Morz 屈服准则中的材料参数进行了分析.在Pietruszczak等[8- 9]的研究基础上，Lade[13- 15]推导了一个基于 Lade 强度准则的各向异性强度准则，并展开了相应的理论及实验研究.然而上述研究多关注的是对各向异性破坏准则的描述，是各向异性对一个物质点破坏的影响.关于各向异性对结构承载力及破坏模式影响的数值模拟仍相对较少.Chu 等[16]在Pietruszczak等[8- 9]提出的各向异性屈服准则的基础上，考虑内摩擦角为组构张量、应力状态及材料主方向的函数，对Drucker-Prager屈服准则进行了修正，以横观各向同性材料为例研究了各向异性对应变局部化的影响，然而不足之处是其基于经典连续体理论，在伴随应变软化出现应变局部化现象的同时会产生网格依赖性.为克服这一问题，文中基于Cosserat连续体理论[17- 19]，推导了修正的Drucker-Prager屈服准则的一致性映射返回算法及切线模量矩阵，重点分析了材料主方向和各向异性程度对横观各向同性材料极限承载力及破坏模式的影响，最后通过与经典连续体计算结果进行比较来验证解决网格依赖性问题的有效性.1 基于Cosserat连续体的横观各向同性弹塑性模型1.1 Cosserat连续体的基本方程与经典连续体相比，Cosserat连续体中引入了独立的转动自由度，平面应变情况下，其位移矢量u可表示为[17]u=[ux uy ωz]T其中，ux、uy为平动位移，ωz为独立的转动位移.应变矢量ε除了与经典意义上所对应的正应变和剪应变有关以外，还包含两个微曲率，如下所示：ε=[εx εy εz εxy εyx lcκxz lcκyz]T其中：可以看到这里引入了一个内部长度参数lc，保证了Cosserat连续体的控制方程在计算过程中的正定性，可有效克服网格依赖性问题.几何方程可表示如下：ε=Lu其中，Cosserat连续体的应力矢量σ引入了两个偶应力mxz和myz(如图1所示)：弹性应力-应变关系可以表示为σ=Deε其中本构矩阵为这里,和G是传统意义上的拉梅常数，Gc则为Cosserat剪切模量.准静态条件下，若忽略体积力及体力偶的影响，则平衡方程可表示为LTσ=01.2 修正的Drucker-Prager准则适用于各向同性材料的Drucker-Prager准则以应力不变量的形式表示为F=q+Aφp+Bφ=0其中,p和q可分别用应力第一不变量及偏应力第二不变量表示如下:其中，P矩阵可表示为式(10)中参数Aφ、Bφ为内摩擦角φ和黏聚力c的函数，若Drucker-Prager准则为Mohr-Coulomb准则的内接圆，则可将其分别表示为通常对黏聚力考虑线性硬化或软化，即其中，c0为初始黏聚力，为硬化(软化)参数，p为等效塑性应变.基于应力不变量与组构张量，Pietruszczak与Mroz给出了一个能够描述各向异性材料的屈服准则[8]：η=η0(1+Ωijlilj)其中：η为描述屈服函数的材料参数，与应力状态有关；η0为η的平均值；Ωij 为描述η偏离η0程度的偏量.对正交各向异性材料Ωij有两个不等的特征值，对横观各向同性材料，其可用一个标量描述，对各向同性材料，Ωij为零，此时η=η0.li、lj为相对于材料主轴的加载方向.式(16)也可用Ωij的主值表示为对横观各向同性材料按如下方式建立坐标系，以竖直方向为y轴，以水平面为x-z 平面，建立整体坐标系Oxyz.令局部坐标系与材料主方向共轴，即各向同性面为1-3面，其法向为e(2)方向，其中方向e(3)与z轴重合，如图2所示.此时，Ω1=Ω3，且Ω1+Ω2+Ω3=0以及=1，式(17)可进一步表示为加载向量按如下方式计算：式中,表示应力张量，N=e(2)⊗e(2).为了简化，仅考虑内摩擦角的各向异性，即令则有式(10)、式(3)、式(1)构成了描述横观各向同性材料的屈服函数.对于岩土材料，通常采用非关联流动法则，塑性势可取为G=q+Aψp+Bψ类似地：其中，ψ为膨胀角.由一致性算法最终得到更新后的应力[17- 18]，σ=CαsE+(pE-KAψΔ)m其中：Δ其中FE=qE+AφpE+BE.以及一致性弹塑性切线模量矩阵：其中:可以验证，当采用关联流动法则时，该矩阵具有主对称性.以Abaqus软件的UEL子程序接口为平台，对文中发展的单元类型和弹塑性本构模型进行了数值实现[20].以平板压缩模型为例，单元类型为平面应变8节点减缩单元，模型尺寸为0.6 m×0.8 m，上下边界约束x方向自由度，左右边界自由，上下边界施加对称均布位移荷载uy，如图3(a)所示.计算中所采用的材料参数如表1所示.1)E为弹性模量；υ为泊松比；β为材料主方向角.Pietruszczak等[8]基于他们提出的各向异性屈服准则对横观各向同性材料物质点的理论分析显示，对于不同的加载方向，材料轴向强度随Ω1的演化规律不同，存在一个转折点，该点材料的轴向强度将独立于各向异性程度，且Zhong等[12]证明了此转折点所对应的材料主方向与材料固有属性有关.这里，为与已有理论结果进行对比，选取了结构的两个代表性单元，标号为Ⅰ和Ⅱ，分别位于板上边缘中间位置和剪切带附近，如图4所示，给出了两个单元1号积分点(见图3(b))上等效应力随加载位移的变化曲线，其中参数取β=30°，Ω1=-0.1.易知峰值qmax代表了该点的强度，图5(a)和5(b)则给出了积分点(物质点)上归一化的材料强度随各向异性程度及材料主方向的变化趋势，可见其趋势与已有理论结果[8]较为吻合(见图5(c),其中法向强度ζ=Δσ1/(2σ0))，从而验证了本程序开发的正确性.下面主要就材料主方向β及各向异性程度Ω1对结构极限承载力及破坏模式的影响做详细分析. 2.1 材料主方向β的影响由图2可知，β为各向同性面相对x轴逆时针转动的角度，描述了材料主轴相对整体坐标系的偏移，β的改变也可理解为加载方向的改变.图6为Ω1=-0.10，加载位移为0.05 m时等效塑性应变分布云图.值得注意的是，Ω1取负值意味着各向同性面法向(e(2)方向)有较大的摩擦强度[14].可以看出剪切带呈现了X型分布，当β=0°和β=90°时，由于加载方向和材料主方向保持一致或垂直，剪切带呈对称分布，如图6(a)和6(e)所示.当β≠0°和β≠90°时，剪切带两个分支呈现不对称性，沿各向同性面方向的剪切带较宽，等效塑性应变值较高，文中称之为强剪切带，与之共轭的另一条剪切带称为弱剪切带，如图6(b)和6(h)所示，与姚仰平等[4]的理论分析相比，其破坏面的位置可能就是强剪切带出现的位置，数值模拟显示还存在一个弱剪切带，只是在破坏时强剪切带占据了主导地位.由图6(b)-6(d)可见，当0°<β<90°时剪切带呈“/”型分布，且随着β的增大“/”方向剪切带逐渐变宽.反之，当90°<β<180°时，剪切带呈“\”型分布，且随着β的增大“\”方向剪切带逐渐变窄.对于一对互补的β角，如β=30°和β=150°，等效塑性应变峰值相等，但剪切带分布模式相反.随着β从0°演化到180°，强剪切带上等效塑性应变峰值出现了增大-降低-增大-降低的对称过程，弱剪切带则具有互补的变化过程，如图7所示.图8给出了随β变化材料的承载力-位移曲线,可以清楚地看到当各向同性面位于水平位置，即β=0°时，结构承载力最大，随着β的增大承载力逐渐降低，90°时最低，这与横观各向同性材料屈服准则中摩擦强度参数各向异性的分布规律一致. 2.2 各向异性参数Ω1的影响如前所述，Ωij为描述η偏离η0程度的偏张量，对于横观各向同性材料可以用一个标量Ω1来描述.需指出文中各向异性参数取值来自于文献[8]，依次取Ω1=0，-0.05，-0.1,-0.15,β=30°,加载位移为0.05 m，其他参数同表1.图9给出了等效塑性应变随各向异性程度变化的结果.可以看出，各向异性程度越大，X型剪切带的非对称性越强，强剪切带等效塑性应变的峰值越大，弱剪切带峰值则越小.图10为承载力-位移曲线，随着Ω1绝对值的增大结构承载力逐渐增大.这里同时给出了归一化的承载力随各向异性程度及材料主方向的综合变化规律，如图11所示，其中，纵轴表示各承载力峰值除以Ω1=0时的承载力峰值，可见结构整体响应亦呈现出了与图5类似的规律，即在β约为47°的位置出现转折点，β小于该值时，结构承载力随Ω1绝对值的增大而增大，大于该值时则相反.2.3 网格依赖性调查基于经典连续体理论，采用有限单元法计算弹塑性问题，当应变软化和局部化现象发生时，不可避免地会出现网格依赖性问题.本节通过对经典连续体和Cosserat连续体两种理论框架下的计算结果进行比较，调查了有限元计算的网格依赖性问题. 取材料主方向β=90°，Ω1=-0.1，其他材料参数同表1，网格密度分别取6×8，12×16，18×24，24×32，30×40，60×80进行比较.图12给出了基于经典连续体等效塑性应变分布图，可以看到随着网格加密，等效塑性应变的峰值逐渐增大，且剪切带宽度明显变窄.图13显示结构极限承载力随网格密度增大而逐渐降低，且模拟软化段的能力亦逐渐降低.图14则为基于Cosserat连续体得到的等效塑性应变分布图，随着网格加密，等效塑性应变峰值有所增大但剪切带宽度基本保持不变.图15显示虽然结构承载力有些许降低，但其变化越来越小，计算结果最终收敛于较为接近的值.需说明的是，图中网格为6×8时计算结果与其他结果相比差异较大，这是由于网格稀疏到一定程度，有限元计算本身的误差所致.经过比较可见引入Cosserat连续体后，网格依赖性问题得到了有效的解决.基于Pietruszczak 和 Mroz提出的各向异性屈服准则，考虑摩擦角为材料主方向及应力张量与微结构张量联合不变量的函数，提出了一个适用于横观各向同性岩土材料的修正的Drucker-Prager屈服准则，并基于Cosserat连续体推导了其一致性映射返回算法和切线模量矩阵.数值算例重点分析了材料主方向及各向异性程度对结构极限承载力及破坏模式的影响，结果表明：(1)材料主方向β=0°或β=90°时剪切带呈对称分布，前者结构承载力最大，后者最低.0°<β<90°时，随β增大，剪切带逐渐变宽，结构承载力逐渐降低，90°<β<180°则相反.β从0°到180°，强剪切带等效塑性应变峰值经历了增大-降低-增大-降低的对称过程，弱剪切带则相反.β互补时剪切带呈现相反的分布模式. (2)Ω1绝对值越大，材料各向同性面的强度越低；对于β≠0°、β≠90°的非对称情况，剪切带呈现的X型分布的非对称性更强.承载力变化趋势存在一个转折点，约为β=47°，β小于该值时，结构承载力随Ω1绝对值的增大而增大，β大于该值时则相反.(3)与经典连续体的模拟结果比较表明，文中基于Cosserat连续体的数值方法较好地解决了网格依赖性问题.† 通信作者: 徐远杰(1956-)，男，博士，教授，主要从事工程结构破坏数值仿真研究.E-mail:***************【相关文献】[1] ARTHUR J R F,Menzies B.Inherent anisotropy in a sand[J].Geotechnique,1972,22(1):115- 128.[2] ARTHUR J R F,CHUA K S,DUNSTAN T.Induced anisotropy in a sand[J].Geotechnique,1979,27(1):13- 30.[3] ODA M.Inherent and induced anisotropy in plasticity theory of granular soils[J].Mechanics of Materials,1993,16(1/2):35- 45.[4] 姚仰平,祝恩阳.横观各向同性土的简明破坏机制解释 [J].岩土力学,2014,35(2):328- 333. YAO Yang-ping,ZHU En-yang.Concise interpretation of damage mechanism for cross-anisotropic soil [J].Journal of Geotechnical Mechanics,2014,35(2):328- 333.[5] 贾乃文.岩土工程中横观各向同性Drucker屈服准则 [J].华南理工大学学报(自然科学版)，1993,21(2):49- 54. JIA Nai-wen.Drucker yield criterion of horizontally isotropic rock and soil [J].Journal of South China University of Technology(Natural ScienceEdition),1993,21(2):49- 54.[6] DUVEAU G,SHAO J F,HENRY J P.Assessment of some failure criteria for strongly anisotropic materials [J].Mechanics of Cohesive-Frictional Materials,1998,3(1):1- 26.[7] GAO Z W,ZHAO J D,YAO Y P.A generalized anisotropic failure criterion for geomaterials [J].International Journal of Solids and Structures,2010,47(22/23):3166- 3185.[8] PIETRUSZCZAK S,MROZ Z.Formulation of anisotropic failure criteria incorporating a microstructure tensor [J].Computer and Geotechnics,2000,26(2):105- 112.[9] PIETRUSZCZAK S,MROZ Z.On failure criteria for anisotropic cohesive-frictional materials [J].International Journal for Numerical and Analytical Methods in Geomechanics,2001,25(5):509- 524.[10] 余天堂,卢应发,SHAO J F,等.沉积岩的一种各向异性模型 [J].岩土力学,2002,23(1):47- 50. YU Tian-tang,LU Ying-fa,SHAO J F,et al.An anisotropic model of sedimentary rocks [J].Journal of Geotechnical Mechanics,2002,23(1):47- 50.[11] 余天堂.岩土材料固有各向异性的模拟 [J].岩石力学与工程学报,2004,23(10):1604- 1607. YU Tian-tang.Modeling of inherent anisotropy for geotechnical material [J].Chinese Journal of Rock Mechanics and Engineering,2004,23(10):1604- 1607.[12] ZHONG S Y,XU W Y,LING D S.Influence of the parameters in the Pietruszczak-Mroz anisotropic failure criterion [J].International Journal of Rock Mechanics and Mining Sciences,2011,48(6):1034- 1037.[13] LADE P V.Modeling failure in cross-anisotropic frictional materials [J].International Journal of Solids and Structures,2007,44(16):5146- 5162.[14] LADE P V.Failure criterion for cross-anisotropic soils [J].Journal of Geotechnical and Geoenvironmental Engineering,2008,134(1):117- 124.[15] LADE P V.Three-dimensional failure in geomaterials:experimentation and modeling [M]∥Constitutive Modeling of Geomaterials.[S.l.]:Springer Berl in Heidelberg,2013:47- 58.[16] CHU X H,CHANG J F,JIANG Q H,et al.Numerical simulation of progressive failure in transversely isotropic geomaterials[C]∥Advance in Heterogeneous Material ncaster:Destech Publication Inc.,2011:965- 968.[17] DE BORST R.Simulation of strain localization:a reappraisal of the Cosserat continuum [J].Engineering Computations,1991,8(4):317- 332.[18] LI X K,TANG H X.A consistent return mapping algorithm for pressure-dependent elastoplastic Cosserat continua and modeling of strain localization [J].Computers & Structures,2005,83(1):1- 10.[19] 余村,楚锡华,唐洪祥,等.基于Cosserat连续体的颗粒破碎影响研究 [J].岩土力学,2013,34(增刊1):67- 79. YU Cun,CHU Xi-hua,TANG Hong-xiang,et al.Study of effect of particle breakage based on Cosserat continuum [J].Journal of Geotechnical Mechanics,2013,34(Sup 1):67- 79.[20] ARSLAN H,STURE S.Finite element analysis of localization and micro-macro structure relation in granular materials(Part I):Formulation [J].Acta Mechanica,2008,197(3/4):135- 152.。
基于Cosserat理论的平面等参元
第19卷第1期2010年3月计算机辅助工程Co m puter A ided Eng i n eeri n gV o.l19N o.1M ar.2010工程数值仿真与CAE算法Num erica l Si m u lati o n of Eng ineeri n g and C AE A l g orithm文章编号:1006-0871(2010)01 0012 05基于Cosserat理论的平面等参元齐 磊, 张若京(同济大学航空航天与力学学院,上海 200092)摘 要:为解释材料在微尺度下的尺度效应,基于C osserat理论,从势能泛函驻值条件出发提出构造8节点Serendipity平面等参元,并建立平面有限元法.每个节点拥有3个独立节点自由度,分别为2个方向的线位移和1个逆时针方向的角位移.用该方法分析含中心小孔的无限平板在单轴拉伸情况下的应力集中问题.数值计算结果与Cosserat理论的解析解非常符合,表明应力集中因数k 受泊松比 ,常数c及a/l值的影响很大;由于偶应力的存在,小孔周围的应力分布明显小于经典弹性力学理论的预测.通过对材料常数c的调节可以将该方法推广应用于基于M i n d li n偶应力理论的数值分析中.关键词:Cosserat理论;偶应力;平面有限元法;等参元;尺度效应;应力集中中图分类号:O241.82 文献标志码:APlane isopara metric ele m ent based on Cosserat theoryQ I Le,i Z HANG Ruo ji n g(Co lleg e of A erospace Eng.&M echanics,TongjiU n i v.,Shanghai200092,China)Abst ract:To expla i n the sca le effect o f m aterials under the m icroscale,based on Cossera t theory and deduced fro m the stationary conditi o n o f potential energy functiona,l an eight node Serendipity p lane isopara m etric e le m ent is proposed and a fi n ite ele m ent m ethod for p lane is developed.There are threei n dependent node freedo m s at each node,i n clud i n g li n ear displace m ent in t w o d irecti o ns andcounterc l o ckw ise angular d isp lace m en.t The stress concentration prob le m o f a i n fi n ite flat p late w ith a centra l s m all hole is ana l y zed by t h e m ethod i n the case of un iax ial tensi o n.The num erical calculation resu lts are i n good agree m en t w ith the analytica l so lution based on Cosserat theory,wh ich i n d icate tha t the stress concen trati o n factor k is strong l y infl u enced by Po i s son rati o ,constant c and the value o f a/l;the stress d istribution ar ound the little ho le is sign ificantly less than the prediction of classical elastic ity theory due to the ex istence o f couple stress.The m ethod can be applied to the num erica l ana l y sis based on M indlin couple stress theory by t h e regu lation o f constant c.K ey w ords:Cosserat theory;couple stress;finite ele m ent m ethod for p lane;isopara m etetric ele m en;tsca le effec;t stress concen trati o n收稿日期:2009 07 20 修回日期:2009 10 10基金项目:国家自然科学基金(10772136);教育部博士点基金(20060247016);上海市重点学科建设项目(B302)作者简介:齐 磊(1983 ),男,河北廊坊人,硕士研究生,研究方向为有限元数值分析和计算,(E m ail)fayak@163.co m;张若京(1946 ),男,北京人,教授,博导,博士,研究方向为固体力学,(E ma il)zhangr@j t ongj.i 0 引 言经典的连续介质力学理论已被广泛应用于许多工程领域,如航空航天、车船、机械和土木等.事实证明,这种理论适合应用于具有宏观结构力学性能的分析.在微尺度下,由于尺度效应的影响,材料的部分力学性能不同于经典连续介质力学理论的预测.近年来的大量实验研究也表明,当金属试件的尺寸在微米或亚微米量级时,材料具有很强的尺度效应.如在细铜丝扭转实验中,FLECK等[1]发现,当铜丝直径从170 m下降到12 m时,无量纲的扭转硬化增加至3倍;在N i材质的薄梁弯曲试验中, STOLKEN等[2]观察到当梁的厚度从100 m下降到12.5 m时,无量纲的弯曲硬化显著增加,而在单轴拉伸情况下这种尺度效应不存在.在微米量级尺度下,微观硬度试验与颗粒增强金属基复合材料中也观察到尺度效应.当压痕深度从10 m减至1 m时,金属硬度增加1倍[3 6];对于以碳化硅颗粒加强的A l Si基复合材料,LLOYD[7]等观察到当保持颗粒体积比为15%的条件下,将颗粒直径从16 m减为7.5 m时,复合材料的强度显著增加.以上这些实验现象用经典弹性力学理论都无法解释,因为在其本构模型中不存在任何尺度参数,无法预测材料的尺度效应.因此,很多学者致力于建立考虑包含尺度影响的连续介质模型,并取得丰硕成果,这其中,COSSERAT两兄弟[8]于1909年提出Cosserat理论.TOUPI N[9],M I N DLI N[10]和KO I TER[11]等对Cosserat理论加以简化又提出各自的偶应力理论.1 Cosserat理论和M i nd li n偶应力理论简介在Cosserat理论中(见图1),由于材料点处有偶应力m 的存在,使得Cauchy应力t 为非对称,图1 偶应力微元体因此将不对称的Cauchy应力t 分解为对称的! 和反对称∀ 的两部分:!=(t +t )/2∀ =(t -t )/2(1)以p 表示体力,q表示体力偶,则平衡方程可表示为! , +∀ , +p =0(2)m , +∀12-∀21+q=0(3)以s 表示应变,则几何变形方程为s11=u1,1, s12=u2,1-#, s22=u2,2s21=u1,2+#, ∃1=#,1, ∃2=#,2(4)式中:u=[u,v]T为线位移向量;#为微元体的微转动;∃1和∃2分别为与m1和m2相对应的形变,称为曲率.在Cosserat理论中,#为与线位移u和v相互独立的自由度,因此应变张量与应力张量一样都是非对称的.将应变张量分解为对称的% 和反对称的& 2部分:% =(s +s )/2& =(s -s )/2(5)式中:% 和& 是分别与对称应力! 和反对称应力∀ 相对应的形变.本构方程为! =D s &∋%&∋∀ =D a &∋&&∋, m =E ∃(6)对于平面应力问题D s &∋=2G∋ &∋∋+1-∋ ∋&∋(7)D a &∋=2cG∋ &∋∋,E =4G l2∋ (8)式中:G为剪切模量; 为泊松比;∋为K ronecker符号;l为表征材料尺度效应的特征长度;c为无量纲的常数.对于平面应变问题,则可令泊松比为 / (1- ).边界条件可以表示为u =u- , #=#-(9)T- =t n=(! +∀ )nM-=m n(10)式中:u- 和#-分别为边界线位移和角位移;T- 和M-分别为作用于边界上的面力和面力偶;n为边界上的单位外法线矢量.势能泛函为((u ,#)=12!A(! % +∀ & +m ∃ )d A-!A(u p +#q)d A-13第1期齐 磊,等:基于Cosserat理论的平面等参元!)(u T - +#M -)d )(11)在M i n dli n 的偶应力理论[12]中,微转动被约束,从而导致应变张量对称.又因为曲率∃包含线位移的2阶导数,这就要求位移形函数必须满足C 1连续,这无疑给单元构造带来很大困难.另外,在M i n dli n 的偶应力理论中,应力张量是非对称的,应变张量是对称的,因此本构方程无法确定所有应力分量与应变分量的关系.本文采用Cosserat 理论,避免C 1连续的要求,通过对材料常数c 的调节可以把Cosserat 理论的数值方法推广应用于M i n dli n 的偶应力理论中.[13]2 有限元法将上述表达式写成列阵形式:!={!x ,!y ,!xy }T∀={∀xy }, m ={m x ,m y }T(12)其中!x =!11, !y =!22, !xy =!12∀xy =∀12, m x =m 1, m y =m 2%={%x ,%y ,%xy }T&={&xy }, ∃={∃x ,∃y }T(13)其中%x =%11, %y =%22, %xy =2%12&xy =2∃12, ∃x =∃1, ∃y =∃2对于平面应力问题D s=G 21- 21-2 1- 21- 0001D a=c G [1], E =4Gl200l2(14)对于平面应变问题,可令泊松比为 /(1- ),则本构关系为!=D s%∀=D a &, m =E ∃(15)势能泛函为(=12!A(%TD s%+&TD a&+∃TE ∃)d A -!A(u Tp+#q )d A -!)(u TT -+#M -)d )(16)如图2所示,采用8节点Serend i p ity 单元,节点自由度分别为2个方向的线位移和1个逆时针方向的角位移U ={u,v,#}T(17)图2 8节点Serend i p ity 单元应力应变列向量分别为!*={!x ,!y ,!xy ,∀xy ,m x ,m y }T%*={%x ,%y ,%xy ,&xy ,∃x ,∃y }T(18)对位移进行差值U e=N e r e=[N 1I 3,N 2I 3,N 3I 3,N 4I 3,N 5I 3,N 6I 3,N 7I 3,N 8I 3]re(19)式中:I 3为3阶单位矩阵;N i 为Serendipity 单元双二次函数;r e为单元位移向量r e={u 1,v 1,#1,u 2,v 2,#2,∀,u 8,v 8,#8}T(20)将几何方程式(4)代入式(19)得单元应变列向量%*=LU e=x0 y -y 00y x x 000-2xyTu v #(21)将式(19)~(21)代入式(16),并令∋(=0得K e r e=Pe(22)其中K e=!A eB TDB d A (23)B =LNe(24)D =D s3#3DaE 2#26#6(25)P e=P eA +P e)=!A eN T {p x,p y,q }Td A +!)eN T{T -x ,T -y ,M -}Td )(26)3 数值算例考虑无限大弹性平板的平面应变问题.圆孔半径为a,远处作用x 方向的均匀分布力为p 0.由于问题的对称性,取1/4平板进行分析,见图3.取计算域为30a,平板厚度取单位长度1,共划分单元54414计 算 机 辅 助 工 程 2010年个,节点1743个(含边中节点).由于关注的是小孔周围的应力集中情况,故在小孔周围局部进行网格细化,以得到更好的结果,见图4.图3 有限元计算网格 图4 小孔周围局部网格细化计算结果见图5和表1,可知由于有偶应力影响,小孔周围的应力集中因数k与经典弹性力学理论解相比有所降低,且k的大小受泊松比 影响较大. 越大,k就越大.本文解与文献[14]中基于Cosserat理论的解析解非常符合,误差不足千分之一.需要说明的是,对于平面应变问题,将不得不面对不可压问题( 接近极限值0.5时).本文所采取的措施是将式(16)中体积分的第1项应变能拆分成斜偏应变能和体积应变能两部分,并分别进行精确积分和2#2减缩积分,以保证在不可压状态下体积应变能所产生的单刚奇异性及整个单刚的非奇异.表1括号内数据均为采用精确积分的结果.图5 应力集中因数k随 变化曲线表1 应力集中因数k随 变化(常数c=0.5,材料特征长度l=a)00.10.20.30.40.490.4990.49990.49999解析解[14] 2.3055 2.35252.40302.4573 2.51612.57322.57912.57982.5798本文解 2.3074 2.35432.40472.4591 2.51782.57482.5819(2.5845)2.5867(2.6819)2.5947(4.6269)误差/%0.080.080.070.070.070.060.110.270.58为便于观察,将应力结果数据转化为极坐标形式.其中,经典解是指基于经典弹性力学理论的解析解!∗/p0沿∗=90∃的r方向分布,见表2和图6.表2 !∗/p沿∗=90∃的r方向分布(常数c=0.5,材料特征长度l=a,泊松比 =0.3) r/a1 1.03 1.051.081.12 1.15 1.181.221.251.331.40经典解3.000 2.803 2.6312.4812.3492.232 2.129 2.0181.9221.7651.642解析解[12]2.457 2.341 2.2382.1462.0641.990 1.924 1.8521.7881.6801.592本文解2.459 2.239 2.1522.3492.0661.996 1.926 1.8561.7871.6871.587误差/%0.080.340.040.280.100.300.100.22-0.050.42-0.31注:经典解是指基于经典弹性理论的解析解.图6 !∗/p沿∗=90∃的r方向分布从表2及图6可知,在∗=90∃时,!∗/p0的值沿r方向逐渐减小,当半径在1.4a处时,!∗/p0的偶应力解(包括Cosserat理论的解析解[14]和本文的数值解)已经与经典解十分接近,偶应力影响已经非常小.因此,在进行网格细化时,只对r<1.4a的单元进行细化.图6中还有2条数据是对不可压问题的计算结果,可以看出平均误差在0.5%以内,表明此单元完全能够用于不可压材料的计算.表2中第2,4,6,8和10这5列数据的本文解都是边中节点数据,且9,10和11这3个节点为单元网格加密过程中过渡单元中的节点,因此误差略微偏大.从图7可以发现,k随着c的增大而减小.当c在0到10变化时,其对k的影响较大.c越接近0,小孔周围的应力分布情况越接近于经典弹性理论的解析解;c值越大,结果越接近于M indlin偶应力理论的解析解,这与文献[13]中的结论一致.当c>15第1期齐 磊,等:基于Cosserat理论的平面等参元100时,其结果可认为是M i n d li n 偶应力理论的数值解.从能量泛函中也不难找到其理论依据:在M i n dli n 偶应力理论中宏观转动与微转动相等,即#=(u 2,1-u 1,2)/2.将此式化为u 2,1-u 1,2-2#=0并将罚函数的方法引入到位能泛函中,修改后的位能泛函与式(11)完全一致,式中体积分第2项即为引入的约束条件.其中,& =u , -u , -2#e ,∀ =G c & ,c 即可认为是罚参数.因此,当c %&时,Cosserat 理论的解即为M indli n 偶应力理论的解,而一般情况下c 并不会取得很大,因此在式(21)中的&xy 项虽然出现形函数N i 及其1阶导数耦合的情况,但不会产生类似剪切锁死现象.然而,为了使此单元同时也能应用于M i n d li n 偶应力理论的数值计算,c 的取值很大,这将会遇到剪切锁死.如前面所述,可将此积分项认为是引入的罚函数积分项,在进行计算时可将其单独提出来进行2#2减缩积分,以保证所计算出来的单刚为奇异,与讨论不可压时对体积应变能的积分处理方法相同,而本文的实际计算过程是将斜偏应变能项与弯曲应变能项合并、体积应变能项与反剪切应变能项合并,分别以精确积分和2#2减缩积分计算单刚,最后再相加形成完整单刚.图7 !∗/p 0沿∗=90∃的r 方向分布从表3及图8发现,应力集中因数k 受a /l 的影响较大.当a /l 值小于6时,k 的值随a /l 值的增大而迅速增大;当a /l 值大于6时,k 的值受a /l 值的影响渐渐减小;当a /l 值大于25时,可近似认为应力k 已经不受偶应力的影响.表3 应力集中因数k 随a /l 变化(常数c =0.5,泊松比 =0.3)a /l 0.112346810152025解析解[14]2.366 2.457 2.5792.6792.7522.844 2.895 2.9252.9622.9772.985本文解2.366 2.459 2.5822.6812.7552.847 2.898 2.9282.9652.9802.988误差/%0.020.080.120.070.110.110.100.100.100.100.10图8 k 随a /l 变化曲线4 结 论基于Cosserat 理论,构造出性能较好的8节点Serendipity 单元,数值计算精度较高,能够很好地反映出偶应力在作用域中的应力分布.通过对带中心小孔受单轴拉伸的无限平板进行分析,发现应力集中因数k 受泊松比 ,常数c 及a /l 值的影响规律,得到小孔周围应力分布情况,能与基于Cosserat 理论的解析解很好符合,同时也可以应用于M i n d li n 偶应力理论的求解.参考文献:[1] FLECK N A ,MULLER G M,AS H BY M F .S trai n grad ient p l asti city :theory and exp eri m en t[J].Act a M etall u rgicaM ateriali a ,1994,42(2):475 487.[2] STOLKEN J S ,EVANS A G .A m i crobend test m et hod for m eas uri ng t he plas tici ty lengt h scale[J].A cta M at eri ali a ,1998,46(14):51095115.[3] NI X W D .M ec h an i ca l prop erties of t h i n fil m s[J].M etall u rgical&M ateri als T ran s A,1989,20(11):2217 2245.[4] STEL M AS H ENKO N A,W ALLS M G,BRO W N L M ,et a l .M icro i ndentations on W and M o oriented s i ngle crystals :an ST M st udy[J].ActaM et all urg i ca M ateri alia ,1993,41(10):2855 2865.[5] M A Q i ng ,CLARKE D R .S i ze depend ent h ardness of silver s i ngle crystals[J ].J M ater Res ,1995,10(4):853 863.[6] POOLE W J ,ASHBY M F ,FLECK N A .M i cro hardnes s of anneal ed and w ork h arden ed copper pol ycrystals [J ].S cri p t a M etall u rgicaM at eri ali a ,1996,34(4):559 564.(下转第26页)参考文献:[1] BEL IS J,DEP AUW J,CALLEWAERT D,et a l.Fail u re m echanis m s and residua l capacit y of annealed gl ass/SGP la m i nat ed bea m s at roo mte mp erature[J].E ng Fail ure An a,l2009,16(6):1866 1875.[2] 张其林.玻璃幕墙结构设计[M].上海:同济大学出版社,2007:57 58.[3] 马眷荣,臧曙光,丁丽梅.夹层玻璃力学模型的探讨[J].航空材料学报,1998,18(3):57 60.[4] J G J102 2003 玻璃幕墙工程技术规范[S].北京:中国建筑工业出版社,2003:22 33.[5] 贺瑞丽,郑百林,武秀根,等.用APDL对悬索桥结构进行有限元分析[J].计算机辅助工程,2008,17(2):30 33.[6] VALLABHAN C V G,DAS Y C,RA M ASAMUDRA M an j un atha.Properti es of PVB i n terlayer u sed i n la m i n ated gl ass[J].J M at C i vil E ng,1992,4(1):71 76.(编辑 廖粤新)(上接第16页)[7] LLOYD D J.Particl e rei n f orced al um i nu m and m agnesi um matri x co m posites[J].IntM at er Rev i e w s,1994,39(1):1 23.[8] COSSERAT E,COSSERAT F.Theorie de corp def or m ab l es[M].E diti ons A H er m ann,Pari s:1909.[9] TOUPI N R A.E l astic m ateri als w it h couple stresses[J].A rch Rati onalM ech&An a,l1962,11(1):385 414.[10] M I NDLI N R D.M i cro s tr u cture i n li near el asti city[J].A rch Rati onalM ech&An a,l1964,16(1):51 78.[11] KOI TER W T.C oup l e stresses i n the t heory of elasti city,I and II[C]//P roc K N ed Ak ad W et B,1964:17 44[12] M I NDLI N R D.In fl u ence of coup le stress on stress con cen trati ons[J].E xp M ech,1962,3(1):1 7.[13] 李雷,吴长春.基于Coss erat理论的应变梯度非协调数值研究[J].工程力学,2004,21(5):166 171.[14] 刘俊,黄铭,葛修润.考虑偶应力影响的应力集中问题求解[J].上海交通大学学报,2001,35(10):1481 1485.(编辑 廖粤新)(上接第21页)参考文献:[1] 王瑞,陈海霞,王广峰.ANSYS有限元网格划分浅析[J].天津工业大学学报,2002,21(4):8 11.[2] 江雄心,万平荣.整体规划技术在三维有限元网格生成中的应用[J].锻压技术,2006:119 122.[3] 吴淑芳,付沛福,李运兴,等.平面区域四边形网格自动生成技术[J].计算结构力学及其应用,1996,13(2):208 216.[4] 李华,李笑牛,程耿东,等.1种全四边形网格生成方法 改进模板法[J].计算力学学报,2002,19(1):16 191.[5] 严蔚敏,吴伟民.数据结构(C语言版)[M].北京:清华大学出版社,2002:216 222.[6] 李华,程耿东,顾元宪.1种新的全四边形网格快速生成方法 模板法[J].计算结构力学及其应用,1996,13(1):25 39.[7] 金晶,吴新跃.有限元网格划分相关问题分析研究[J].计算机辅助工程,2005,14(2):75 78.(编辑 廖粤新)。
带静电斥力弹性杆的Cosserat方程及其数值模拟
Ｖｏ １ ． ２ ６ ＮＯ ． ４
ＮＯ Ｖ ．２ ０ １ ３
文章编号 ： １ ０ ０ ６ —１ ０ ３ ７ （ ２ ０ １ ３ ） ０ ４—０ ０ ０ ４ —０ ６ ｄ ｏ ｉ： １ ０ ． ３ ９ ６ ９ ／ ｊ ． ｉ ｓ ｓ ｎ ． １ ０ ０ ６ —１ ０ ３ ７ ． ２ Ｏ １ ３ ． １ １ ． ０ ２
动 力学 。 ］ 分析 对 Ｄ ＮＡ运 动 的动力 学机理 进行 了探索 。 在 弹性杆 动力 学模 型 的研 究 中 ， 接触 和 自接 触 问题是 近年来 研究 的重要 内容 和难 点 。 目前研 究 的 问
题 如发 型 、 网和 ＤＮＡ 的缠绕 等 问题 都 涉 及到 接 触 和 自接 触 。由于 接 触＿ ６ 的过 程 改变 弹 性杆 系统 的 拓扑 结 构， 给数 值模 拟 带来 很 大 困难 。为 了解决 这一 问题 ， 许多 学者提 出一些 处 理方法 ］ ， 但都 不 能解决一 般 的 问题 。生物学 家研 究 Ｄ ＮＡ 大分 子的避 免 接触 和穿 越 的方 法 给我 们 重要 的 启示 。ＤＮＡ 大 分子 本 身 带有 同
的 弹性杆 的 Ｃ ｏ ｓ ｓ ｅ ｒ ａ ｔ 方程 。最后利 用数值 方法 进行数 值模 拟 。 关键 词 ： 弹性杆 ； 接触 ； 静 电斥力 ； Ｃ ｏ ｓ ｓ ｅ ｒ ａ ｔ 方程； 欧拉参数 中图分 类号 ： ０１ ７ ４ ． ２
主题 分 类号 ： ４ ２ Ｂ ２ ０
文 献标识 码 ： Ａ
真 研究 中 ， 发型、 网和结 等复杂 的 弹性杆模 型 的拓扑 结构 和力 学 性质 成 为弹 性杆 研 究 的新 的热点 。另 外 ， 弹
性杆 模 型在 Ｄ ＮＡ 等 大分子 的研究 中得 到应 用 。通 过复 杂的生 物实 验技术 ， 人们 已经 确定 了 Ｄ ＮＡ分 子 可 以 模 型化 为弹性 细杆 ， 并测 出了相 应 的力学 参数 如杨 氏模 量 、 剪 切模 量 和 泊松 比等… ， 并 通 过 对 弹性 杆模 型 的
预应力曲杆的Cosserat动力学模型
预应力曲杆的Cosserat动力学模型曹登庆;宋敉淘【摘要】To dynamically model slender structures conveniently by Cosserat rod theory, the dynamic equations for pre-stressed curved Cosserat rods are established within the framework of this rod theory by defining variables to describe their initial configurations, their pre-stresses and their initial deformations, which move along with their cross sections. Based on their corresponding specific assumptions, the dynamic equations for a cable with initial sag and a circular arch beam, which have wide applications in bridge structures, are explored within the framework of pre-stressed curved Cosserat rods, respectively. The results are in agreement with their corresponding equations obtained by other methods in references. Since the Cosserat theory is geometrically exact, the dynamic equations for pre-stressed curved Cosserat rods derived here are not only retaining all geometric nonlinear characteristics but also of universality, which can be effectively and efficiently exploited to nonlinear dynamically model for slender structures in practical engineering situations.%为了更方便地使用Cosserat杆模型对各种细长结构进行动力学建模，在该杆模型的框架下通过引入描述固连在杆横截面随其一起运动的初始变形、初始横截面转动以及横截面间的初始接触力和力矩等变量，建立了考虑预应力的细长曲杆的Cosserat动力学模型．在得到的模型的框架下，基于相应的物理假设，从数学演绎的角度推导了桥梁建设中广泛应用的具有初始垂度的拉索以及圆拱梁的动力学方程，结果与文献中采用其他方法得到的相应的动力学方程一致．由于Cosserat理论采用了精确的几何构形，这里导出的细长预应力曲杆动力学模型不仅保留了所有的几何非线性特征，而且具有很好的普适性，通过它可以容易地导出实际工程问题中细长结构的非线性动力学模型．【期刊名称】《哈尔滨工业大学学报》【年(卷),期】2012(044)011【总页数】6页(P1-6)【关键词】Cosserat杆模型;预应力曲杆;具有初始垂度的拉索;圆拱梁【作者】曹登庆;宋敉淘【作者单位】哈尔滨工业大学航天学院,哈尔滨150001;哈尔滨工业大学航天学院,哈尔滨150001【正文语种】中文【中图分类】O326;O34118世纪的Bernoulli和Euler认为梁之所以能抗弯是梁在运动过程中轴向伸缩的结果，且梁横截面受到的弯矩正比于梁轴线的曲率.Euler-Bernoulli梁模型用梁轴线的挠度描述梁的运动，梁横截面的转动用挠度对空间坐标的导数来描述.Coulomb于1776年首次将力和力矩的平衡方程用于对梁微元进行力学分析，建立了较Euler-Bernoulli梁更为精确的动力学方程，并开始考虑梁承受扭矩的情况［1］.Kirchhoff于1859年计算了杆单元的应变，从能量的角度得到了杆的动力学方程，并提出了基于 Kirchhoff动力学比拟［2－3］的 Kirchhoff弹性杆模型.该模型忽略杆轴线的伸缩变形和横截面的剪切变形，并将杆横截面视为刚体，以描述杆轴线的3个线位移参数以及描述各横截面的3个姿态参数为未知变量建立杆的动力学方程.Rayleigh于1877年引入梁横截面惯性矩的影响，考虑轴向和横向振动对梁进行建模，并指出在高频振动时该修正是很必要的［1］.Timoshen-ko ［4－5］引入梁横截面间剪切力的影响对梁进行建模，并指出 Euler-Bernoulli梁和 Rayleigh梁为Timoshenko梁的两种特殊情况.更加精确的杆的动力学模型的问世随着工程设计和分析的需要而显得越来越必要和迫切.近年来，基于Cosserat兄弟和其他人的早期工作，形成了关于细长杆的Cosserat 介质理论.Antman［6］在他的著作中运用该理论的基本思想建立了柔性杆和薄壳的Cosserat模型.Tucker和Wang［7］和 Cao 等［8－9］进一步完善了Cosserat杆理论，并分别应用于对钻柱的非线性动力学研究和MEMS中复杂结构三维非线性动力学模拟.Burton和Tucker［10］回顾了Cosserat杆模型的发展以及它近年来在各个领域中的应用.国内学者在Cosserat杆模型上也做了很多重要的工作.刘延柱等［11－12］从更偏重于力学分析的角度得到了任意初始构形圆截面Cosserat杆的动力学方程，在其框架下推导了Rayleigh梁、Kirchhoff杆、以及Timoshenko梁的动力学方程，并对基于精确的Cosserat杆模型螺旋杆进行了稳定性分析.Cosserat杆模型在Kirchhoff杆模型的基础上考虑了杆轴向的伸缩变形以及横截面间剪切力引起的横截面的转动，由于所作的假设更少，Cosserat杆模型较以往的杆模型都更为精确.鉴于Antman所建立的 Cosserat杆模型［6］几何形式的直观性，本文的工作将在其框架下进行.本质上Cosserat杆模型适用于对具有任何初始构形和预应力的杆进行建模，但为了讨论的方便，以往文献中给出Cosserat杆动力学方程的同时，在Cosserat杆的本构关系中实际上都只将杆的由预应力引起的初始变形用1个抽象的符号表示，或者直接只考虑无预应力直杆(Non-Prestressed Straight Rod，NPSR)的情况［7－9］，而且在 Cosserat杆的动力学方程中，他们并没有对初始量和动力学意义上的量进行区分.如何在Cosserat杆动力学方程以及其本构关系中加入杆的初始构形和预应力项以描述具有预应力的曲杆(Prestressed Curved Rod，PCR)的动力学行为是关键.本文的目的便是依据文献中NPSR的Cosserat模型的基本思想，用Cosserat杆模型对PCR进行建模，并称得到的模型为PCR的Cosserat动力学模型.在该模型的框架下，再根据具体的假定，便可以很明确地从数学演绎的角度推导各种特殊情况下的一维细长结构的动力学模型.本文中，矢量用小写斜体黑体字母表示，如u，v和w;张量用大写斜体黑体字母表示，如K，I和J;重复下标为哑指标，表示采用爱因斯坦求和约定;运算符⊗为两张量间的张量积运算.考虑一根细长杆，在初始时刻t0时各横截面质心组成的曲线(质心线)上定义拉格朗日弧坐标s，s∈［0，L0］，L0为杆的初始长度.如图1所示，假设运动过程中t时刻杆的质心线为r(s，t)，分别以r(s，t)上的各点为原点定义局部的右手笛卡尔直角坐标系s ∈ ［0，L0］|→{d1(s，t)，d2(s，t)，d3(s，t)}，并将坐标系固定在相应位置处杆的横截面上.将d1(s，t)垂直于相应位置处杆的横截面，将d2(s，t)和d3(s，t)置于相应位置处杆横截面内.忽略杆运动和变形过程中密度和横截面面积的变化，于是Cosserat杆模型的动力学方程可以由牛顿第二定律得到:式中:ρ(s)为杆的密度;A(s)为杆的横截面面积;n(s，t)和m(s，t)分别为杆横截面间接触力和接触力矩;f(s，t)和l(s，t)分别为外力密度和外力矩密度(可以包括重力，阻尼力，电磁力，气动力等);I(s，t)为截面的惯性矩张量.由于Cosserat杆模型考虑了剪切力引起的杆横截面的转动，于是r'(s，t)不再始终垂直于相应位置处的横截面，令将t0时刻的di(s，t)(i=1，2，3)记为即(s)=di(s，t0).若杆为NPSR，即在t0时刻横截面与横截面间接触力和接触力矩为零的直杆，则(s)实际上并不依赖于变量s，记为.同时，沿着s的方向将保持不变，即和(s)也独立于变量s，分别记为.于是，可以通过固连在横截面上的局部坐标系定义矢量u(s，t)和w(s，t)以分别度量杆横截面的转动随s和t的变化，式(3)和(4)分别意味着(参见文献［6－9］)于是，Cosserat杆的本构关系可以表示为如下形式:式(7)中的张量 K(s，t)和 J(s，t)以及式(1b)中的张量I(s，t)的分量形式分别为式中:E(s)为杆的弹性模量;G(s)为杆的剪切模量;dA=dξdη，η和ξ分别为坐标系基矢量d2(s，t)和d3(s，t)对应坐标轴的坐标;Ψ(ξ，η)为横截面的翘曲函数，由圣维南原理知其对于所有横截面都相同［13］;κ(s)为 Timoshenko剪切系数，现有不同的理论对其进行物理意义上的解释，Shames和Dym［14］系统探讨过这个问题，但不论哪种解释，κ(s)的值都取决于横截面的形状.张量I(s，t)和J(s，t)可以通过选择 d2(s，t)和 d3(s，t)的方向达到对角化，以使得它们的形式变得更为简洁.下面阐述式(7)中NPSR的Cosserat模型本构关系背后的力学依据.式(7a)中出现的d1(s，t)表示t0时刻杆的质心线没有剪切变形，只有沿着杆相应处横截面法向的变形.由于采用了杆初始构形上的弧坐标，于是初始时刻质心线对弧坐标的微分总是为d.杆初始构形上任意一微元ds便可以用ds d作为其长度和方向的度量，它随着杆的整体运动而运动，作为微元的初始状态始终垂直于相应位置处的横截面，于是变为ds d1(s，t).同时，变形后的微元为v(s，t)ds.于是微元的正应变(沿d1(s，t)的方向)以及剪切应变(分别沿着d2(s，t)和 d3(s，t)的方向)分别为在 d1(s，t)，d2(s，t)和d3(s，t)上的投影.至此，便不难得到式(7a).式(7b)表示横截面间的接触力矩由横截面间相对转角来度量，对其分析可类似于对式(7a)所作的分析那样，通过定义随着杆整体运动的杆横截面的初始转角而实现.由于在前面定义的初始时刻局部坐标系间没有相对转角，于是初始时刻横截面间无相对转角.由上面分析可见，本构关系(7)的获得，关键在于定义随着杆横截面一起运动的各种初始量，这给PCR的本构关系的建立及分析预应力对动力学方程的影响提供了思路.假如细长杆为PCR，即杆的横截面与横截面间接触力n(s，t)和接触力矩m(s，t)在t=t0时不为零(即 n0(s)=n(s，t0)≠ 0，m0(s)=m(s，t0)≠0)，则在静平衡位置处杆在预应力和初始外力f0=f(s，t0)和初始外力矩l0=l(s，t0)作用下保持平衡，由牛顿第二定律得其中r0(s)=r(s，t0)为静平衡位置时杆的质心线.杆质心线上各质点的线位移分为静位移r0(s)和动位移rd(s，t)两部分，即在运动过程中，杆的各微元一方面有着自身的变形和运动，一方面又随着杆的整体运动而运动.杆横截面间的初始接触力、初始接触力矩、初始相对转角以及杆的初始变形固连在定义在相应位置处横截面上的局部坐标系上，它们随着杆的整体运动而运动.它们相对于局部坐标系不随时间变化，然而相对于全局坐标系却是不断变化的，也即预应力矢量和初始变形矢量是随时间变化的变量，于是可以将杆质心线的初始变形以及杆横截面间的初始相对转角、初始接触力、初始接触力矩分别记为于是式(11)和(12)中式(13)中v(s)和u(s)分别为在初始局部坐标系下的分量.ξij(s，t)=d(s)·dj(s，t).由第1节的内容知道u(s)由d(s)的定义决定.为了保证类似于式(5)的定义有意义，应将d(s)取为s个连续的单位矢量函数.曲杆无法做到像直杆那样将d(s)独立于s.不考虑杆运动过程中产生裂缝的情况，d(s)取为垂直于相应位置处的横截面能保证d(s)连续.同时，取d(s)和d(s)为任意的连续单位矢量函数都不会影响到最后的结果.因为d(s)的选取影响着u(s)的同时，也影响着 di(s，t)，进而导致 u(s，t)－ (s，t)的值不依赖于d(s)的选取.将外力(外力矩)分为初始力(初始力矩)和动态力(动态力矩)两部分，即将式(10)－(16)代入式(1)中，并将得到的两方程分别减去式(9)中的两个方程得到PCR的Cosserat模型的动力学方程为当n0、m0和r0为零时，式(17)退化为式(1)的形式.值得注意的是，拉索可以视为杆忽略弯扭刚度以及横截面间剪切刚度的特殊情况，于是具有初始垂度的拉索的动力学方程可以由方程(17a)来描述，而且此时张量K需要退化成如下形式:另外需要注意的是，前面的推导是在假定曲杆的横截面、初始质心线、横截面间初始接触力以及初始接触力矩为已知的前提下进行的.实际在一般情况下，对于曲杆，特别是具有预应力的曲杆，如何选择杆横截面并不是一件容易的事情.同时，求解杆的预应力也非易事.不过，这并非我们所关心的主要问题.而且，对于很多特殊情况，确定杆的上述初始量并不困难，作为例子我们将在第3节中采用本节建立的PCR的Cosserat模型推导具有初始垂度的空间拉索以及圆拱梁的动力学方程的过程中可以看到.如图2所示，在全局右手笛卡尔直角坐标系{o;e1，e2，e3}下，具有初始垂度的拉索运动过程中的形状为其中 W1(s，t)，W2(s，t)和 W3(s，t)分别为拉索各质点在e1，e2和e3方向的线位移.由于不考虑拉索横截面间剪切力的影响，于是d1(s，t)沿着拉索的切线方向，即将式(18)代入式(2)得全局坐标系下张量K写成矩阵形式为将式(19)－(21)代入式(12a)，并根据式(13)得另一方面拉索的动应变e(s，t)为于是由式(22)和(23)得若拉索初始张力为H，则将式(28)进行泰勒级数展开，忽略高于一阶的小量并整理得将阻尼从外力中分离出来，即其中c为阻尼系数.将式(18)，(29)和(30)代入式(17a)中得到拉索动力学方程为其中分别为拉索在方向e1，e2和e3所受到的时变分布力.假如拉索在平面内运动，取式(31)的平面形式，可以发现其和文献［15］中用广义Hamilton原理得到的结果一致.如图3所示，初始时刻圆拱梁圆心为o.如果不考虑圆拱梁的预应力，则可以通过如下方式选取圆拱梁的横截面:通过o点，取垂直于圆拱梁所在平面的平面，并让它与圆拱梁相截，则截出的圆拱梁的截面即可选为圆拱梁的横截面.在圆拱梁初始构形上对质心线每个质点s建立右手笛卡尔直角坐标系{o;eθ(s)，er(s)，ez(s)}，其中eθ(s)沿着质心线在s处的切线方向并且指向s增加的方向，er(s)从o指向s，ez(s)=eθ(s)×er(s).可见，ez(s)实际上不随s而变化.圆拱梁初始时刻质心线为圆拱梁在运动过程中质心线的形状为其中Wθ(s，t)和 Wr(s，t)分别为 t时刻质点 s在eθ(s)和er(s)方向的线位移.取记圆拱梁横截面相较于初始构形转过的角度为ψ(s，t)，并以逆时针转动为正，于是将式(32)和(33)分别代入式(2)，将式(34)代入式(6)，将式(35)分别代入式(5)和(6)，并根据式(13)－(15)得由于不考虑圆拱梁的预应力，即二维情况下，张量 K在坐标系{o;eθ(s)，er(s)，ez(s)}中的分量形式为张量J和 I在坐标系{o;eθ(s)，er(s)，ez(s)}中的分量形式分别为同样地将阻尼从外力中分离出来，并将式(33)，(36)－(39)代入式(17)中得到圆拱梁的动力学方程为其中分别为圆拱梁在eθ和 er方向上所受到的分布力，lz为圆拱梁在ez方向上所受到的分布弯矩，并且若圆拱梁为小变形，则可以用式(40)的线性化形式来描述圆拱梁的运动，此时式(41)退化为将式(42)和(8)代入式(40)，并忽略非线性项得若忽略阻尼和外力，式(43)和文献［16］中用广义Hamilton原理得到的圆拱梁的线性动力学方程一致.1)建立了初始局部坐标系下杆各初始量(质心线的初始变形，横截面的初始转角，横截面间的初始接触力和接触力矩)与运动过程中它们反映在PCR的Cosserat动力学方程和本构关系式中的相应的初始变量的具体关系.2)将初始变量耦合进Cosserat杆动力学模型及其本构关系中建立了PCR的Cosserat动力学模型及其本构关系.3)作为特例，桥梁中广泛使用的具有初始垂度的空间拉索以及圆拱梁的动力学方程在PCR的Cosserat动力学模型的框架下作了推导.其过程体现了Cosserat杆模型包容性强以及数学推导的透明性和一贯性，简化了繁杂的物理和力学层面上的分析过程.PCR的Cosserat动力学模型的建立使得各种特殊的一维细长结构的动力学建模变成了一个很纯粹的数学推导过程，极大地方便了使用Cosserat杆模型对各种特殊的一维细长结构的建模.【相关文献】［1］ARAFAT H N.Nonlinear response of cantilever beams［D］.Blacksburg，Virginia:Virginia Polytechnic Institute and State University，1999:17－18.［2］LOVE A E H.A treatise on the mathematical theory of elasticity［M］.New York:Dover，1944:399 －400.［3］刘延柱.弹性细杆的非线性力学［M］.北京:清华大学出版社，2006:14－18.［4］TIMOSHENKO S P.On the correction for shear of the differential equation for transverse vibrations of prismatic bars［J］.Philosophical Magazine，1921，41(1):744 －747.［5］TIMOSHENKO S P.On the transverse vibrations of bars of uniform cross section ［J］.Philosophical Magazine，1922，43(1):125－131.［6］ANTMAN S S.Nonlinear problems in elasticity［M］.New York:Springer-Verlag，1995:13－359.［7］TUCKER R W，WANG C.Torsional vibration control and Cosserat dynamics of a drill-rig assembly［J］.Meccanica，2003，38(1):143－159.［8］CAO Dengqing，LIU Dongsheng，WANG C H－T.Nonlinear dynamic modelling for MEMS components via the Cosserat rod element approach［J］.Journal of Micromechanics and Microengineering，2005，15(6):1334 －1343.［9］CAO Dengqing，TUCKER R W.Nonlinear dynamics of elastic rods using the Cosserat theory:modelling and simulation［J］.International Journal of Solids and Structures，2008，45(2):460 －477.［10］BURTON D A，TUCKER R W.Practical applications of simple Cosserat methods［C］//Mechanics of Generalized Continua.In:Advances in Mechanics andMathematics.Berlin:Springer，2010，21:87 －98.［11］LIU Yan-zhu.On dynamics of elastic rod based on exact Cosserat model［J］.Chinese Physics B，2009，18(1):1 －8.［12］刘延柱，薛纭.基于精确Cosserat模型的螺旋杆稳定性分析［J］.应用数学和力学，2011，32(5):570－578.［13］TIMOSHENKO S P，GOODIER J N.Theory of elasticity［M］.New York:McGraw-Hill，1951:275 －313.［14］SHAMES I H，DYM C L.Energy and finite element methods in structural mechanics ［M］.Chaman Offset，Delhi，India:New Age International Publishers，1991:197 －203. ［15］SOFI A，MUSCOLINO G.Dynamic analysis of suspended cables carrying moving oscillators［J］.International Journal of Solids and Structures，2007，44(21):6725 －6743. ［16］JAVID F，ESMAILZADEH E，YOUNESIAN D.An investigation into the vehicle-curved bridge dynamic interaction［J］.International Journal of Automotive Engineering，2011，1(3):235－243.。
地基渐进破坏及极限承载力的Cosserat连续体有限元分析
第28卷第11期 岩 土 力 学 V ol.28 No.11 2007年11月 Rock and Soil Mechanics Nov. 2007收稿日期：2005-11-09 修改稿收到日期：2006-04-30基金项目：国家自然科学基金资助项目（No. 50278012, No. 19832010）；国家973项目（No. 2002CB412709）资助。
作者简介：唐洪祥，男，1973年生，博士，主要从事岩土力学与工程的理论与应用研究工作。
E-mail ：tanghx@文章编号：1000－7598－(2007) 11－2259－06地基渐进破坏及极限承载力的 Cosserat 连续体有限元分析唐洪祥1，李锡夔2（1.大连理工大学 海岸与近海工程国家重点试验室，大连 116023；2.大连理工大学 工程力学系，大连 116023）摘 要：利用Cosserat 连续体理论和所发展的有限元数值方法，模拟了地基由应变软化引起的以应变局部化为特征的渐进破坏过程，并从等价塑性应变的发展变化，阐述了渐进破坏过程对所能发挥的极限承载能力的影响。
结果表明，Cosserat 连续体模型能有效地模拟由应变软化引起以应变局部化为特征的渐进破坏现象，对地基等土工结构物有必要进行渐进破坏分析。
同时指出，在求解软化型土体地基的极限承载力时，如果仍按传统的极限平衡或极限分析理论进行分析，可能得出偏于危险的结果。
关 键 词：Cosserat 连续体；有限元法；应变软化；应变局部化；渐进破坏分析；极限承载力 中图分类号：TB 115 文献标识码：AFinite element analysis of Cosserat continuum for progressive failureand limit bearing capacity of soil foundationTANG Hong-xiang 1, LI Xi-kui 2(1.State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116023, China;2. Department of Engineering Mechanics, Dalian University of Technology, Dalian 116023, China)Abstract: Based on the derived finite element formulations of Cosserat continuum, progressive failure phenomena of the soil foundation, characterized by strain localization due to strain softening, are numerically analyzed. The influence of progressive failure process on the limit bearing capacity is illustrated by the developments of equivalent plastic strain. Numerical results demonstrate the effectiveness of the Cosserat continuum finite elements in preserving the well-posedness of the localization problem and simulating the progressive failure phenomena characterized by strain localization due to strain softening, necessity of the progressive failure analysis for earth structures such as the soil foundation etc.. Simultaneously it is pointed out that the unsafe results are to be obtained as the limit bearing capacity of strain softening soil foundation is analyzed by traditional theories of limit equilibrium and limit analysis.Key words: Cosserat continuum; finite element method; strain softening; strain localization; progressive failure analysis; limit bearing capacity1 引 言超固结黏土和密实砂的三轴剪切等实验室试验结果，以及对现场土工结构物如边坡、路堤与地基的滑动破坏所观察到的现象表明，土体会出现剪胀以及在峰值之后会出现应变软化行为。
Cosserat弹塑性模型在ABAQUS中的数值实施
Cosserat弹塑性模型在ABAQUS中的数值实施彭从文;董龙彬;吴群【摘要】The subroutine UEL embedded in commercial FEM software ABAQUS is utilized to develop a user element based on the material of pressure-dependent elasto-plastic Cosserat continuum model. The user element adopt plane eight-node isoparametric unit which include three degree of freedom (two translation and one rotation) at each node. The Drucker-Prager material model with associate flow rule is used, and the stress integration algorithm is four-order Runge-Kutta method. The user element is used to analyze the influence of mesh density and the material characteristic length in the material localization problem. The results show that the thickness of shear band and the equivalent plastic strain are independent of mesh density, with the increase of the material characteristic length, the thickness of shear band and the equivalent plastic strain decrease constantly.%利用大型有限元软件ABAQUS提供的接口程序UEL,开发了压力相关弹塑性Cosserat连续体材料的用户单元.采用平面八节点等参单元,包括平动与转动三个自由度,考虑相关联流动的Drycjer-Prager材料模型,应力计算采用四阶龙格-库塔法.利用该单元分析了材料局部化问题中网格密度与材料特征长度的影响.结果表明,网格密度对材料剪切带厚度与等效塑性应影响很小；随着特征长度增大,剪切带厚度增大,等效塑性应变峰值减小.【期刊名称】《武汉工程大学学报》【年(卷),期】2011(033)006【总页数】5页(P102-106)【关键词】Cosserat模型；FEM；ABAQUS；UEL；局部化【作者】彭从文;董龙彬;吴群【作者单位】长江大学城市建设学院,湖北荆州434023;荆州市城市规划设计研究院,湖北荆州434000;深圳中广核工程设计有限公司,广东深圳518031【正文语种】中文【中图分类】TB120 引言偶应力理论是微极理论的一个特例，Cosserat兄弟最先提出了完整的偶应力理论[1]，Toupin[2]，Mindlin等[3]对该理论作了进一步的发展和完善.Cosserat理论引入了旋转自由度和相应的微曲率，引入了与微曲率能量共轭的偶应力、以及具有“特征长度”意义的尺度参数.该理论可以较好地处理网格敏感性和控制方程失去椭圆性的问题，近年来，由于细观力学、非均质力学的发展，Cosserat理论重新受到关注，逐渐成为研究热点之一.数值方法是重要的研究手段之一，为了提高计算精度与效率，基于大型通用数值计算平台的二次开发方法得到了广泛的应用.目前,许多数值计算平台没有内嵌Cosserat计算模型，关于Cosserat模型的二次开发还不多见[4-6]. ABAQUS是目前最流行、功能最强的商用有限元软件之一，该软件可以进行结构静、动力分析，具有强大的非线性计算能力、丰富的材料库及良好的扩充功能，自1997年进入我国以来，越来越多的国内企业和研究机构采用ABAQUS作为产品研发和科学研究的工具.本文采用ABAQUS的用户接口程序，研究压力相关弹塑性Cosserat连续体模型的用户子程序UEL实施方法.1 Cosserat连续体模型考虑Cosserat连续体平面问题，每个材料点有三个自由度.u=应力、应变分别定义为：S=E=几何方程为：ui,j=uj,i-eijkωk(1)κij=ωj,i(2)静力平衡方程为：σij,j+fi=0(3)mij,i+eijkσik+qj=0(4)式(1)～(4)中，fi、qj分别为体积力与体积力偶；eijk为排列算子；ux,uy,ωz分别是平面内平移与转动自由度；mxz,myz偶应力；κxz,κyz为微曲率.对于弹性材料，其本构关系为[7]S=DeE(5)式(5)中，G,v,a,l分别是材料的剪切模量、泊松比、COSSERAT材料参数及特征长度.对于弹塑性Cosserat材料，采用基于Drucker-Prager屈服准则的弹塑性Cosserate连续体模型，其屈服函数与流动势函数分别为[8](6)(7)式(7)中,其中，α,β,A,B为材料参数，σ,s分别为应力与偏应力张量.Cosserat连续体弹塑性本构关系推导方法同经典连续介质力学，应力应变增量关系为(8)式(8)中，内变量采用等效塑性应变,为硬化项定义为2 UEL实现方法及计算流程2.1 子程序编写注意事项(1)ABAQUS提拱了两类单元自定义方法.一类是线性单元，可以通过结果文件或INP文件直接给出，不需要编写UEL；另一类就是通用单元，通过UEL子程序定义；(2)UEL子程序中更新变量与分析问题类别有关.同一个模型中可能遇到不同的分析步，如地应力平衡、静力分析、摄动步分析等，因此，编写UEL时要区别处理；(3)UEL允许自定义荷载.包括集中荷载、均布荷载及弯矩等.其中，对于均布荷载，须定义荷载标志号；(4)自定义单元在ABAQUS/CAE中不可见.若想在ABAQUS/CAE中显示自定义单元变形图，可以将ABAQUS标准单元与自定义单元绑定，同时将标准单元材料参数设为小值.UEL子程序中所有输出变量均通过SDV写入结果文件(.fil、.dat)，其分量在ABAQUS/CAE中不可见.2.2 AMATRX与RHS计算UEL界面与ABAQUS内核主要通过AMATRX与RHS等变量进行数值传递.本文采用八节点等参单元，设单元节点位移、插值函数与位移应变转换矩阵分别为d、N与B，单元内任一点位移u及应变ε为u=∑Nidi=Nd(9)ε=Bd(10)式(9～10)中,B=将式(9)、(10)代入Cosserat介质虚功方程(11)进行方程离散.(11)式(11)中，f=，t=对于线性问题，结合材料本构关系，得到式(12).(12)式(12)中，将式(12)简记为Kd=f′(13)对于非线性问题，采用Newton-Raphson方法，将式(13)改写为ψ(d)=K(d)d-f′(14)设ψ(d)为具有一阶导数的连续导数，初始近似值为d(0)，第i次迭代的近似值为d(i).将函数ψ(d)在d(i)处展开，保留线性项，忽略高阶项得:(15)d(i+1)=d(i)+Δd(i)(16)式(15)中，与Ψ(d)分别为UEL中需更新变量AMATRX与RHS.2.3 计算流程计算流程如图1所示.图1 UEL流程图Fig.1 Flow of UEL3 性状分析3.1 有限元模型模型几何尺寸10×5 m2，采用三种不同网格密度，单元数分别是15×30、20×40、25×50.边界条件为底端竖向固定，左侧水平向固定.材料弹性模量25 GPa，泊松比0.3，内摩擦角35°，粘聚力1.5 MPa，软化模量15 MPa.采用相关联流动法则.顶部采用位移加载方式，加载量为20 mm，加载方向向下.为了触发局剪切带，对左下角单元弱化处理.采用高斯完全积分，四阶龙格-库塔显式应力积分方法.3.2 计算结果计算模型分析了不同网格密度及特征长度的影响，计算结果如图2～6所示.(1)局部化带的客观性.图2为经典连续介质理论得到的等效塑性应变云图，图3与图4分别是采用Cosserat理论计算得到的等效塑性应变云图与应力应变曲线.由图可知，采用Cosserat理论计算时，随着网格密度增加，剪切带厚度与等效塑性应变峰值基本不变.当采用经典连续介质理论计算时，计算结果有明显的网格依赖性，随着网格密度增加，软化带逐渐变窄，等效塑性应变峰值也不断增大,计算收敛趋于弱化.图2 不同网格密度等效塑性应变云图Fig.2 Contour of equivalent plastic strain with different mesh density注：(a)15×30, (b)20×40,(c)25×50图3 不同网格密度等效塑性应变云图(特征长度0.15)Fig.3 Contour of equivalent plastic strain with different mesh density (characteristic length:0.15)注：(a)15×30, (b)20×40,(c)25×50图4 不同网格密度下应力位移曲线Fig.4 Stress-displacement curve with different mesh density(2)特征长度的影响.Cosserat理论引入特征长度作为正则化机制，特征长度决定Cosserat连续体模型模拟应变局部化问题的能力并影响局部化剪切带宽度大小.图5与图6分别为不同特征长度下采用Cosserat理论计算得到的等效塑性应变云图和应力位移曲线.由图2～6可知，随着特征长度增大，剪切带厚度增大，等效塑性应变峰值减小，材料软化模量降低.图5 不同特征长度下等效塑性应变云图(网格20×40)Fig.5 contour of equivalent plastic strain with different characteristic length (mesh density: 20×40)注：(a)0.2 (b)0.15 (c)0.05图6 不同特征长度的影响Fig.6 Effect of different characteristic length4 结语基于ABAQUS接口程序UEL，开发了压力相关弹塑性Cosserat连续体材料的用户单元，并采用该单元分析了有限元网格密度及材料特征长度对材料局部化的影响.结果表明，采用Cosserat理论计算时网格密度对材料剪切带厚度、等效塑性应变影响很小，这也在一定程度上说明本文方法的正确性.要特别说明的是，基于ABAQUS平台进行二次开发能有效地利用现有程序代码，减小开发工作量，缩短有限元程序开发周期，极大地提高科研工作效率.参考文献：[1] Cosserat E, Cosserat F. Theorie des Corps Deformables [M].Paris: Herman et Files, 1909.[2] Toupin R A. Elastic materials with couple stresses [J].Archive Rational Mechanics and analysis, 1962(11): 385-414.[3] Mindlin R D, Tiersten H F. Effects of couple stresses in linear elasticity [J]. Archive Rational Mechanics and analysis, 1962(11): 415-448.[4] 杨乐, 吴德伦, 许年春. 偶应力理论的层状岩体洞室数值模拟[J]. 重庆建筑大学学报, 2008, 30(3):73-77.[5] 尹雪英, 杨春和, 李银平. 层状盐岩体三维Cosserat介质扩展本构模型的程序实现[J]. 岩土力学, 2007,28(7):1415-1420,1426.[6] 朱珍德, 秦天昊, 王士宏, 等. 基于Cosserat理论的柱状节理岩体各向异性本构模型研究[J]. 岩石力学与工程学报,2010,29(增2): 4068-4076.[7] Sharbati E, Naghdabadi R. computational aspects of the cosserat finite element analysis of localization phenomenon[J]. Computational materials science,2006(38):303-315.[8] HAYDAR ARSLAN. Localization analysis of granular materials in cosserat elastoplasticity -formulation and finite element Implementation [D]. Colorado: University of Colorado, 2006, 85-88.。
Cosserat连续体模型中本构参数对应变局部化模拟结果影响的数值分析
值分析 , 令 h 0 = - 5 0
0 1h0 , 0 2h0 , 0 4h0 , h0 , 2 0h0 , 10 0h 0 , Gc = 0 5 G, 其他参数与 3 节中所采用的相同。 图 6( a) 和图 6( b) 分别表示所得相应的荷载 位移曲线和剪切带宽 度。 可见 , 当 hp = 0 1 h0 时, 荷载 位移曲线在达到峰 值后没有软化段, 类似于理想弹塑性材料的行为, 没 有明显的剪切带显示与应变局部化现象; 随着软化 模量绝对值的增大, 荷载 位移曲线在达到峰值后 出现明显的软化段且软化段越来越陡 , 并有明显的 剪切带显示且剪切带的宽度越来越小; 当取 hp = 10 0h0 时, 荷载 位移曲线在达到峰值后迅速下降, 也就没能得到现剪切带与应变局部化现象。 再对 3 节中图 3 所示的二维平板压缩问题进 行数值 分析。 令 h0 = - 1 5 105 N/ m 2 , 分别 取 hp = 0 1 h0 , 0 5 h0 , h 0 , 2 0h 0 , 5 0 h0 , 10 0 h0 , 其他
唐洪祥* 1 ,
摘
李锡夔2
( 1. 大连理工大学 海岸与近海工程国家重点实验室 , 大连 116024; 2. 大连理工大学 工程力学系 , 大连 116024) 要 : 基于所发展的压力相关弹塑性 Co sser at 连续体模型及相应的数值方法 , 以一维剪切层 及二维平板 压缩问
题为例 , 数值分析了 Co sserat 连续体模型中的本构参数 Cosserat 剪模 、 软化模量及内部长度参 数对应变局 部化数 值模拟结果的影响 。 结果表明在一定取值范围内 , Co sser at 剪模对 数值模拟结 果几乎没有 影响 , 并给出了具 体数 值计算时的取值范围 ; 软化模量绝对 值越大 , 后破 坏段的 荷载 位移 曲线越 陡 , 计 算得到 的剪切 带宽度 越窄 ; 内部 长度参数越大 , 后破坏段的荷载 位移曲 线越平缓 , 计算得到的剪切带越宽 。 关键词 : Cosserat 连续体模型 ; Cosserat 剪模 ; 软化模量 ; 内部长度参数 ; 数值分析 中图分类号 : O 241 文献标识码 : A
基于Cosserat弹性细杆模型的电缆形态仿真
火炮发射过程中身管温度场及弯曲度的有限元计算
火炮发射过程中身管温度场及弯曲度的有限元计算
赵金辉;何忠波;傅建平;王江涛
【期刊名称】《火力与指挥控制》
【年(卷),期】2011(036)005
【摘要】利用ANSYS有限元软件,建立了某火炮身管的实体模型.对火炮发射过程中膛内的热交换进行模拟仿真,得出火炮发射后身管复杂的温度场分布.并综合考虑身管自重以及热作用的影响,对其进行热结构耦合分析,求出身管发射后不同时刻的弯曲量,为连发射击时的射击修正提供重要参考.
【总页数】4页(P106-109)
【作者】赵金辉;何忠波;傅建平;王江涛
【作者单位】军械工程学院,石家庄,050003;军械工程学院,石家庄,050003;军械工程学院,石家庄,050003;解放军75610部队,广东,东莞,518110
【正文语种】中文
【中图分类】TJ304
【相关文献】
1.混凝土箱梁桥施工过程中的水化热温度场实例及有限元计算 [J], 唐云清;桂玉枝
2.复合材料身管压力场/温度场全等效应力有限元计算 [J], 王惠源;薄玉成;佘建生
3.速射身管发射状态下的温度场及热应力的有限元分析 [J], 胡志刚;赵建波
4.LOCA条件下环境介质影响感应加热过程中锆合金内部温度场的有限元计算 [J], 席航;谢东升;张瑷月;康锐;吴璐;潘虎成
5.LOCA条件下环境介质影响感应加热过程中锆合金内部温度场的有限元计算 [J], 席航;谢东升;张瑷月;康锐;吴璐;潘虎成
因版权原因，仅展示原文概要，查看原文内容请购买。
基于Cosserat理论的微梁振动特性的尺度效应
摘要 不少微观实验已经证实, 微尺度领域 材料的力学性能存在尺度效应。采用偶应力理论( 又称 Cosserat 理论) 研 究微梁振动特 性( 主要是固有频率) 的尺度效应。文中首 先对偶 应力理 论进行简 介, 然后采用 Hamilton 变分原 理推导 基 于 Cosserat 理论的微 梁无阻尼自由振动的微分方程, 分析微梁固有频率对微 尺度的依赖性。结果表 明, 当 微梁的厚度 减 小到可以和材料的本征长度相比时, 微梁的固有 频率将显著增大。
002
机
械
强
度
2007 年
在微电子机械系统中( micro electro mechanical sys tem, MEMS) 中, 大量出现各种微梁结构, 如微加速度计 惯性元件的支撑梁、振动 式微机械陀螺的折 叠梁等。 这些微梁的厚度大多在微 米量级或亚微米 量级。显 然, 沿用经典的弹塑性理论预测其力学性能有一定的 局限性。目前, 国内外对微梁力学性能的尺度效应的 研究已有不少文献报道, 但这些研究主要集中在静态 方面, 对微梁动态方面的尺度效应的研究则还未见报 道。
康 新等: 基于 Cosser at 理论的微梁振动特性的尺度效应
003
弹性理论可设细长微梁的位移为
SIZE EFFECT ON THE DYNAMIC CHARACTERISTIC OF A MICRO BEAM BASED ON COSSERAT THEORY
康 新 1 席占稳2 ( 1. 南京理工大学 力学与工程科学系, 南京 210094)
( 2. 南京理工大学 机械工程系, 南京 210094) KANG Xin1 XI ZhanWen2
1 引言
经典弹塑性理论是建立在材料的连续、均匀假设 之上的, 忽略了材料内部微结构对材料性能的影响, 这 对于非均匀变形场在宏观尺度范围内是足够精确的。 但在微尺度下( 一般指构件特征尺度在微米或亚微米 量级) , 材料内部微结构的非均匀性将使材料的力学性 能出现异常, 表现出明显的尺度效应, 即微材料的力学 性能依赖于微尺度的大小, 这一点已被不少微观实验 所证实[ 1~ 4] 。而经典弹塑性理论的本构关系中不包含 任何与材料微结构有关的本征长度参数, 因此经典弹 塑性理论不能预测微尺度下材料力学 性能的尺度效 应。
钢筋弯曲过程的有限元模拟
钢筋弯曲过程的有限元模拟
孟欣佳;张立香
【期刊名称】《河北工程大学学报：自然科学版》
【年(卷),期】2011(28)2
【摘要】本文针对钢筋在弯曲过程中的大变形及非线性的特点,运用有限元分析软件ANSYS/LS-DYNA建立了钢筋弯箍系统的有限元模型,对钢筋的弯曲过程进行了数值模拟,并分析了钢筋在不同时刻的弯曲应力变化情况。
当钢筋弯曲角度接近于零时,钢筋产生的最大应力发生在与弯箍转套接触部位;随着时间的增大,与箍筋模接触部位的应力迅速增大,且应力最大值逐渐由钢筋的边缘部位向中心部位靠近。
【总页数】4页(P30-32)
【关键词】钢筋弯曲;弯箍;数值模拟;弯曲应力;ANSYS/LS-DYNA
【作者】孟欣佳;张立香
【作者单位】山西机电职业技术学院;中国农业大学工学院
【正文语种】中文
【中图分类】TG333.26
【相关文献】
1.碳纤维增强复合材料桥梁弯曲破坏过程有限元模拟 [J], 张威;汪金花;景晓颖
2.双凸模差速挤压方形弯曲管件工艺过程有限元模拟 [J], 宋继顺;杜德恒;马叙;尹康;崔宏祥
3.铝合金三维弯曲管件差速挤压过程有限元模拟 [J], 杜德恒;宋继顺;马叙;尹康;崔
宏祥
4.管材三维无芯弯曲过程有限元模拟 [J], 李秋;王华君;田梦芸;孔祥志;彭春宇
5.基于流动速度场控制的弯曲管件双凸模挤压过程有限元模拟 [J], 宋继顺;杜德恒;崔宏祥;马叙;尹康
因版权原因，仅展示原文概要，查看原文内容请购买。
压力相关弹塑性Cosserat连续体模型与应变局部化有限元模拟
3
(1)
压力相关弹塑性 Cosserat 连续体模 型
为描写压力相关弹塑性本构行为，具体考虑如
式中：u x，u y 和 ω z 分别为平面内的平移和旋转自由 度。相应地，应变和应力向量分别定义为
下形式的 Drucker-Prager 屈服准则： (2) (3)
F = q + Aϕ σ h + B = 0
收稿日期：2003–12–03；修回日期：2004–02–24 基金项目：国家自然科学基金资助项目(50278012，19832010)；国家重点基础研究发展规划(973)项目(2002CB412709) 作者简介：李锡夔(1940–)，男，博士，1962 年毕业于大连工学院机械系，现任教授、博士生导师，主要从事计算力学、岩土与岩土环境力学方面的 教学与研究工作。E-mail：xikuili@。
σ = Deε e
其中，
De =
(7)
0 0 0 0 λ λ λ + 2G λ 0 0 0 0 λ + 2G λ λ 0 0 0 0 λ λ + 2G 0 0 0 G + GC G − GC 0 0 0 0 0 0 G − GC G + GC 0 0 0 0 0 2G 0 0 0 0 0 0 0 0 2G
PRESSURE-DEPENDENT ELASTO-PLASTIC COSSERAT CONTINUUM MODEL AND FINITE ELEMENT SIMULATION OF STRAIN LOCALIZATION
LI Xi-kui1，TANG Hong-xiang2
(1. State Key Laboratory of Structural Analyses for Industrial Equipment，Dalian University of Technology，Dalian 116024，China； 2. School of Civil and Hydraulic Engineering，Dalian University of Technology，Dalian 116024，China)
层状岩体地下洞室的Cosserat 理论有限元分析
由于单元体中两层弹模、泊松比不同，为保持 各层曲率相同，应力分布沿整个截面不再连续和等 梯度（图 3） 。
第3期
杨 乐等：层状岩体地下洞室的 Cosserat 理论有限元分析
983
展模型退化成仅含一种介质的情况，此时所得系数 与文献[5]推导系数近似一致。 2.3 Cosserat 扩展介质的屈服条件 对于互层岩体而言，由于考虑两种岩层的体积 含量不同而破坏形式不同。当某一种岩层体积含量 过小，则可忽略其厚度，仅考虑它的法向和切向刚 度，视为文献[2－7]中所研究的各层力学特性相同 的层状岩体。因此，这种情况是互层岩体的一种特 例。笔者在文中研究的是体积含量较均衡、物理特 性不同的互层岩体。若两种岩层的弹模相差较大， 相对软岩层会比硬岩层更早出现压剪破坏及弯曲破 坏。应用 Mohr-Coulomb 屈服准则，剪切破坏的屈 服条件：
，但均没有考虑岩层抵抗弯曲变形
的能力；引入偶应力的 Cosserat 理论却能体现这种 特性。 20 世纪 80 年代以来， Cosserat 连续介质理论 逐渐被引入岩土工程。国内外学者为验证该理论的 有效性、适用场合在数值分析的程序实现方面做了 － 大量的研究工作[3 6]。结果证明，该理论对于层状 岩体具有很好的适用性，能较好地模拟层状岩体顺
2.2
T
;
（1）
Cosserat 扩展介质的弹性本构方程[8] Cosserat 介质扩展模型是指单元体包含了物理
性质不同、厚度不同的材料，且满足假设条件：两 种材料的界面黏合性好，无滑动、张开或嵌入（图 2） 。考虑相邻两岩层间的细观位移协调，建立宏观 平均意义的应力： 1 h/2 c 11 h / 2 11dx2 A 11A B 11B h 1 h/ 2 c 12 h / 2 12 dx2 A 12 A B 12 B h 偶应力：
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
要 ：由于 传 统 连续 介 质 理论 的本 构 关 系 中 不 包 含 任 何 长 度 相 关 的量 纲 ， 此 它 难 以 预 测 微 电 子 工 业 设 因
计 与 制 造 中常 常 遇 到 的 微 米 尺 度 下 的尺 寸 效 应 问题 ． 过 建 立 Ｃ ｓｅａ 理 论 的变 分 原 理 和 有 限 元 基 本 方 程 ， 通 ｏｓｒｔ 进 而 编 制 了 相 应 的 有 限 元 程 序 ＧＦ Ｍ２ 模 拟 了微 米 尺 度 下 悬 臂 梁 的 尺 寸 效 应 ． 究 结 果 表 明 ， 着 内禀 长 Ｅ Ｄ， 研 随 度 的 增 大 ， 的抗 弯 刚 度 也 增 大 ． 究 结 果 证 明 了 Ｃ ｓｅａ 介 质 理 论 分 析 材 料 尺 寸 效 应 的可 行 性 ． 梁 研 ｏ ｓｒｔ
ａ ｃ ｌ ￡ ｅ ｔｓｍ ｕｌ ｔｏ ｆｍ ｉ ｒ — ｅ ｍ ｎｄ ｓ ａ ｅ ｅ ｔ ｃ ｉ ａｉｎ ｏ ｃ ｏｂ ａ
ＺＨ ＡＯ ｎｇ ，ＺＨ ＥＮＧ ｎ — ｅ ，ＣＨ ＥＮ ｕ ・ Ｂｉ Ｙｉ ｇ ｒ ｎ Ｃｈ ｎ。。
（ ． ｐ ｒ ｍｅ ｔｏ ｃ ｉ ｃ ｕ ａ ｇ ｎ ｅ ｉｇ，Ｌ ｇ ｓ ｉ Ｅｎ ｉ ｅ ｒｎ ｉ ｅ ｓｔ ｆＰＬＡ， １ Ｄｅ ａ ｔ ｎ ｆＡｒ ｈ ｔ ｔ ｒ ｌＥｎ ｉ ｅ ｒｎ ｅ ｏ ｉｔ ｇ ｎ ｅ ｉ ｇ Ｕｎ ｖ ｒ ｉ ｏ ｃ ｙ
ｃ ｏ ｅ ｅ ｔ ｏ ｃ ｎ ｓ ｒ ．Ｂａ ｅ ｎ ｔ ｅ ｇｅ ｒ ｌＣｏ ｐｌ— ｔ ｅ ｓｔ ｏ ｙ，ｔ ａ ｉ ｔｏｎ ｌｐｒｎｃｐ ｅ ｒ — ｌｃ ｒ ｎｉｓ ｉ ｄｕ ｔ ｙ ｓ ｄ ｏ ｈ ｎｅ ａ ｕ ｅ ｓ ｒ ｓ ｈｅ ｒ ｈｅ ｖ ｒａ ｉ ａ ｉ ｉ ｌ ａ ｄ ｆｎｉｅ ｅ ｅ ｅ ｔ ｂ ｓ ｃ ｅ ａ ｉ ｎ ｏ ｓ ｅ ａ ｈｅ ｒ ｓ ｅ ｔ ｂｌｓ ｄ，ｔ ｏ ｒ ｓ ｏｎ ｎ ｉ ｎ ｉ ｔ ｌ ｍ ｎ ａ ｉ ｑｕ ｔｏ ｆ Ｃｏ ｓ ｒ ｔ ｔ ｏ ｙ ｗａ ｓ ａ ｉｈｅ ｈｅ ｃ ｒ ｅ ｐ ｄｉ ｇ ｆ— ｎｉｅ ｅｅ ｅ ｏ ｒ ｍ ｔ ｌｍ ｎｔ ｐｒ ｇ ａ ＧＦＥＭ ２ ｗａ ｓｇ ｄ ｉ ｒ ｒ ｎ９ Ｄ ｓ ｄｅ ｉ ｎｅ ｎ Ｆｏ ｔａ ０，ｔ ｉ ｅ ｅ ｆ ｃ ｆ ａ ａ ｉｅ ｅ ｈｅ ｓｚ ｆ ｅ ｔ ｏ ｃ ｎｔｌ ｖ ｒ ｂｅ ｍ ｎ ｍ ｉｒ ｎ ｓ ａ ｅ ｗａ ｉ ｕ ａ ｅ ． ａ ｉ ｃ ｏ — ｃ ｌ ｓ ｓｍ ｌ ｔ ｄ Ｋｅ ｒ ｓ：Ｃｏ ｓ ｒ ｔ ｔ ｏ ｙ；ｖ ｒａ ｉ ｎａ ｉ ｉ ｌ ｙ ｗｏ ｄ ｓ ｅ ａ ｈｅ ｒ ａ ｉｔｏ ｌｐｒｎｃｐ ｅ；ｆｎｉｅ ｅ ｅ ｎｔ ｍ ｅ ｈｏｄ；ｓ ｚ ｆｅ ｔ ｉ ｔ ｌ ｍｅ ｔ ｉｅ ｅ ｆ ｃ ；ｎｕ—
Ｃｈ ｎ ｑ ｎ ０ ０ １ ｏ ｇ ｉ ｇ ４ ０ ４ ，Ｃｈ ｎ ；２ Ｓ ｈ ｏ ｆＣｉ ｉ Ｅｎ ｉ ｅ ｒｎ ｎ ｃ ｉ ｃ ｕ ｅ Ｃｈ ｎ ｓ ａ ｉ ａ ． ｃ ｏ ｌ ｖ ｌ ｇ ｎ ｅ ｉ ｇ ａ ｄ Ａｒ ｈ ｔ ｔ ｒ ， ａ ｇ ｈ ｏ ｅ Ｕｎ ｖ ｒ ｉｙ ｏ ｃｅ ｃ ｎ ｃ ｎ ｌ ｇ ｉ ｅ ｓｔ ｆＳ ｉｎ ｅ ａ ｄ Ｔｅ ｈ ｏ ｏ ｙ，Ｃｈ ｎ ｓ ａ ４ ０ ０ ａ ｇ ｈ １ ０ ４，Ｃｈ ｎ ；３ Ｙｕ ｆ ｇ ｗａ ｉ ａ ． ｎ ｕ Ｈｉ ｈ ｙ Ｐｌ ｎ ｉ ｇ ａ ｄ Ｄｅ ｉ ｎ Ｉ ｓ ｉｕ ｅ ａ ｎ ｎ ｎ ｓｇ ｎ ｔｔ ｔ ，Ｙｕ ｆ ０ ０ ｎ ｕ ７ ０ ０ ，Ｃｈ ｎ ） １ ｉ ａ
（． 放 军 后 勤 工 程 学 院 建 筑 工 程 系 ， 庆 １解 重 湖南 长沙 ４０４ ； ． 沙理工 大学 土木与建筑学 院 ， ００ １ ２长 ７０ ０ ） １ ０ ０ ４ ０ ０ ； ． 东 省 云 浮 市 公 路 勘 察 规 划 设 计 院 ， 东 云浮 １０ ４ ３ 广 广
摘
Ｖｏ ． ． １ ７ Ｎｏ ４
Ｄｅ ． ｃ ２Ｏ１ Ｏ
文 章 编 号 ： ６ ２ ９ ３ （０ ０ ０ － ０ ３ － ０ １ ７ — ３ １２ １ ）４ ０ ７ ６
Ｃ ｓｅａ 理 论 的有 限元 实现 及 ｏｓｒｔ 微 梁 弯 曲的尺 寸 效 应 模 拟
赵 冰 ｋ ，郑 颖 人 ，陈 淳
Ａｂ ｔａ ｔ ｓ ｒ ｃ ：Ｍ ｉｒ ｎ ｓ ａ ｅ ｓｚ ｆｅ ｔｉ ｆ ｅ ｎ ｏ ｎ ｅ ｅ ｎ ｔ ｅ ｄ ｓｇ ｎ ｎ ｆ ｃ ｕ ｅ ｏ — ｃ ｏ － ｃ ｌ ｉｅ ｅ ｆ ｃ ｓ ｏ ｔ ｎ ｅ ｃ ｕ ｔ ｒ ｄ ｉ ｈ ｅ ｉ ｎ ａ ｄ ｍａ ｕ ａ ｔ ｒ ｆｍｉ
关 键 词 ： ｏｓｒｔ 论 ； 分 原 理 ；有 限 元 ；尺 寸 效 应 ； 值 模 拟 Ｃ ｓｅａ 理 变 数
中图分类号 ： ＴＵ４ ２ ５ 史献 标 识 码 ： Ａ
Ｆｉ ｔ ｌｍ ｅ ｅ ｈ ｄ ｉ ｐｌ ｍ ｅ ｔ ｔｏ ｆ Ｃｏ ｓ ｒ ｔｔ ｅ ｒ ｎｉｅ ｅ ｅ ｎｔｍ ｔ ｏ ｍ ｅ ｎ ａ ｉ ｎ ｏ ｓ ｅ ａ ｈ ｏ ｙ
第 ７ 第 ４期 卷
２１ ０ ０年 １ ２月
长 沙 理 工 大 学 学 报 ห้องสมุดไป่ตู้自 然 科 学 版 ）
Ｊ ｕ ｎ ｌｏ ａ ｇｈ ｉｅｓｔ ｆＳ ｉｎ ｅａ ｄＴｅｈ ｏｏ ｙ Ｎａｕ ａ ｃｅ ｃ ） ｏ ｒ ａ ｆＣｈ ｎ ｓ ａＵｎ ｖｒｉ ｏ ｃｅ ｃ ｎ ｃ ｎ ｌｇ （ ｔｒ ｌＳ ｉｎ ｅ ｙ