考虑内部缺陷的FRP 拉索力学行为仿真分析

2022-10-09朱德举沈琰李晟郭帅成

湖南大学学报（自然科学版） 2022年1期

朱德举，沈琰，李晟，郭帅成

（1.湖南大学土木工程学院，湖南长沙 410082；2.绿色先进土木工程材料及应用技术湖南省重点实验室（湖南大学），湖南长沙 410082）

传统钢拉索材料存在自重大、耐疲劳和耐腐蚀性能较差等缺点，采用轻质、高强且耐腐蚀性能好的纤维增强树脂基复合材料（FRP）代替钢材制备大跨径斜拉桥的拉索组件可以有效弥补上述缺陷.众多学者［1-5］通过试验对各类FRP 材料的力学性能进行了充分测试，包括微观尺度纤维丝性能、宏观尺度单筋、拉索力学性能以及FRP-混凝土结构性能方面的表征.目前对于FRP 材料的研究通常是从宏观角度进行，忽略了材料内部复杂的原始缺陷构造，因而难以从本质上揭示材料变形与破坏的机制.为了获得大拉力FRP 拉索的力学性能评价和提升方法，需要采用多尺度力学理论，并利用仿真模拟结合试验测试的方法来综合分析.

FRP 拉索微观尺度上的缺陷主要是纤维丝内部孔洞缺陷.Eshelby 等效夹杂理论［6］、自洽理论、Mori-Tanaka方法［7］等力学模型结合有限元计算可以高效、高精度地预测含孔洞纤维丝的等效弹性模量.Sun 等［8］采用代表性体积元（Representative Volume Elements，RVE）预测单向纤维复合材料力学性能.常岩军等［9］通过扫描电子显微镜（SEM）对2.5D-C/SiC 复合材料微观结构进行观测，并利用M-K 方法计算了含不同孔洞形状和体积的复合材料的弹性模量，结果表明孔洞含量对弹性模量有显著影响.汪海滨等［10］通过SEM观察到的气相渗透工艺制备下的C/C-SiC 微观孔隙和微裂纹几何信息，建立了含有4 种缺陷的有限元模型，并计算了其等效弹性模量.任明法等［11］在纤维-基体尺度和复合材料单层板尺度上进行等效处理，分析了孔洞含量和孔洞形状对单向纤维增强复合材料单层板的等效弹性模量的影响.汪海滨等［12］分别建立了纤维丝、纤维束和宏观构件尺度有限元模型，并结合能量法研究了不同尺度下的关联参数对C/C-SiC 复合材料宏观力学性能的影响.

FRP 拉索宏观尺度的力学性能也是当前复合材料领域研究的热点.彭哲琦等［13］建立了以纤维浸胶纱为基本组成单元的单筋和FRP 拉索多尺度模型，并进行单筋和整索的强度分布规律预测.汪昕等［14］基于功能原理推导了FRP 锚固体系的锚固力计算公式，并基于参数分析发现，变刚度荷载传递介质的锚固体系可以有效避免大吨位FRP 拉索中的应力集中问题.Peng 等［15］进一步基于剪滞理论，以纤维浸胶纱和树脂基体为基本单元建立了FRP 拉索的随机有限元模型，模拟分析结果表明，其在强度预测和破坏形态分析方面的预测精度高于传统以纤维单丝为基本单元的FRP拉索有限元模型.

本文通过扫描电子显微镜（SEM）观测FRP 在微观尺度上的原始缺陷形态，对含原始缺陷的FRP 拉索进行了多尺度随机破坏有限元模拟.利用Python-ABAQUS 二次开发功能，针对含孔洞缺陷的纤维丝（微观）建立了代表性体积元有限元模型，基于多尺度（微观-细观-宏观）同一本构模型与渐进损伤理论，通过自定义子程序建立了宏观尺度FRP 单筋和拉索的有限元模型.

1 含内部缺陷FRP拉索有限元模型

