复合材料雷击放电效应三维数值模拟

王富生，姬尧尧，刘志强，岳珠峰，张庆茂

（1．西北工业大学力学与土木建筑学院，西安 710129；2．中国航空工业集团公司成都飞机设计研究所，成都 610041）

复合材料雷击放电效应三维数值模拟

王富生1，姬尧尧1，刘志强1，岳珠峰1，张庆茂2

（1．西北工业大学力学与土木建筑学院，西安 710129；2．中国航空工业集团公司成都飞机设计研究所，成都 610041）

建立长空气间隙及复合材料层合板三维有限元模型，采用棒－板长间隙正极性流注生长概率模型为基于经典流注理论的随机放电模型。计算空气间隙区域内各网格点电位，列出电极周围所有流注待发展路径，计算每条流注待发展路径的生长时间，确定具有最小生长时间路径为优先选择路径。流注到达复合材料层合板后选流注与板交叉点作为雷击附着点进行复合材料层合板雷击直接效应分析，获得复合材料层合板在雷电压作用下电势、温度及热应力分布。计算结果表明，流注瞬间产生的高电势、温度及热应力主要沿顶层电导率最大方向对称扩展，所研究的初步结果可作为定性描述；为准确模拟复合材料的雷击放电效应过程，需进一步考虑复合材料的雷击烧蚀及热力学等破坏机理。

雷击；复合材料层合板；电势；温度；热应力

复合材料因力学性能、减重效果及防腐性能优良广泛用于飞机结构设计，但与金属材料如铝合金、钛合金相比电导性能差、对雷电较敏感，遭遇雷击可能性较大，尤其在雷击效应直接作用下更易遭受损伤，大幅降低复合材料结构的刚度及强度［1－3］。

雷电对飞机的初始附着过程普遍认为是长间隙放电［4］，Renardieres等［5－6］结合试验观测大量放电现象，对典型电极配置的正、负极性冲击电压作用长间隙放电基本物理过程有一定认识。棒－板长间隙正放电击穿的物理过程经历初始电晕、先导发展、最后跃变及主放电阶段；负放电过程较正放电复杂，在先导流注前方出现空心柱及向两边同时发展的正、负空间先导及流注，当负极性流注抵达平板后会形成向上发展的正极性迎面先导。基于文献［5－6］提出的长间隙放电模型中以Hutzler模型［7］、Gallimberti I模型［8］及Cooray V模型［9］最具代表性，其共性为均需计算空间电荷。Hutzler模型计算空间电荷时基于半经验公式，一定程度上限制该模型的应用；Cooray V模型用模拟电荷法计算空间电荷时因有限个模拟电荷难以反映空间电荷的复杂分布，较大误差亦限制该模型的应用；Gallimberti I模型结合经典气体放电理论与近代等离子体理论，能在一定程度上符合试验结果，为较好的长间隙放电物理模型。通过对长间隙放电击中点概率分布及影响因素试验研究［10］表明，长间隙放电路径具有一定随机性，而此特征Gallimberti I模型较难描述。因此，棒－板长间隙正极性流注生长概率模型［11］获得应用。该模型为结合经典流注理论的随机放电物理模型；但其研究仅停留在二维基础上，未用于飞机结构的雷击放电仿真，亟待建立理论三维模型用于飞机复合材料结构的雷击效应评估。虽已有研究从雷击热效应角度对复合材料的温度分布进行数值分析［12］，但雷击放电载荷强制施加于复合材料结构，忽略棒－板长间隙的放电过程，具有一定误差。

1 合材料放电流注生长仿真

1.1 仿真流程

