本发明属于cae仿真、工程结构力学计算领域，具体是指一种基于应力杂交化后处理的板壳边缘效应的计算方法。
背景技术：
板壳结构是一类非常重要的结构形式，被广泛地应用于航空航天、土建、船舶等工程领域。板壳结构边缘效应是一种典型的局部应力现象，是指在某些工况下自由边界邻近非常狭小的区域内，应力分布会出现非常剧烈的梯度变化。而这一狭小的区域称为板壳的边界层，尺度与板壳厚度同量级。在一些工程应用领域中，例如复合材料层合板(vidalp,polito,d′ottaviom,etal.assessmentoftherefinedsinusplatefiniteelement:freeedgeeffectandmeyer-pieningsandwichtest.finiteelementsinanalysisanddesign,2014,92(0):60-71.)、超大型浮式海洋结构(wangcm,xiangy,utsunomiyat,etal.evaluationofmodalstressresultantsinfreelyvibratingplates.internationaljournalofsolidsandstructures,2001,38(36-37):6525-6558.)，该区域内的应力峰值极有可能会导致板壳结构的破坏或失效，因此准确的求解边缘效应现象非常重要。
在绝大多数的工程应用中，板壳结构分析均会使用有限元方法进行求解。但是常规的有限元方法是一种基于分片插值的纯数值方法。受限于单元的插值阶次不足，其在处理边缘效应这类局部应力问题时表现十分乏力。针对边缘效应问题的现有技术大致可以分为两类。第一类是通过网格自适应技术来对板壳结构的边界层区域进行处理，通过减少边界层区域的网格尺寸来提高有限元计算的精度(schwabc,surim.thepandhpversionsofthefiniteelementmethodforproblemswithboundarylayers.mathematicsofcomputation,1996,65(216):1403-1429.)。这是比较常用的方式。但是由于边缘效应的作用范围远小于结构的整体尺寸，自适应计算的收敛速度会受到很大的影响。特别是在分析显式动力学问题时，边界层内过小的单元尺寸会严重降低整体计算的稳定时间步长，从而导致计算成本的显著增加。第二类则是使用特殊构造的有限元模型来处理边界层区域。例如ramesh等(rameshss,wangcm,reddyjn,etal.computationofstressresultantsinplatebendingproblemsusinghigher-ordertriangularelements.engineeringstructures,2008,30(10):2687-2706.)通过大幅度增加节点自由度数量来提高插值阶次。但是这种做法同样会导致高昂的计算成本。
本发明人也曾通过构造特殊的有限元模型来有效地求解该类问题(shangy,cens,licf,etal.aneffectivehybriddisplacementfunctionelementmethodforsolvingtheedgeeffectofmindlin-reissnerplate.internationaljournalfornumericalmethodsinengineering,2015,102(8):1449-1487.)(shangy,cens,liz,etal.improvedhybriddisplacementfunction(ihdf)elementschemeforanalysisofmindlin-reissnerplatewithedgeeffect.internationaljournalfornumericalmethodsinengineering,2017,111(12):1120-1169.)，但是其依赖于特定的有限元分析程序，其在通用性以及求解效率上均有明显不足。此外，该方法只能够用于分析各向同性材料的平板结构边缘效应问题。
技术实现要素：
本发明针对现有技术中存在的问题，公开了一种基于应力杂交化后处理的板壳边缘效应的计算方法，该方法是首先基于有限元软件强大的求解能力，利用求解结果对有限元网格模型进行应力杂交化后处理，本发明方法计算效率更高，能够求解更广泛的工程结构，具有更好的灵活性与适用性。解决了现有技术中存在的问题。
本发明是这样实现的：
一种基于应力杂交化后处理的板壳边缘效应的计算方法，其特征在于，步骤如下：
步骤一：建立的板壳结构模型：基于mindlin-reissner理论对实际的板壳结构进行简化，建立相应的力学模型；m-r理论考虑了板壳结构的横向剪切变形，能够有效地解释边缘效应存在的原因；因此在本方法中，具体选用mindlin-reissner板壳理论。
步骤二：针对步骤一建立的板壳结构力学模型；使用有限元软件对其进行求解，得到结构整体的位移、应变能、应力、应变；本发明中使用的有限元软件中商业的有限元软件均可，为方便也可以指定为具体的软件，如abaqus。)商业有限元软件具备强大的求解能力，该方法充分利用其强大的求解能力，计算效率更高，能够求解更广泛的工程结构，具有更好的灵活性与适用性。
步骤三：针对板壳结构发生边缘效应现象的部分，即自由边界的临近区域，提取出由步骤二计算得到的该区域的位移；然后对该区域进行应力杂交化后处理，从而获得边界层区域内的应力分布，并提取应力峰值与所出现的位置作为评估参数；
步骤四：对板壳结构的有限元网格模型减少网格尺寸，(可以将网格的特征尺寸减到一半或者根据实际的用户需要减少网格尺寸)增加网格密度，然后重复步骤二以及步骤三，得到一组新的应力峰值与位置；评估前后两次计算所得的结果是否满足用户预先设定的收敛准则以及或所需的计算精度；是否满足的判据是两次计算的相对误差低于某一个值。这个值通常由用户根据自己的要求确定，例如可取1％或更低。值越低，计算成本越高。
4.1，如若不满足，则进一步通过减少网格尺寸、增加网格密度来细化有限元网格模型，并重复步骤二到步骤四，然后再次进行收敛性判断，直至满足收敛准则以及所需的计算精度；
4.2，如若满足，即得到符合要求的板壳结构边缘效应的内力结果。
进一步，所述的步骤一具体为：
1.1，m-r理论考虑了板壳结构的横向剪切变形，能够有效地解释边缘效应存在的原因。使用mindlin-reissner板壳理论可以对板壳结构的厚度方向给予几何建模，从而将结构简化为一个二维结构，使用其中性层的二维拓扑来进行表征。
在mindlin-reissner理论中，板壳结构每一个点有6个自由度，包括3个平动自由度和3个转动自由度。将这些自由度可以表示为下述矢量形式：
q＝[uvwθxθyθz]t(1)
式中，q为自由度矢量，u,v,w为三个坐标方向的平动自由度，θx，θy，θz为三个方向的转动自由度。
根据变形模式特征，这6个自由度又可以分为两组，分别对应着面内变形部分：
qin＝[uvθz]t(2)
以及横向弯曲变形部分：
qout＝[wθxθy]t(3)
1.2，通过理论分析可知边缘效应局部应力现象主要需要考虑其弯曲变形模式，mindlin-reissner理论中与弯曲变形有关的广义应变分量和应力分量为别为
ε＝[κxκyκxyγxzγyz]t(4)
r＝[mxmymxytxty]t(5)
式中，e为广义应变分量，κx,κx,κxy为弯曲曲率，γxz,γyz为横向剪切应变；r为应力分量，mx,mx,mxy为弯矩和扭矩，tx,ty为横向剪力。
其中，扭矩与剪力分量因为边缘效应会在很小的范围内发生剧烈的梯度变化，应变与位移之间的得到：
κx＝-θx,x,κy＝-θy,y,κxy＝-θx,y-θy,x(6)
γx＝w,x-θx,γy＝w,y-θy(7)
应力需满足平衡方程：
mx,x+mxy,y-tx＝0,mxy,x+my,y-ty＝0,tx,x+ty,y+q＝0(8)
式中，q是横向分布载荷；
应变与应力之间的本构关系可以表示为
d是材料的弹性矩阵，具体表达需要根据板壳结构的具体材料决定,db是与弯曲变形相关的子矩阵，ds是与横向剪切变形相关的子矩阵。
进一步，所述的步骤二具体为：
2.1，使用有限元软件对板壳结构模型进行分析求解，得到整个结构的位移与应力结果。
2.2，根据结构的边界约束类型，甄别出其自由边界部分，确定需要进行应力杂交化后处理的区域；
2.3，提取2.2中所述区域内的网格信息以及节点位移结果，并从中分离出与离面弯曲变形有关的位移分量。
进一步，为了求解边缘效应，应力杂交化后处理主要是对弯曲变形部分进行。所述的步骤三具体为对弯曲变形部分进行应力杂交化后处理，过程如下：
3.1，应力杂交化后处理是对弯曲变形部分进行的，其核心思想是由杂交元技术衍生而来。首先设定一组弯曲应力场，即：
r＝sβ(10)
β是一组待定参数，s是一组应力试探函数；
3.2，根据步骤二得到的结构整体的位移qout设定边界位移模式：
d＝nqout(11)
选用梁的变形模式来近似模拟板壳结构的边界变形，由此可得n的具体表达；
3.3，上述弯曲应力场与边界位移在每个单元中需要满足最小余能原理，即
式中第一项为域内积分，ω表示控制域；第二项为边界积分，γ表示控制域的边界；
3.4，为了准确反应出自由边界处部分内力为零的特殊性，需要使用lagrange乘子法对式(13)进行修正，即将其变换为：
式中，
式中m、
m＝∫∫ωstd-1sdω(16)
hλ＝(γm-1γt)-1γm-1h(18)
h＝∫γstltndγ(19)
γ＝∫slsds(20)
l是方向余弦矩阵；
3.5，在求解式(15)后再将β代入到式(10)中即可得到后处理之后的应力结果；此外，对式(10)求导可以进一步得到各应力分量的峰值以及其所出现的位置；
3.6，边缘效应是一种典型的局部应力现象，在很小的尺度范围内，应力会出现非常显著的梯度变化。因此式(10)中的试探函数s必须具备有效模拟该现象的能力。本方法的一个关键点就是将具有下述指数函数形式的应力试函数扩充到式(10)中：
式中m，n是在切平面内该自由边界的方向余弦；k是与结构厚度以及材料属性相关的参数；
本发明与现有技术相比的有益效果在于：
相比较于现有技术中使用纯数值的有限元方法以及依赖于特定的有限元模型、有限元分析程序，本发明首先利用成熟的有限元分析软件进行整体结构的求解，再针对板壳结构发生边缘效应现象的部分进行杂交化后处理，根据是否满足收敛性需求进行精细化处理，最终得到符合要求的板壳结构边缘效应的内力结果。首先在分析程序求解边缘效应问题上，本发明的方法计算效率更高，能够求解更广泛的工程结构，具有更好的灵活性与适用性；此外，与通过网格细分求解边缘效应问题的常规做法相比，本发明的方法有效的克服了细分网格导致计算规模过大的缺点，显著地提高了计算效率与计算精度。
附图说明
图1为本发明的方法原理流程图；
图2为运用本发明方法求解边缘效应的实例示意图；
图3为运用本发明方法求解实例所得到的扭矩分量和剪力分量的分布曲线示意图。
具体实施方式
为使本发明的目的、技术方案及效果更加清楚，明确，以下参照附图并举实例对本发明进一步详细说明。应当指出此处所描述的具体实施仅用以解释本发明，并不用于限定本发明。
如图1所示，图1为一种基于应力杂交化后处理的板壳边缘效应的计算方法，主要分为4个部分：
第一部分：基于mindlin-reissner理论对实际的板壳结构进行简化，建立相应的力学模型用于求解；具体步骤为：
使用mindlin-reissner板壳理论可以对板壳结构的厚度方向不予几何建模，从而将结构简化为一个二维结构，使用其中性层的二维拓扑来进行表征。在m-r理论中，每一个节点有6个自由度，包括3个平动自由度和3个转动自由度。将这些自由度可以表示为下述矢量形式：
q＝[uvwθxθyθz]t(1)
式中，q为自由度矢量，u,v,w为三个坐标方向的平动自由度，θx，θy，θz为三个方向的转动自由度；
根据变形模式，这6个自由度又可以分为两组，分别对应着面内变形部分与离面弯曲变形部分，即
qin＝[uvθz]t(2)
qout＝[wθxθy]t(3)
边缘效应局部应力现象主要需要考虑其弯曲变形模式。根据m-r理论，与弯曲变形有关的广义应变分量和应力分量为别为
ε＝[κxκyκxyγxzγyz]t(4)
r＝[mxmymxytxty]t(5)
式中，e为广义应变分量，κx,κx,κxy为弯曲曲率，γxz,γyz为横向剪切应变；r为应力分量，mx,mx,mxy为弯矩和扭矩，tx,ty为横向剪力；
其中，扭矩与剪力分量因为边缘效应会在很小的范围内发生剧烈的梯度变化。应变与位移之间的得到：
κx＝-θx,x,κy＝-θy,y,κxy＝-θx,y-θy,x(6)
γx＝w,x-θx,γy＝w,y-θy(7)
应力需满足平衡方程：
mx,x+mxy,y-tx＝0,mxy,x+my,y-ty＝0,tx,x+ty,y+q＝0(8)
应变与应力之间的本构关系可以表示为
式中弹性矩阵的具体表达需要根据板壳结构具体材料与铺层来决定。
第二部分：对于步骤一中建立的结构模型，使用有限元软件软件，例如abaqus，对其进行整体求解，得到整个结构的位移。根据结构的边界约束类型，甄别出其自由边界部分，确定需要应力杂交化后处理的区域。提取该部分区域内的网格信息，以及从节点位移中提取出与离面弯曲变形有关的位移分量，即式(3)，以用于进一步的后处理。
第三部分：针对板壳结构发生边缘效应现象的部分，即自由边界的临近区域，提取出由步骤二计算得到的该区域的位移；然后对该区域进行应力杂交化后处理，从而获得边界层区域内的应力分布，并提取应力峰值与所出现的位置作为评估参数。具体步骤为：
应力杂交化后处理是对弯曲变形部分进行的，其核心思想是由杂交元技术衍生而来。首先假设一组符合要求的弯曲应力场，即
r＝sβ(10)
β是一组待定参数，s是一组应力试探函数。然后根据步骤二得到的位移qout假设一个合理的边界位移模式
d＝nqout(11)
具体选用梁的变形模式来近似模拟板壳结构的边界变形，由此可以的n的具体表达。上述假设的应力场与边界位移在每个单元中需要满足最小余能原理，即
为了准确反应出自由边界处部分内力为零的特殊性，需要使用lagrange乘子法对式(13)进行修正，即将其变换为：
通过推导可以得到待定参数β与步骤二所得位移qout之间的关系
式中m、
m＝∫∫ωstd-1sdω(16)
hλ＝(γm-1γt)-1γm-1h(18)
h＝∫γstltndγ(19)
γ＝∫slsds(20)
式中，l是方向余弦矩阵。在求解式(15)后再将β代入到式(10)中即可得到后处理之后的应力结果。此外，对式(10)求导可以进一步得到各应力分量的峰值以及其所出现的位置。
边缘效应是一种典型的局部应力现象，具有强烈的梯度变化。因此式(10)中的试探函数s必须具备有效模拟该现象的能力。对此需要将具有下述指数函数形式的应力试函数扩充到式(10)中：
式中m，n是在切平面内该自由边界的方向余弦；k是与结构厚度以及材料属性相关的参数；
第四部分具体为：
将板壳结构的有限元网格模型进行细化，减少网格尺寸(将网格的特征尺寸减到一半，或者根据用户实际需要减少网格尺寸)，增加网格密度，然后重复步骤二与步骤三，获得边缘效应一组新的内力峰值与峰值所在位置。评估前后两次计算结果的收敛性是否满足精度的要求。
若不满足要求，则再次对有限元网格模型进行精细化处理，并重复步骤二到步骤四的分析计算。若满足要求，则输出最终次的计算结果。
图2是一个边缘效应简单算例的示意图：边长为a的一个矩形方板承受横向均布载荷作用。左右边界为自由边，上下边界为简支约束。利用对称性，可以只取四分之一结构，即阴影部分，进行建模并分析。
如图3所示，由应力杂交化后处理得到的边缘效应扭矩分量和剪力分量分布曲线示意图，3a所示的是扭矩分量沿ab边分布；3b所示的是剪力分量沿ab边分布。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。