由于纤维丝强度的随机差异性以及加工误差，FRP 的微观结构表现出随机性，其中包括：1）不同纤维丝的几何尺寸和缺陷数目存在差异，力学性能不完全相同；2）不同区域纤维丝的分布和间距具有随机性；3）不同区域纤维丝与界面的粘结性能存在差异.本文采用场发射扫描电子显微镜（SEM）对玄武岩纤维增强复合材料（BFRP）的横截面进行了表征（如图1 所示）.可以发现BFRP 微观尺度上的内部缺陷主要有两种形式：1）纤维丝内部的微孔洞缺陷；2）纤维丝与基体之间的界面缺陷.通过大量试样的电镜观察和统计分析可发现FRP 在内部不同区域孔洞等缺陷的体积分数不超过15%（通常2%左右，取决于样品的制备品质），而不同位置的纤维丝与基体之间界面裂隙缺陷面积差异较大.本研究中主要研究孔洞缺陷对于纤维弹性模量的影响规律.

图1 玄武岩纤维增强树脂基复合材料（BFRP）微观结构Fig.1 Micro-structure of basalt fiber reinforced polymer（BFRP）composites

1.1 含孔洞缺陷的纤维丝微观有限元模型

针对BFRP筋材（如图2（a）所示），本文选取含孔洞缺陷的纤维丝的代表性体积元（Representative Volume Element，RVE）建立微观有限元模型.通过模拟BFRP筋材的单轴拉伸力学行为，探讨其内部孔洞长径比、尺寸以及空间分布位置对纤维丝等效弹性模量的影响.

FRP 筋模型中内含孔洞缺陷的纤维丝满足以下假设［10］：1）纤维丝均具备相同圆形截面，周期性排列并且互相平行（如图2（b）所示）；2）FRP 单筋为横观各向同性，且泊松比为零；3）纤维丝内孔洞均为单方向排列.

基于ABAQUS 软件通过Python 语言进行二次开发建立上述几何和单元模型（如图2（c）和（d）所示），定义孔洞长径比为α=c/a，c为椭球的长轴，a为椭球的短轴.本研究选取了3 种代表性孔洞形状：1）球状孔，α=1；2）扁椭球状孔，α＜1；3）长椭球状孔，α＞1.BFRP 筋的直径为4.45 mm，密度为2.11 g/cm3；玄武岩纤维丝截面为圆形，直径约为13 μm.BFRP 筋的纤维体积分数约为60 %，纤维丝间距约为16.6 μm.Li［16］的研究表明，相对于矩形和棱形的RVE，截面为六边形的RVE模型更能代表大多数复合材料的统计特征，尤其是具有横观各向同性特点的材料.基于上述研究成果，本研究中选用了六边形RVE 模型进行进一步分析，模型基本尺寸参数如表1 所示，模型单元类型为C3D8R，采用六面体单元对网格进行划分，单根纤维丝含有9 600 个单元，基体划分为50 786 个单元.其中玄武岩纤维单丝和树脂基体的基本力学参数依据欧云福［17］的试验结果确定，如表2所示.单轴拉伸加载时，将RVE 模型一端固定，另一端沿长度方向进行位移加载（0.02 mm），共设置20个加载步.

图2 含孔洞缺陷BFRP筋的微观尺度有限元模型图Fig.2 Micro-scale finite element model of BFRP tendon with hole defects

表1 BFRP筋的RVE模型基本尺寸参数Tab.1 Basic dimensions parameters of the RVE model for BFRP tendon

表2 玄武岩纤维单丝和树脂基体的基本力学参数表Tab.2 Basic mechanical parameters of the basalt fiber and epoxy resin

1.2 多尺度渐进损伤统一本构模型

通过给微观RVE 模型施加位移荷载，得到单丝的应力-应变关系，通过均质化模型计算得到细观纤维束模型的平均应力、应变和弹性模量；细观模型同理可以将信息传递到宏观模型，为预测单筋和拉索的力学性能提供参数支持.

通过微观RVE模型和边界条件可以计算出同时含孔洞缺陷的微观应力和微观应变的分布［18］，将纤维束的细观应力与细观应变定义为［19］：

细观本构关系可以定义为：