建立经典流注理论的三维棒－板间隙放电物理模型见图1，基于经典的流注起始判据计算棒－板间隙流注起始电压及空气间隙区域内各网格点电位，列出电极周围所有流注待发展路径，计算每条流注待发展路径的生长时间，确定具有最小生长时间路径为优先选择路径。当流注到达复合材料层合板后选出雷击附着点进行复合材料层合板雷击直接效应分析。具体实现过程为：据边界条件求解空气间隙各网格点电位分布，棒电极施加电压波形，电晕起始条件采用Peek判据，当棒电极头部电场强度大于起始临界电场强度时判定为放电电子崩起始，取起始临界电场场强Ec＝28．5 kV／cm［8］；当空间电场具有足够支持二次电子崩发展能力时，部分电子崩将转化为等离子体通道后初始电晕起始，起始条件可表示为N≥Nc。其中N为正离子数目，Nc为初始电晕起始时的正离子数目，Nc＝0．55×108［8］。在建立的三维放电模型中规定棒电极头部正下方点为初始放电点，遍历已发展为流注空间点，对其周围26个方向满足局部场强大于临界场强且电离区中出现引发电子崩的二次电子点作为所选待发展点，见图2。

图1 棒－板间隙放电模型Fig．1 Dischargemodel of rod-plate gap

图2 待发展点选择示意图Fig．2 Sketch map of developed point

流注形成时间表达为

式中：ξ为均匀分布的0～1间随机数；r（E）＝k（E／Ec）η为流注生长概率函数；E为发展点电场强度；Ec为空气临界电场强度；k＝5×105为时间常数；η＝3为发展概率指数。

计算所有待发展点的流注形成时间，选具有最小时间待发展点作为流注将要发展的空间点。若周围有两个以上旧流注点，则选与新流注点之间场强最大的流注点作为新流注点的父节点，判断流注是否到达板电极，若未到则重新循环，直至流注发展到复合材料板判定为放电结束。流注到达复合材料层合板后选出流注与板的交叉点作为雷击附着点。整个流注生长分析流程见图3，可以此编写ANSYS的APDL程序用于复合材料结构的三维雷击放电数值模拟。

图3 流注生长分析流程图Fig．3 Analysis process of streamer growth

1.2 有限元模型

用ANSYS有限元软件中热电耦合实体单元SOLID5模拟复合材料层合板、实体单元SOLID122模拟空气间隙，施加高电压通过空气击中复合材料层合板进行雷电效应模拟，在给定的电边界条件及热边界条件下获得复合材料板的电势分布，通过瞬态热分析获得复合材料板的温度及热应力分布。

复合材料板放电有限元模型及复合材料板有限元模型见图4、图5。棒电极施加雷电操作波，具体模拟时高电压载荷施加在上方空气中心。棒电极下部与板间之距为200 mm，复合材料板与地相接，电势为0。空气间隙顶面及侧面电势为0，相对电导率为1。复合材料板顶面、侧面采用热传导第三类边界条件，规定层合板边界与周围流体间换热系数及周围流体温度；底面绝热采用热传导第二类边界条件，规定边界上热流密度值为0。考虑热辐射作用，热辐射率为0．9。空气温度为25℃，底面、侧面电势为0。

图4 复合材料板放电有限元模型Fig．4 Discharge FEmodel of composite plate

图5 复合材料板有限元模型Fig．5 FEmodel of composite plate

图6 电压A波形Fig．6 A voltage wave

雷电操作波用飞机结构雷击直接效应电压A波形［13－14］，见图6。由图6看出，波形上升率为1 000 kV／μs（±50%），增加的幅值直到击穿中间空气介质，导致试验件被击穿或闪络滑过，此时电压跌至零。击穿未发生时（雷击电压发生器为开环电压）电压跌落率及电压衰减时间不作规定。

采用碳纤维／环氧树脂基复合材料层合板IM600／133［12］，尺寸为400 mm×400 mm×5 mm，层数为8，铺层方式为［45／90／－45／0］s，材料属性见表1。

表1 材料属性Tab.1 Com positematerial proper ties

2 计算结果和分析

2.1 放电路径

