一种基于格构模型的混凝土材料损伤演变过程的模拟方法

1.本发明涉及一种基于格构模型的混凝土材料损伤演变过程的模拟方法，属于混凝土材料数值模拟领域。


背景技术：

2.作为土木工程领域中最常见的建筑材料，混凝土凭借其优良特性被广泛应用于民用建筑、桥梁、公路、隧道、人防工程以及水利工程等各种基础和民用工程建设当中。混凝土材料在浇筑过程中，由于水化反应放热导致温度升高，骨料的不规则性使混凝土材料在冷却成型之后内部就含有空隙、毛细孔等缺陷。在外部条件作用下，缺陷附近区域将出现应力集中，从而导致该区域损伤开始演化，微裂缝逐渐形成。在外荷载持续作用下，这些微裂纹将逐渐扩展、贯通直至最后形成宏观裂缝，最终导致混凝土构件失效破坏。因此，混凝土材料的失效破坏始于材料内部的损伤演化，研究该损伤演化过程有利于了解混凝土材料的失效破坏机理，对提高混凝土结构的强度具有重要意义。
3.目前针对混凝土材料损伤演化过程的数值模拟，多采用实体单元对混凝土试件进行数值建模，所建立的数值模型具有较多的自由度，因此具有较大的计算量，往往会导致计算效率偏低。


技术实现要素：