对RVE 模型施加位移载荷，即可得到微观尺度上的应力和应变，再通过式（1）可以得到纤维束模型的细观应变和细观应力，其中式（1）的计算采用数值积分的方式进行：

式中：e表示RVE 单元号；Ne表示RVE 单元总数；和是RVE 单元e的平均应力和平均应变；Ve表示RVE单元e的体积.

基于上述假设，含孔洞纤维丝在横截面上各向力学性能相同，纤维丝看作是横观各向同性材料，同时纤维平行排布不考虑互捻，因此宏观的FRP 单筋和拉索同样为横观各向同性.式（2）中的本构模型为：

本构模型共有6 个独立的弹性常数，泊松比为0，纤维丝受力的方向为1 方向，纤维单丝横截面上相互垂直的两个方向为2和3方向.

构成FRP 单筋和拉索的纤维束的抗拉强度均服从Weibull 分布，即当单筋或拉索受到拉伸荷载作用时，抗拉强度较低的纤维束首先发生断裂破坏并丧失承载能力，所有荷载由余下纤维束继续承担，直至所有纤维束均达到自身抗拉强度破坏后，单筋和拉索不再承担荷载作用.

根据二参数Weibull 分布确定的纤维束强度分布如下：

通过ANSYS 自定义USERMAT 建立式（2）所示的本构模型，其中相关材料参数通过实验测得.微观尺度的孔洞分布会影响细观尺度RVE 和宏观尺度FRP 拉索的力学性能，进而影响Weibull 参数的取值.首先，将尺度参数（σ0）和形状参数（m）代入式（5）得到满足Weibull 分布的材料属性并赋给纤维束；然后在USERMAT 中实现刚度折减.具体作法为：通过TB.USER命令读入主程序中的单元类型、当前应力σ、当前应变ε以及根据指定位移和荷载步换算得到的应变增量Δε等基本参数，计算纤维束当前应变下的应力增量Δσ=D×Δε，并且更新应力σ=σ+Δσ；当纤维束的当前应力达到基于Weibull 分布赋予的抗拉强度后，其弹性模量折减为原数值的万分之一，纤维束失效破坏［20］.用户子程序和主程序之间的数据传递和运算流程如图3所示［21］.

图3 渐进损伤计算流程图Fig 3 Flow chart for the calculation of the progressive damage

1.3 传统M-T基本理论[22]

在M-T 方法中，将纤维丝定义为“0”相，孔洞定义为“1”相.现设纤维丝D域内有许多孔洞，c0为纤维丝的体积分数，孔洞的体积分数为c1，弹性模量为E1，受均匀应力σ作用，同时引入一个形状相同的均质无孔洞纤维丝作为对比，其弹性模量为E0，即无孔洞纤维丝的弹性模量为E0.基于M-T 基本理论和Eshelby 变换，可以得到带缺陷纤维丝与均值无孔洞纤维丝弹性模量之间的关系：

1.4 基于随机破坏的含缺陷单筋和拉索的宏观有限元模型

FRP 拉索由37 根FRP 单筋组合而成，FRP 单筋由若干FRP 纤维束构成，FRP 纤维束可以看作是由六边形RVE 按一定周期排列组合而成.RVE 有纤维丝、基体、界面相和内部孔洞.利用八节点实体单元Solid185 建立圆柱体模型来模拟BFRP 单筋（如图4所示）.单筋模型的长度为100 mm，玄武岩纤维单丝的直径为13 μm，横截面积为132.67 μm2；而BFRP单筋的有效横截面积为15.55 mm2，大约含有70 680 根纤维丝.为了提高计算效率，将单筋划分为48 根纤维束，其等效半径和等效横截面积分别为0.32 mm和0.32 mm2，每根纤维束内含有2 454 根纤维丝.纤维束的弹性模量由随机缺陷RVE 模型计算得到.单筋模型的加载过程为一端固定，另一端施加轴向拉伸位移荷载（0.02 mm），设定100个加载步.

图4 纤维束尺度的BFRP单筋有限元模型Fig.4 FE model of the single BFRP tendon at the scale of single yarn

2 结果与讨论