典型时刻放电形态放大见图7，放电形态的整体效果见图8。仿真获得初始电晕起始电压为4 500 kV，初始电晕起始时刻为t＝4．5μs。初始电晕起始后随电极电压增加注入流注区域的空间电荷亦增加，流注进入持续发展阶段；t＝9．7μs时流注出现跃变到达复合材料板电极，此时截断电压为9 700 kV。雷电流到达复合材料层合板与其交叉点即为雷击附着点，恰为有限元节点9 583处，坐标为（0．21m，－0．006 7 m，0．2 m）。高电压击穿空气间隙后放电通道变为等离子通道，电导率极大提高，对复合材料结构产生直接效应。

2.2 复合材料板电势分布

复合材料板整体电势分布云图见图9，每层电势分布云图见图10。由图10（a）看出，雷电流沿面内45°方向传导后向边缘释放。层合板电势分布区域呈轴对称形式，表明雷电流传导沿顶层电导率最大方向，其它两方向电导率对雷电流传导影响较小。由图10（b）看出，电导率最大方向为90°，因厚度方向电阻率大、电流传导时间短雷电流在顶层瞬间扩散，导致该层电势较小。同理可知其它六层电导率最大方向各不同，导致各层电势分布不同；由于沿板厚方向电导率较小导致电势逐层减小。

图7 典型时刻流注生长形态Fig．7 Streamer growth form at typical time

图8 流注生长形态整体效果Fig．8 The whole effectmap of streamer growth form

图9 复合材料层合板整体电势分布Fig．9 The whole potential

图10 复合材料层合板各层电势分布Fig．10 Potential distribution of each layer for composite lam inate

复合材料板各层中心节点电势随时间变化曲线见图11。由图11看出，在第一、三、五、七层即奇数层中电势均为正值，而在第二、四、六、八层即偶数层中电势均为负值；随时间推移正、负电势绝对值均呈线性增长趋势。复合材料各层中心节点电势对比见图12，复合材料各层所有节点中最大电势对比见图13。从中可以看到：雷电流主要影响复合材料层合板第一和第二层，尤其对第一层影响最大。

图11 复合材料各层中心节点电势随时间变化曲线Fig．11 Time history of potential at the center point of each layer for composite laminate

图12 最后时刻板各层中心节点的电势Fig．12 Potential at the center point of each layer at last time

图13 最后时刻板各层所有节点中的最大电势Fig．13 Themaximum potential for all points of each layer

图14 复合材料层合板整体温度分布Fig．14 The whole temperature distribution of composite laminate

图15 复合材料层合板各层温度分布Fig．15 Temperature distribution of each layer for composite laminate

图16 复合材料各层中心节点温度随时间变化曲线Fig．16 Time history of temperature at the center point of each layer for composite laminate

2.3 复合材料板温度分布

复合材料层合板整体温度分布云图见图14、各层温度分布云图见图15。由图14看出，因流注向下作用时间极短、峰值电压较高，雷电流瞬间产生的高热量有少部分向板厚度方向传导，在雷电压附着点附近会产生板厚度方向的温度集中区域；由于焦耳热与电压平方呈正比，故温度分布与电势分布一致。由图15看出，层合板第一层温度集中区域沿电导率最大的45°方向分布，由于其它两方向电导率较小，致温度梯度变化较大；由于厚度方向电阻率较大又受热导率影响，致电流的热效应不明显，虽第二层90°方向与第一层45°方向温度分布近似，但数值相差较大，第三～八层温度变化更小，各层温度分布均近似呈椭圆形。

复合材料板各层中心节点温度随时间变化曲线见图16。由图16看出，随时间的推移温度不断增加，表现出较强非线性，且曲线变化趋势一致，均呈凹型。复合材料各层中心节点最大温度对比见图17、各层所有节点最大温度对比见图18。由图18看出，雷电流主要影响复合材料层合板的第一、二层，第一层影响最大，其它六层均较小；对复合材料各层中心节点，影响最大仍为第一层，但第三层影响较第二层大。

一般树脂基复合材料在超过200℃时刚度下降显著，而碳纤维／环氧树脂基复合材料的碳纤维升华温度为3316℃［15］，本文未考虑复合材料烧蚀、熔融及汽化所致能量消耗，致复合材料表层附着点附近温度较高。若反映实际情况，需考虑复合材料烧蚀等所致损伤问题。