4.本发明提供一种基于格构模型的混凝土材料损伤演变过程的模拟方法，模拟整个混凝土材料的损伤演化过程，提高了计算效率。
5.本发明解决其技术问题所采用的技术方案是：
6.一种基于格构模型的混凝土材料损伤演变过程的模拟方法，具体包括以下步骤：
7.步骤s1：构建混凝土材料细观格构模型；
8.步骤s2：基于timoshenko梁理论建立混凝土材料细观格构模型的有限元方程；
9.步骤s3：基于混凝土材料的特性，采用mazars损伤本构模型模拟混凝土材料力学性能的劣化；
10.步骤s4：按位移加载的牛顿迭代法对混凝土材料细观格构模型进行计算求解，将计算结果经过处理对混凝土材料的损伤演化过程进行仿真分析；
11.作为本发明的进一步优选，步骤s1具体包括：
12.步骤s11：确定混凝土材料的尺寸，在待构建的混凝土材料细观格构模型边界上生成固定单元节点，在混凝土材料细观格构模型区域内随机生成单元节点，遍历所有单元节点，并按照由下至上、由左至右的顺序，将单元节点全部存储于节点信息矩阵内；
13.步骤s12：继续遍历所有单元节点，并按照由下至上、由左至右的顺序将所有单元节点依次相连，形成连接后的单元节点编号并存放在单元信息矩阵内，完成格构模型的几何搭建；
14.步骤s13：确定各级配骨料的粒径范围及个数，在混凝土材料细观格构模型区域内
生成各骨料级配的随机骨料，沿着骨料边界向外等厚度拓展界面过渡层，将生成的骨料信息存入骨料信息矩阵中；
15.步骤s14：结合步骤s12中搭建的格构模型，将单元中点的位置对步骤s13中的骨料进行识别，对于不同相材料的单元赋予相应的材料参数，即若单元中点的位置位于骨料区域内，则将此单元判定为骨料单元，若单元中点的位置位于骨料边界与界面过渡层之间，则将此单元判定为界面过渡层，若单元中点的位置位于界面过渡层区域外，则将此单元判定为砂浆单元；
16.作为本发明的进一步优选，步骤s2具体包括：
17.步骤s21：构建的混凝土材料细观格构模型选用timoshenko梁单元，并根据有限元理论推导出剪切刚度矩阵和弯曲刚度矩阵；
18.步骤s22：剪切刚度矩阵采用缩减积分代替精确积分，将剪切刚度矩阵和弯曲刚度矩阵叠加，获取单元刚度矩阵；
19.步骤s23：通过坐标转换矩阵将单元坐标系下的单元刚度矩阵转化为整体坐标系下的单元刚度矩阵，再将各个整体坐标系下的单元刚度矩阵中的元素按照编码规则叠加至总体刚度矩阵中，得到混凝土材料细观格构模型的有限元方程，其中总体刚度矩阵为变量，其与混凝土材料细观格构模型中单元的损伤程度有关联；
20.作为本发明的进一步优选，步骤s23中得到的混凝土材料细观格构模型的有限元方程为
21.k
·
q＝p
ꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
22.公式(1)中，k为总体刚度矩阵，q为位移向量，p为荷载向量；
23.作为本发明的进一步优选，步骤s3具体包括：
24.步骤s31：步骤s2中得到的混凝土材料细观格构模型有限元方程中，总体刚度矩阵与混凝土材料细观格构模型中单元的损伤程度有关联，采用mazars损伤本构模型模拟混凝土材料力学性能的劣化，通过损伤变量d来表示混凝土材料细观格构模型中单元的损伤程度；
25.步骤s32：基于lemaitre应变等效原理，得到混凝土材料细观格构模型中单元的刚度削减方程，通过刚度的不断削减体现混凝土材料整体力学性能的劣化；
26.作为本发明的进一步优选，步骤s31中，当d＝0，表示混凝土材料细观格构模型中单元无损伤；
27.当单元等效应变大于应变阈值时，混凝土材料损伤开始演化；
28.当d＝1，表示混凝土材料细观格构模型中单元失效破坏；
29.作为本发明的进一步优选，损伤变量d的表达公式为
[0030][0031]
公式(2)中，ε
p
为应变阈值，ε为单元等效应变，a
t
和b
t
为混凝土材料常数，对于timoshenko梁单元，其单元等效应变取为单元的最大主拉应变；
[0032]
混凝土材料细观格构模型中单元的刚度削减方程为
[0033]
[0034]
公式(3)中，e为单元无损状态下的弹性模量，为单元受损之后的弹性模量；
[0035]
作为本发明的进一步优选，步骤s4具体包括：
[0036]
步骤s41：对混凝土材料细观格构模型施加边界条件，在加载点处施加位移荷载并将其分为k个荷载步；
[0037]
步骤s42：在第k个荷载步内，给定初始位移增量dk，可计算出各格构单元应变εk[0038]
εk＝b
·dk
ꢀꢀꢀ
(4)
[0039]
公式(4)中，b为单元应变矩阵；将位移增量与单元应变分别存储于位移增量列阵与单元应变矩阵内；
[0040]
步骤s43：在第k个荷载步内，将各格构的单元等效应变与应变阈值进行比对，若单元等效应变大于应变阈值ε
p
，则代入公式(2)，得到单元发生损伤并更新损伤变量dk[0041][0042]
此时在第k个荷载步内，单元应力为
[0043]
σk＝(1-dk)
·d·
εkꢀꢀꢀ
(5)
[0044]
公式(5)中，d为单元无损状态下的单元弹性矩阵，εk为格构单元应变，且εk＝b
·dk
；将更新后的损伤变量dk与单元应力σk分别存储于损伤变量列阵与单元应力矩阵中；
[0045]
步骤s44：根据第k个荷载步内的单元节点内力与外力是否平衡来判断整个混凝土材料是否处于平衡状态，若结构处于平衡状态，则该荷载步计算收敛，更新所有单元的历史变量并进入第k+1个荷载步的迭代计算过程；
[0046]
此时由内应力产生的节点力必须与所施加的节点外力大小相等且方向相反，整个结构平衡的弱形式为
[0047][0048]
公式(6)中，为满足边界条件的任意虚位移，tk为第k个荷载步中的面力，fk为第k个荷载步中的体力；
[0049]
若整体结构未处于平衡状态，则节点内力与节点外力之间的残余力向量rk表示为
[0050][0051]
公式(7)中，n为单元形函数；
[0052]
步骤s45：在第k个荷载步中进行多次的循环迭代，当残余力向量小于预设的容差值时，则认为该荷载步内的迭代计算收敛并进入第k+1个荷载步的迭代计算过程，其中每一步的迭代公式为
[0053][0054]
公式(8)中，为第k个荷载步内结构的切线刚度矩阵，为该荷载步内的位移增量，i表示该荷载步中的第i次迭代；
[0055]
步骤s46：将所有荷载步的计算结果存储在结果信息矩阵中，计算完成后对结果信息矩阵内的元素进行后处理，模拟出混凝土材料的损伤演变过程。
[0056]
通过以上技术方案，相对于现有技术，本发明具有以下有益效果：
[0057]
1、本发明提供的模拟方法，建立格构模型，通过离散的timoshenko梁单元对混凝
土材料进行力学等效，对处于不同相材料位置的单元赋予相应的材料特性，直接体现混凝土各相材料的力学特性；
[0058]
2、本发明建立的混凝土材料细观格构模型将一个复杂的多轴问题转化为多个简单的单轴问题，即将混凝土材料划分成梁单元，保证模拟有效性的同时大大提升了计算效率。
附图说明
[0059]
下面结合附图和实施例对本发明进一步说明。
[0060]
图1是本发明提供的模拟方法全过程流程图；
[0061]
图2是本发明提供的实施例关于混凝土材料细观格构模型；
[0062]
图3是本发明提供的mazars损伤本构模型；
[0063]
图4是通过本发明建立的混凝土材料细观格构模型进行求解过程流程图；
[0064]
图5a-图5d是通过本发明的模拟方法进行混凝土材料损伤演化过程。
具体实施方式
[0065]
现在结合附图对本发明作进一步详细的说明。本技术的描述中，需要理解的是，术语“左侧”、“右侧”、“上部”、“下部”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，“第一”、“第二”等并不表示零部件的重要程度，因此不能理解为对本发明的限制。本实施例中采用的具体尺寸只是为了举例说明技术方案，并不限制本发明的保护范围。
[0066]
如背景技术中阐述的，目前针对混凝土材料损伤演化过程所建立的数值模型具有较多的自由度，导致整个计算量非常大，不仅准确率降低，而且计算效率也偏低。因此亟需研究出一种新的模拟方法，在提高计算速率以及准确率的同时，充分了解混凝土材料整个损伤演化过程。
[0067]
图1所示是本技术提供的基于格构模型的混凝土材料损伤演变过程的模拟方法，首先依据构件(这里即为待试验的混凝土材料)尺寸生成单元节点，通过单元节点生成单元信息矩阵，完成数值模型的搭建后设置材料单元以及材料参数，生成赋予材料特性的混凝土材料细观格构模型；接着，基于timoshenko梁理论建立混凝土材料细观格构模型的有限元方程，包括计算单元刚度矩阵以及组装总刚矩阵；由于混凝土材料的特殊性能，需要采用mazars损伤本构模型模拟混凝土材料力学性能的劣化；接着按位移加载的牛顿迭代法对构建的混凝土材料细观格构模型进行求解，计算过程大概包括位移求解、计算损伤变量以及计算单元刚度削减方程，最后对计算结果进行损伤分析，得到混凝土材料损伤演化过程的仿真模拟；这里需要阐述的是，位移求解、损伤变量的计算以及单元刚度削减方程的计算，均是因为混凝土这种材料本身的特殊性，在将混凝土材料划分为梁单元后，考虑剪切变形的影响，分析梁单元力学性能的劣化，以此使得建立的混凝土材料细观格构模型的有限元方程在进行求解时更加全面以及准确。
[0068]
具体的，关于构建混凝土材料细观格构模型，包括步骤s11：确定混凝土材料的尺寸，在待构建的混凝土材料细观格构模型边界上生成固定单元节点，在混凝土材料细观格
构模型区域内随机生成单元节点，遍历所有单元节点，并按照由下至上、由左至右的顺序，将单元节点全部存储于节点信息矩阵内；
[0069]
步骤s12：继续遍历所有单元节点，并按照由下至上、由左至右的顺序将所有单元节点依次相连，形成连接后的单元节点编号并存放在单元信息矩阵内，完成格构模型的几何搭建；
[0070]
步骤s13：确定各级配骨料的粒径范围及个数，在混凝土材料细观格构模型区域内生成各骨料级配的随机骨料，沿着骨料边界向外等厚度拓展界面过渡层，将生成的骨料信息存入骨料信息矩阵中；
[0071]
步骤s14：结合步骤s12中搭建的格构模型，将单元中点的位置对步骤s13中的骨料进行识别，对于不同相材料的单元赋予相应的材料参数，即若单元中点的位置位于骨料区域内，则将此单元判定为骨料单元，若单元中点的位置位于骨料边界与界面过渡层之间，则将此单元判定为界面过渡层，若单元中点的位置位于界面过渡层区域外，则将此单元判定为砂浆单元。
[0072]
基于timoshenko梁理论建立混凝土材料细观格构模型的有限元方程包括步骤s21：构建的混凝土材料细观格构模型选用timoshenko梁单元，并根据有限元理论推导出剪切刚度矩阵和弯曲刚度矩阵；
[0073]
步骤s22：剪切刚度矩阵采用缩减积分代替精确积分，将剪切刚度矩阵和弯曲刚度矩阵叠加，获取单元刚度矩阵；
[0074]
步骤s23：通过坐标转换矩阵将单元坐标系下的单元刚度矩阵转化为整体坐标系下的单元刚度矩阵，再将各个整体坐标系下的单元刚度矩阵中的元素按照编码规则“对号入座”叠加至总体刚度矩阵中，得到混凝土材料细观格构模型的有限元方程，其中总体刚度矩阵为变量，其与混凝土材料细观格构模型中单元的损伤程度有关联；
[0075]
混凝土材料细观格构模型的有限元方程为
[0076]k·
q＝p
ꢀꢀꢀ
(1)
[0077]
公式(1)中，k为总体刚度矩阵，q为位移向量，p为荷载向量。
[0078]
采用mazars损伤本构模型模拟混凝土材料力学性能的劣化包括：步骤s31：得到的混凝土材料细观格构模型有限元方程中，总体刚度矩阵与混凝土材料细观格构模型中单元的损伤程度有关联，采用mazars损伤本构模型模拟混凝土材料力学性能的劣化，通过损伤变量d来表示混凝土材料细观格构模型中单元的损伤程度；
[0079]
当d＝0，表示混凝土材料细观格构模型中单元无损伤；
[0080]
当单元等效应变大于应变阈值时，混凝土材料损伤开始演化；
[0081]
当d＝1，表示混凝土材料细观格构模型中单元失效破坏；
[0082]
其中，损伤变量d的表达公式为
[0083][0084]
公式(2)中，ε
p
为应变阈值，ε为单元等效应变，a
t
和b
t
为混凝土材料常数，对于timoshenko梁单元，其单元等效应变取为单元的最大主拉应变；
[0085]
步骤s32：基于lemaitre应变等效原理，得到混凝土材料细观格构模型中单元的刚度削减方程，通过刚度的不断削减体现混凝土材料整体力学性能的劣化；
[0086]
混凝土材料细观格构模型中单元的刚度削减方程为
[0087][0088]
公式(3)中，e为单元无损状态下的弹性模量，为单元受损之后的弹性模量，从公式(3)中明显可知格构单元的弹性模量随着单元损伤值会发生变化。
[0089]
图4所示，按位移加载的牛顿迭代法对混凝土材料细观格构模型进行计算求解，将计算结果经过处理对混凝土材料的损伤演化过程进行仿真分析，包括步骤s41：对混凝土材料细观格构模型施加边界条件，在加载点处施加位移荷载并将其分为k个荷载步；
[0090]
步骤s42：在第k个荷载步内，给定初始位移增量dk，可计算出各格构单元应变εk[0091]
εk＝b
·dk
ꢀꢀꢀ
(4)
[0092]
公式(4)中，b为单元应变矩阵；将位移增量与单元应变分别存储于位移增量列阵与单元应变矩阵内；
[0093]
步骤s43：在第k个荷载步内，将各格构的单元等效应变与应变阈值进行比对，若单元等效应变大于应变阈值ε
p
，则代入公式(2)，得到单元发生损伤并更新损伤变量dk[0094][0095]
此时在第k个荷载步内，单元应力为
[0096]
σk＝(1-dk)
·d·
εkꢀꢀꢀ
(5)
[0097]
公式(5)中，d为单元无损状态下的单元弹性矩阵，εk为格构单元应变，且εk＝b
·dk
；将更新后的损伤变量dk与单元应力σk分别存储于损伤变量列阵与单元应力矩阵中；上述即是根据单元的应变结合削减后弹性模量更新单元节点力；
[0098]
步骤s44：根据第k个荷载步内的单元节点内力与外力是否平衡来判断整个混凝土材料是否处于平衡状态，若结构处于平衡状态，则该荷载步计算收敛，更新所有单元的历史变量并进入第k+1个荷载步的迭代计算过程；
[0099]
此时由内应力产生的节点力必须与所施加的节点外力大小相等且方向相反，整个结构平衡的弱形式为
[0100][0101]
公式(6)中，为满足边界条件的任意虚位移，tk为第k个荷载步中的面力，fk为第k个荷载步中的体力；
[0102]
若整体结构未处于平衡状态，则节点内力与节点外力之间的残余力向量rk表示为
[0103][0104]
公式(7)中，n为单元形函数；
[0105]
步骤s45：在第k个荷载步中进行多次的循环迭代，当残余力向量小于预设的容差值时，则认为该荷载步内的迭代计算收敛并进入第k+1个荷载步的迭代计算过程，其中每一步的迭代公式为
[0106][0107]
公式(8)中，为第k个荷载步内结构的切线刚度矩阵，为该荷载步内的位移增量，i表示该荷载步中的第i次迭代；
[0108]
步骤s46：将所有荷载步的计算结果存储在结果信息矩阵中，计算完成后对结果信息矩阵内的元素进行后处理，模拟出混凝土材料的损伤演变过程。
[0109]
接下来本技术基于上述的模拟方法提供一个实施例，如图2所示是实施例提供的混凝土细观格构模型，图中圆形结构外部的区域为砂浆，圆形结构外围的白色环状部分为界面过渡区，圆形结构内部的区域为骨料，该模型尺寸为50mm
×
50mm，模型包含了10400个节点与30797个单元；本实施例采用三级配的随机圆形骨料投放方式，其中大粒径骨料6个，粒径范围为8-10mm；中等粒径骨料9个，粒径范围为6-8mm；小粒径骨料15个，粒径范围为2-6mm，界面过渡区(itz)厚度取0.5mm；本实施例中骨料单元的弹性模量取45000mpa，假设骨料不发生损伤，砂浆单元的弹性模量选用25000，抗拉强度为5.00mpa，itz单元的弹性模量取为18750mpa，抗拉强度为2.81mpa，三相材料的泊松比均取为0.2。
[0110]
在建立timoshenko梁格构模型的有限元方程时，格构模型的节点数量为10400，而每个节点均有3个自由度，因此所集成的整体刚度矩阵的阶数为31200
×
31200。
[0111]
接下来通过格构单元刚度的不断削减来体现混凝土试件整体力学性能的劣化，在本实施例中，假设骨料单元不发生破坏，因此对其不设置应变阈值，将砂浆单元的应变阈值设为200e-6，itz单元的应变阈值设为150e-6，三相材料的a
t
均设置为0.84，b
t
均设置为22900。
[0112]
遍历模型左边缘的所有节点，约束其竖直方向与水平方向的位移，在模型右边缘的所有节点上施加水平方向的位移荷载，位移加载总量为0.02mm；在格构模型计算过程中需将位移加载总量分成多个位移荷载步。在每一荷载步内，首先根据结构的初始状态计算出相应的节点外荷载增量，并同时由位移增量计算出每一个单元的应变，进而得出各单元的等效应变值。图3所示是mazars损伤本构模型，图示中，(ε
p
，σ
p
)即为应变阈值，若单元等效应变值超过了损伤本构所设置的应变阈值，则该单元产生了损伤。各格构单元的弹性模量将随着单元损伤值发生变化，从而导致整个结构的刚度发生劣化。根据单元的应变与削减后的弹性模量更新单元节点力，进而得到整个结构的节点内力。若结构的节点内力与节点外荷载的差足够小，则所有节点处于近似的平衡状态，此时认为计算收敛并进入下一荷载步；在每一荷载步的数值计算过程中，需要将所有节点的位移、荷载反力与各单元的应变、应力以及损伤值记录在信息矩阵中，以便后续荷载步或后处理过程中提取利用；当所有荷载步计算完毕后，经后处理得到的混凝土损伤演化过程如图5所示，图中黑色区域代表损伤值为0，白色区域表示损伤值为1，由图5a至图5d可以看出，白色区域逐渐扩大，因此模型内部单元的损伤值从0-1逐渐演化。
[0113]
由此可知，本技术提供的模拟方法，可以对混凝土材料的失效破坏激励进行全过程的损伤演化，充分考虑了剪切变形的影响，并且将混凝土材料划分为梁单元，将一个复杂的多轴问题转化为多个简单的单轴问题，保证模拟有效性的同时大大提升了计算效率。
[0114]
本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本技术所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。
[0115]
本技术中所述的“和/或”的含义指的是各自单独存在或两者同时存在的情况均包括在内。
[0116]
本技术中所述的“连接”的含义可以是部件之间的直接连接也可以是部件间通过其它部件的间接连接。
[0117]
以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项发明技术思想的范围内，进行多样的变更以及修改。本项发明的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。