2.1 孔洞缺陷对纤维丝等效弹性模量的影响

由于纤维丝在受到拉伸作用时，故纤维丝等效横截面积的大小决定了该纤维丝强度的大小.定义受力方向上的孔隙率η为：

式中：Ah为垂直于受力方向截面的孔洞面积之和；Af为纤维丝横截面积.

图5 为孔洞尺寸与空间分布位置对纤维丝等效横截面积的影响.当孔隙率一定时，模拟结果的离散性随着孔洞半径的增大而增大，尤其是对含扁椭球状孔洞的纤维丝而言，当缺陷在任一横截面上分布得越离散，η越大，纤维丝等效弹性模量越小.如表3所示，控制孔隙率分别为1%、5%、10%、15%时，纤维丝的内部孔洞形状分别为三种尺寸的椭球状、三种直径的圆球状和四种尺寸的扁椭球状.其中在椭球中，a和b均为赤道半径，c为极半径.为了减小孔洞位置分布随机性造成的影响，取两次模拟结果的平均值，并与M-T 理论分析的结果进行比较（如图6 所示）.可以发现，随着孔隙率的增大，纤维丝的等效弹性模量降低，下降趋势与M-T 理论计算结果相同.由于理论计算的方法没有考虑孔洞在受力方向投影位置上重叠的可能性，因此数值模拟结果偏高.

图5 含不同孔洞尺寸纤维丝的等效弹性模量（孔隙率=10%）Fig.5 Equivalent elastic modulus of fibers containing internal holes with different sizes at 10% porosity

孔洞的横截面积是影响纤维丝等效弹性模量最直接且最重要的因素，图7 对比了在同一孔隙率下，分别改变孔洞横截面积和孔洞体积时，纤维丝等效弹性模量的变化趋势.如图7（a）所示，当孔洞体积不变时，横截面积越小的孔洞，等效弹性模量值越大；随着孔隙率的增加，孔洞横截面积大小对纤维丝等效弹性模量的影响越显著.如图7（b）所示，孔洞横截面积相同的条件下，当孔隙率低于10%时，孔洞体积的大小对等效弹性模量几乎无影响；当孔隙率高于10%时，由于孔洞的体积越大，孔洞个数会相对越少，纤维丝的等效弹性模量越大.

图7 孔洞横截面积与体积对纤维丝等效弹性模量的影响Fig.7 Effects of cross-sectional area and volume of holes on equivalent modulus of elasticity of the single fiber

由M-T 理论可知，长径比的大小会影响弹性模量随孔隙率增大而降低的变化程度.本文选取了5组不同的长径比（0.25、0.5、1、2、4）进行模拟分析，并将模拟结果与M-T 理论计算结果进行比较（如图8所示）.随着长径比的增大，孔洞形状越趋向于“裂缝”形，等效弹性模量下降速度趋缓并接近于直线.长径比越小，孔洞形状越趋向于“裂片”形，在孔隙率较低时，等效弹性模量下降速率很快.随着孔隙率的增加，下降速度趋于缓和.因此，扁椭球状孔洞缺陷的存在对纤维丝的强度影响更为显著.

图8 孔洞长径比和孔隙率对纤维丝弹性模量的影响Fig.8 Effects of slenderness ratio of holes and porosity on the elastic modulus of the single fiber

2.2 含孔洞缺陷纤维丝的应力分布

在相同孔隙率（5%）时，含5 种不同长径比孔洞的纤维丝在受拉过程中von Mises 应力云图如图9（a）所示.可以看出，含孔洞纤维丝在拉伸过程中，孔洞周边区域出现了显著的应力集中，距离孔洞较远处应力分布趋于均匀，应力最小值在孔洞沿受力方向上的两个端点处.当孔洞集中在纤维丝的某一纵向截面时，在此截面处孔洞同时断裂的可能性提高，材料的脆性相应增大.在孔隙率分别为1%，5%，10%和15%时，含半径为0.794 μm 圆球状孔洞的纤维丝拉伸过程中的von Mises 应力云图如图9（b）所示.随着孔隙率的升高，同尺寸的孔洞数量虽然增多，但应力分布反而更均匀；在孔隙率较高时，纤维丝沿着孔洞集中区域断裂的可能性较高.