图17 最后时刻板各层中心节点的最大温度Fig．17 Temperature at the center pointof each layer at last time

图18 最后时刻板各层所有节点中的最大温度Fig．18 Themaximum temperature for all points of each layer

图19 复合材料层合板整体热应力分布Fig．19 The whole heat stress distribution of composite laminate

图20 复合材料层合板各层热应力分布Fig．20 Heat stress distribution of each layer for composite laminate

2.4 复合材料板热应力分布

雷电流温度效应在复合材料内部产生的整体热应力分布云图见图19、各层热应力分布云图见图20。由图19看出，由于流注向下作用时间极短、峰值电压高，瞬间产生的高热应力沿顶层电导率最大方向扩展。由图20看出，除第一、二层分布近似呈椭圆形外，其它各层均表现出不同分布；铺层顺序对热应力分布影响较大，沿厚度方向铺层对称导致热应力分布基本对称。

复合材料板各层中心节点热应力随时间变化曲线见图21。由图21看出，随时间推移热应力不断增加，曲线变化形式各异，均呈台阶式发展，每条曲线均由不同斜率线段组成。复合材料各层中心节点最大热应力对比见图22，各层所有节点中最大热应力对比见图23。由两图看出，雷电流作用下复合材料的热应力主要影响层合板第一、二层，第一层影响最大，而其它层影响相差不大，影响最小为第六层；对复合材料各层中心节点影响最大仍为第一层，但对第二层影响与其它层相差不多，影响最小仍为第六层。由于各层热应力分布不均匀造成复合材料附着点附近热应力集中，易造成复合材料纤维断裂及基体破坏，层间热应力亦会引起复合材料分层。本文计算所得热应力远高于碳纤维／环氧树脂基复合材料极限应力，若需反映实际情况，需考虑复合材料的热应力损伤问题，包括纤维断裂、基体破坏及分层等。

图21 复合材料各层中心节点热应力随时间变化曲线Fig．21 Time history of heat stress at the center point of each layer for composite laminate

图22 最后时刻板各层中心节点最大热应力Fig．22 Heat stress at the center point of each layer at last time

图23 最后时刻板各层所有节点最大热应力Fig．23 Themaximum heat stress for all points of each layer

3 结 论

（1）推广正极性流注生长概率放电模型至三维空间可用于复合材料层合板雷击附着点获取及电势、温度、热应力定性评估。

（2）由于流注向下作用时间极短、峰值电压较高，瞬间产生的高电势、温度、热应力主要沿顶层电导率最大方向对称扩展。在不考虑复合材料烧蚀、熔融、汽化及热应力损伤情况下，雷击高电压主要影响复合材料前几层，对其它层影响较小。

（3）复合材料各层的电势、温度及热应力与时间的关系分别为线性、非线性及台阶式发展趋势。

（4）实际中复合材料达不到本文模拟的温度及热应力时即遭烧蚀、熔融、汽化或受力破坏，故为准确进行复合材料的雷击数值模拟，需通过研究复合材料的雷击破坏机理发展适合的复合材料烧蚀、熔融、汽化及热力学破坏准则。

［1］Baldacim S A，Cristofani N，Junior J L F，et al．Lightning effects in aircraft of the composite material［C］．17CBECIMat，2006：3247－3258．

［2］Paolo F，Mark M，Damage resistance and tolerance of carbon epoxy composite subjected to simulated lightning strike［J］．Composites，Part A，2009，40：954－967．

［3］Yoshiyasu H，Shingo K，Yutaka I，et al，Artificial lighting on graphite／epoxy composite laminate［J］．Composites，Part A，2010，41：1461－1470．

［4］马乃祥．长间隙放电［M］．北京：中国电力出版社，1998．

［5］Renardieres L．Positive discharges in long air gaps at Les Renardieres 1975 results and conclusions［J］．Electra，1977，53：131－152．

［6］Renardieres L．Negative discharges in long air gaps at Les Renardieres 1978 results［J］．Electra，1981，74：167－218．

［7］Hutzler B．Leader propagation model for predetermination of switching surge flashover voltage of large air gaps［J］．IEEE Transaction on Power Apparatus and Systems，1978，96（4）：1087－1096．

［8］Gallimberti L，Bacchiega G，Bondiou-Clergerie A，et al． Fundamental processes in long air gap discharge［J］．Applied Physics，2002，3：1335－1359．

［9］Becerra M，Cooray V．A simplified physical model to determine the lightning upward connecting leader inception［J］．IEEE Transaction on Power Delivery，2006，21（2）：897－908．

［10］D’Alessandro F，Kossmann C J，Gaivoronsky A S，et al，Experimental study of lightning rods using long sparks in air［J］．EEE Transaction on Dielectrics and Electrical Insulation，2004，11（4）：638－649．

［11］贺恒鑫，何俊佳，钱冠军，等，棒－板长间隙正极性流注生长概率模型及应用［J］．高电压技术，2008，34（10）：2047－2053．

HE Heng-xin，HE Jun-jia，QIAN Guan-jun，et al，Stochastic model of streamer growth in air for a rod-plate gap and its application［J］．High Voltage Engineering，2008，34（10）：2047－2053．

［12］Toshio O，Yoshiyasu H，Akinori Y．Coupled thermalelectrical analysis for carbon fiber／epoxy composites exposed to simulated lighting current［J］．Composites：Part A，2010，41：973－983．

［13］SAE－ARP－5412A，Aircraft lightning environment and related testwaveforms［S］．2005．

［14］段泽民，曹凯风，程振革，等，飞机雷电防护实验与波形［J］．高电压技术，2000，26（4）：61－63．

DUAN Ze-m in，CAO Kai-feng，CHENG Zhen-ge，et al．Lightning protection tests and waveforms for aircraft［J］．High Voltage Engineering，2000，26（4）：61－63．

［15］Griffisetal CA，Masumra R A，Chang C I，et al．Response of graphite epoxy composite subjected to rapid heating［J］．Composite Materials，1981，15：427－442．

3D numerical simu lation of com posite discharge effect under lightning strike

WANG Fu-sheng1，JIYao-yao1，LIU Zhi-qiang1，YUE Zhu-feng1，ZHANG Qing-mao2
（1．School of Mechanics Civil Engineering and Architecture，Northwestern Polytechnical University，Xi'an 710129，China；2．Chengdu Aircraft Design and Research Institute of Aviation Industry Corporation of China，Chengdu 610041，China）

A 3D FE model of long air gap and composite laminate was built．The positive streamer growth probabilitymodel of a long rod-plate gap was adopted as a stochastic discharge model based on the traditional streamer theory．The detailed simulation process includes：calculating the electric potential at each grid point in air gap，listing routes of all streamers developed around the pole，calculating the growth time of route developed for each streamer and defining the route with minimum growth time as the prior selected one．When the streamer reaches the composite plate，the cross point of the streamer and plate was regarded as the attachment point of lightning strike．A direct effect analysis was carried out to obtain the potential，temperature and heat stress distribution of composite plate can be obtained under the high voltage wave strike of lightning．The results show that the high potential，temperature and heat stress produced instantaneously by fast streamer mainly spread out symmetrically along the direction with maximum conductivity on top plate．The analysis results，provided as a quality description of the process，are only the preliminary．In order tomore accurately simulate the process of composite discharge effect under lightning strike，the damagemechanism of ablation and the thermodynamics for composite laminate need to be further taken into account．

lightning strike；composite laminate；potential；temperature；heat stress

TB332

：A

10．13465／j．cnki．jvs．2014．22．003

国家自然科学基金（51475369）；航空科学基金（2013ZF53068）；西北工业大学基础研究基金（JC20110257）

2013－07－24 修改稿收到日期：2013－11－21

王富生男，博士，副教授，1979年生邮箱：fswang＠nwpu．edu．cn