图9 含不同长径比孔洞以及不同孔隙率的纤维丝的应力云图Fig.9 The von mises stress contours of fibers containing holes with different slenderness ratios and porosities

2.3 不同纤维束Weibull 参数与弹性模量下的宏观尺度单筋、拉索模拟

由于纤维束的弹性模量受缺陷数量和形式的影响，同时纤维束尺度Weibull 参数的选取也会影响单筋受拉时纤维束断裂的随机性.因此，为了得到符合试验结果的单筋本构模型，表4 列举了单筋模型在不同纤维束Weibull 参数与弹性模量下的模拟结果和Peng 等［15］的试验结果.将不同参数单筋模型拉伸的应力-应变曲线与文献［15］试验结果进行对比（图10 和表4），可以看出，当纤维束弹性模量Ef=55 GPa、Weibull 形状参数mf=40、尺度参数σf=1 500 MPa 时，单筋有限元模拟结果（FEM2）与试验结果吻合最好.单筋模型预测的弹性模量和抗拉强度与试验结果相比，误差分别为5%和2%.

图10 模型预测BFRP单筋应力-应变曲线与试验结果对比Fig.10 Comparison of the model predicted BFRP tendon stress-strain curves with the test results

表4 不同纤维束Weibull参数与弹性模量下的单筋模型模拟结果Tab.4 Simulation results of BFRP tendon model with different fiber bundle Weibull parameters and elastic modulus

以FEM2的纤维束模型为例进行渐进损伤分析，其应力-应变变化过程如图11 所示.刚开始加载时所有纤维束共同承担外部荷载作用，呈现线性应力应变关系；当应变达到0.022 mm/mm 时，出现第一根纤维束断裂，断裂的纤维束不再承担荷载退出工作，单筋整体应力为1 220 MPa，单筋损伤率约为1%，全部荷载由剩下47 根纤维束继续承担；随着应变的进一步增加，纤维束开始渐进破坏，应变为0.027 mm/mm 时，单筋达到抗拉强度极限值1 432 MPa，此时有4根纤维束发生断裂，单筋的损伤率约为7.4%；当应变达到0.030 mm/mm时，断裂纤维束达到45根，单筋的损伤率为93%.

图11 拉索试验结果与模型预测应力-应变曲线Fig.11 Stress-strain curve of FEM results and test results for BFRP cable

FRP 拉索由37 根直径为4 mm 和长度为100 mm的单筋组成，单筋相互平行.拉索模型的单元类型与单元划分、边界条件等与单筋模型相同.拉索模型预测的应力-应变曲线与Wang等［23］的试验结果对比如图10 和表5 所示.可以看出，当纤维束弹性模量Ef=55 MPa、Weibull 形状参数mf=20、尺度参数σf=1 500 MPa 的拉索有限元模拟结果（FEM1）与试验结果吻合良好.拉索模型预测的弹性模量和抗拉强度与试验结果（平均值）相比，误差分别为3.8%与8%.

表5 不同纤维束Weibull参数与弹性模量下的拉索模型模拟结果Tab.5 Simulation results of BFRP cable model with different Weibull parameters and elastic modulus of fiber bundle

以FEM7组拉索模型为例进行渐进损伤分析，其破坏过程为：当应变为0.01 mm/mm 时，拉索内所有纤维束的应力均为595 MPa，所有纤维束共同承担拉伸荷载；当应变达到0.019 mm/mm 时，出现纤维束断裂，此时拉索的应力为1 000 MPa；当应变为0.024 mm/mm 时，拉索的抗拉强度达到1 281 MPa，此时拉索内部损伤率为6.9 %.随着应变进一步增加，纤维束继续发生渐进破坏，在应变为0.026 mm/mm 时，50 %的纤维束断裂；当应变为0.027 mm/mm 时，1 527 束纤维束发生断裂，损伤率达到86%，未断裂的纤维束应力接近1 639 MPa.