CN117592338B - 一种考虑材料非线性的应变beso混凝土构件优化方法 - Google Patents
一种考虑材料非线性的应变beso混凝土构件优化方法 Download PDFInfo
- Publication number
- CN117592338B CN117592338B CN202311661221.9A CN202311661221A CN117592338B CN 117592338 B CN117592338 B CN 117592338B CN 202311661221 A CN202311661221 A CN 202311661221A CN 117592338 B CN117592338 B CN 117592338B
- Authority
- CN
- China
- Prior art keywords
- strain
- sensitivity
- optimization
- unit
- matrix
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000005457 optimization Methods 0.000 title claims abstract description 82
- 238000000034 method Methods 0.000 title claims abstract description 50
- 239000004567 concrete Substances 0.000 title claims abstract description 38
- 239000000463 material Substances 0.000 title claims abstract description 27
- 230000035945 sensitivity Effects 0.000 claims abstract description 72
- 238000013461 design Methods 0.000 claims abstract description 56
- 238000001914 filtration Methods 0.000 claims abstract description 19
- 238000004458 analytical method Methods 0.000 claims abstract description 8
- 238000012935 Averaging Methods 0.000 claims abstract description 6
- 238000013138 pruning Methods 0.000 claims abstract description 3
- 239000011159 matrix material Substances 0.000 claims description 47
- 238000006073 displacement reaction Methods 0.000 claims description 19
- 238000004364 calculation method Methods 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 11
- 239000011150 reinforced concrete Substances 0.000 claims description 11
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 6
- 238000011084 recovery Methods 0.000 claims description 6
- 230000002787 reinforcement Effects 0.000 claims description 6
- 230000007246 mechanism Effects 0.000 claims description 4
- 238000005354 coacervation Methods 0.000 claims description 3
- 239000000109 continuous material Substances 0.000 claims description 3
- 238000012217 deletion Methods 0.000 claims description 3
- 230000037430 deletion Effects 0.000 claims description 3
- 238000013178 mathematical model Methods 0.000 claims description 3
- 239000011343 solid material Substances 0.000 claims description 3
- 238000012916 structural analysis Methods 0.000 claims description 3
- 230000007704 transition Effects 0.000 claims description 3
- 238000012163 sequencing technique Methods 0.000 claims description 2
- 230000003014 reinforcing effect Effects 0.000 claims 1
- 230000006872 improvement Effects 0.000 description 9
- 238000004422 calculation algorithm Methods 0.000 description 5
- 238000005336 cracking Methods 0.000 description 3
- 238000003079 width control Methods 0.000 description 3
- 229910000831 Steel Inorganic materials 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 239000011148 porous material Substances 0.000 description 2
- 239000010959 steel Substances 0.000 description 2
- 229910001294 Reinforcing steel Inorganic materials 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000003592 biomimetic effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000013013 elastic material Substances 0.000 description 1
- 230000009395 genetic defect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 108090000623 proteins and genes Proteins 0.000 description 1
- 230000002040 relaxant effect Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
- Complex Calculations (AREA)
Abstract
本发明涉及一种考虑材料非线性的应变BESO混凝土构件优化方法，其步骤如下：S1：确定设计域；S2：离散设计域；S3：预设优化参数；S4：计算迭代次数；S5：计算目标体积率；S6：有限元分析；S7：求解不平衡力；S8：判断是否收敛；S9：将荷载步一分为二，重新加载，返回执行S6；S10：计算灵敏度；S11：过滤灵敏度；S12：与历史灵敏度平均处理；S13：进行第一次设计变量更新；S14：过滤设计变量；S15：再执行一次S13操作；S16：删减和增添单元；S17：判断剩余体积是否小于目标体积；S18：返回执行S6‑S16；S19：寻找最优拓扑结构；S20：输出。本发明能显著降低结构应变峰值。
Description
技术领域
本发明涉及钢筋混凝土结构拓扑优化领域，尤其涉及一种考虑材料非线性的应变BESO混凝土构件优化方法。
背景技术
世界各国对混凝土结构的现行规范中，都包含了对构件抗裂性或最大裂缝宽度限值的相应要求。那么，在设计阶段对裂缝宽度的预测成为了一项重要的工作，陶慕轩等将裂缝带理论拓宽到配筋混凝土构件，并据此提出基于弥散裂缝模型和开裂应变近似计算结构裂缝宽度的方法。
拓扑优化是一种在给定设计域内，针对指定目标、约束和边界条件，寻找优化的结构拓扑或材料布局的方法。自从它被扩展到须考虑几何非线性、材料非线性，甚至须同时考虑以上两者的混凝土结构设计，除了计算效率可能成为应用的控制因素之一外，潜在的大位移可能会导致低密度单元中的切线刚度矩阵变得不确定甚至负定。因此，需要通过去除低密度元素或放松收敛准则等方法来避免这个问题，这正是在众多拓扑优化方法中，具有离散特性的Evolutionary Structural Optimization-type(ESO-type)方法所擅长完成的工作，而其中最为成熟、应用也最为广泛的就是Bi-directional ESO(BESO)方法。
BESO方法在混凝土结构设计应用中多基于应力开展。Duysinx等认为这样的方式可以有效降低结构中的应力集中现象；Rozvany等提出了一种最优性准则方法continuum-type optimality iteria(COC)来获得降低结构峰值应力的最优拓扑；Harzheim等结合仿生生物生长的算法和最优性准则方法来调整结构的几何形状和荷载分布以降低结构的应力水平。基于这样的优化，通常可以设计出应力水平更低或应力分布更均匀的混凝土结构与构件。
然而，基于应力的拓扑优化是难以实现以裂缝宽度控制为目标的混凝土结构设计的。因为混凝土是一种初始即包含微小缺陷和不均匀性、存在应变软化的材料。首先，基于应力的拓扑优化方法主要建立在线弹性材料力学假设之上，无法准确模拟混凝土材料实际存在的软化行为，继而难以估计这种软化行为的影响；其次，根据第二强度理论——最大伸长线应变理论，混凝土开裂由应变控制，而不是由应力的控制，而且应变软化的存在，会导致应力峰值与应变峰值通常发生在不同的时刻。所以，基于非线性应变开展拓扑优化，对于将拓扑优化引入到混凝土结构或构件裂缝宽度控制中，是一项不可或缺的工作。
发明内容
本发明的目的是提供一种考虑材料非线性的应变BESO混凝土构件优化方法，本方法能显著降低结构应变峰值，从而达到控制结构最大裂缝宽度的目标。
为达到上述目的而采用了一种考虑材料非线性的应变BESO混凝土构件优化方法，其步骤如下：
S1：确定设计域，荷载工况和边界条件；
S2：离散设计域；
S3：预设优化参数；
S4：计算优化所需总迭代次数；
S5：计算当代目标体积率；
S6：进行有限元分析；
S7：求解不平衡力；
S8：判断不平衡力是否满足收敛条件，若否，则执行S9，若是，则执行S10；
S9：将荷载步一分为二，并重新加载，返回执行S6；
S10：计算灵敏度；
S11：过滤灵敏度；
S12：与历史灵敏度平均处理；
S13：进行第一次设计变量更新；
S14：过滤设计变量；
S15：再执行一次S13操作，接着执行S16操作；
S16：删减和增添单元；
S17：判断剩余体积是否小于目标体积，若否，则返回执行S5；若是，则执行S18；
S18：返回执行S6-S16，直至最小化混凝土构件中的最大裂缝宽度；
S19：寻找并确定最优拓扑结构；
S20：输出最优拓扑结构。
作为本发明的进一步改进，基于最大裂缝宽度最小化目标等价于结构峰值应变最小化目标，建立优化的数学模型转化为如下：
寻找：x＝[x1,x2,…,xi,…,xn]
最小化：
服从：R(U,x)＝0
V(x)＝∑xivi＝V*
xi＝0or 1,i＝1,...,n,
式中为结构中Mises应变峰值，εeq,i为第i个单元中心点处的Mises应变；V、vi、V*分别为结构材料总体积、第i个单元的体积和结构目标体积；x为设计变量向量，即结构拓扑的解空间；xi为第i个单元的设计变量值，可取0或1，分别代表已被删除和存留状态；R为节点不平衡力向量，按下式计算：
式中U是位移矩阵；B是几何函数矩阵，用于在已知位移的情况下计算应变；以平面问题为例，假定材料本构与设计变量xi无关，Mises应变εeq按下式计算：
式中εx、εy分别为x、y方向的应变，γxy为剪切应变；ε为单元应变向量，v是材料的泊松比；A是应变系数矩阵，按下式计算：
式中
引入p范数函数近似计算结构峰值应变，将式(1)中的目标函数转化为：
最小化：
式中εPN为应变的凝聚函数；p为人为设定的范数值，当p＝1时，εPN为所有单元应变之和。
作为本发明的进一步改进，S10中，优化灵敏度的步骤如下：
使用伴随方法推导灵敏度，忽略非线性分析状态转换引起的不可微性；
首先，采用线性插值模型建立设计变量与单元材料属性的关系：
Di＝xiD0(6)
式中Di表示第i号单元参与有限元计算时的实际弹性矩阵，D0表示实体材料的弹性矩阵；
再采用连续性材料Mises应力插值模型的Mises应力表达式：
σi＝DiBiui/xi＝D0Biui(7)
式中σi和ui分别为第i号单元的应力向量与节点位移向量；
然后，引入与荷载子步数nload相同数量的拉格朗日乘子λk，每一个λk与未知的位移矩阵U具有相同的维数，并且使任意荷载增量步，第k步与其上一步，第k-1步的不平衡力R的差值ΔRk为0：
式中Li是索引矩阵，用于从整体矩阵中索引局部矩阵，第i号单元的节点位移向量ui＝LiU，第i号单元的节点力向量fi＝LiF，于是，式(5)中的原目标函数等价于：
显然，式(5)与式(9)的这种等价关系，对于单元灵敏度同样成立，即满足而且对于按力加载，/>所以，将式(9)等式两边对xj求偏导，可得：
式中是非线性系统在第k个荷载增量步中达到平衡时的结构全局有限元切线刚度矩阵；
在第k个荷载增量步的结构分析中，近似有：
式中ΔFk和ΔUk分别为第k个荷载增量步时的荷载增量矩阵和结构位移增量矩阵；
将式(11)两边同时对设计变量xj求偏导，可得：
于是，根据式(10)可转化得到：
对式(10)和式(13)取平均，可得：
取此时/>的系数为0，式(14)变换成：
式中为第i号单元在第k个荷载增量步时参与有限元计算的单元切线刚度矩阵，在式(15)上再增乘了一个反号因子“-xi”，得到每个单元的原始灵敏度表达式：
作为本发明的进一步改进，S11包括：
对单元灵敏度进行如下过滤操作：
式中为线性权重因子，计算如下：
式中rsen为灵敏度过滤半径；Δ(i,j)为第i号单元和第j号单元中心点之间的距离。
作为本发明的进一步改进，S12包括：
将过滤后的灵敏度与前两代的历史灵敏度信息进行平均处理以进一步过滤灵敏度：
式中(l)代表优化迭代循环次数。
作为本发明的进一步改进，S13中设计变量更新机制的步骤如下：
得到每个单元灵敏度数值后，按下式执行第1次设计变量更新操作：
式中和/>分别为灵敏度上阈值和下阈值，首先，据下式计算当代的目标体积V(l)：
V(l)＝max{V*,(1-cer)V(l-1)} (21)
式中cer表示进化率，根据V(l)计算当代优化执行完删除和恢复操作后将减少的单元总数目n0，然后，对全部存留单元按灵敏度值从小到大进行排序，将第n0位次上单元的灵敏度值作为阈值设置参考值再据事先预设的单元恢复数目最大比例/>以目标体积V(l)包含单元总数目为基数，计算当代优化中可执行恢复操作的最大单元数目n1，对全部已被删除单元按灵敏度值从大到小进行排序，取第n1位次上单元的灵敏度值/>与/>进行比较，以较小者作为/>的取值，若/>等于/>则/>也等于/>若/>等于/>则按n0+n1计算出当代需执行删除操作的单元数目n2，再取之前按灵敏度值从小到大的排序的全部存留单元中第n2位次上单元的灵敏度值，作为/>的取值。
作为本发明的进一步改进，S14包括：
增加对设计变量的过滤操作，其过滤方程如下：
式中线性权重因子的大小计算：
式中rden为设计变量过滤半径，取值需大于灵敏度过滤半径rsen，以保证优化过程的收敛性和网格独立性；
在执行了以上过滤操作后，设计变量不再只有0和1两种单元状态，而是包含0到1之间的任意值。
作为本发明的进一步改进，S17还包括：在优化至体积约束后，需要继续迭代循环直至满足下式时终止优化：
式中τ和Nitr分别为容许收敛误差和预设差值总数，m为当前代数。
作为本发明的进一步改进，S19包括：
最后，从优化至体积约束后迭代循环所得的拓扑结构中，按目标函数进行评价，取其中最小的结构方案作为最优拓扑，找出最优拓扑结构。
作为本发明的进一步改进，构建钢筋混凝土单元刚度矩阵：
式中DC为混凝土的弹性矩阵，ρx为水平方向的配筋率，ρy为竖直方向的配筋率，υ为混凝土的泊松比；Ec和Es分别为混凝土和钢筋的杨氏模量
按下式定义裂缝应变作为描述裂缝宽度的等效指标：
式中εtol为带筋混凝土的总应变，εe为其中的弹性应变。
本发明运用BESO方法，引入p范数近似估计应变以建立包含材料非线性的结构优化模型；随后利用伴随法推导了相应的灵敏度，开发考虑材料非线性的结构应变拓扑优化方法；再探讨高度非线性应变下的设计变量更新机制；最后通过数值算例开展新方法有效性的验证工作。
附图说明
图1为基于BESO算法的应变优化流程图。
图2为开孔简支深梁的几何模型和边界条件。
图3为简支梁算例的拓扑结构随代数的演化过程(p＝6)及它们的应变分布。
图4为开洞简支梁算例的最优拓扑结构及其应变基本情况。
图5为钢筋混凝土悬臂梁的优化设计域、边界条件和荷载工况。
图6为应变优化方法和应力优化方法对比情况。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制；术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性；此外，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
实施例1
1.基于BESO的混凝土构件应变优化
1.1优化模型
混凝土裂缝通常发生在构件中应变最大的部位，所以最大裂缝宽度最小化目标，等价于结构峰值应变最小化目标，于是，优化的数学模型转化为如下：
寻找：x＝[x1,x2,…,xi,…,xn]
最小化：
服从：R(U,x)＝0 (1)
V(x)＝∑xivi＝V*
xi＝0or 1,i＝1,...,n,
(式中为结构中Mises应变峰值，εeq,i为第i个单元中心点处的Mises应变；V、vi、V*分别为结构材料总体积、第i个单元的体积和结构目标体积；x为设计变量向量，即结构拓扑的解空间；xi为第i个单元的设计变量值，可取0或1，分别代表已被删除和存留状态；R为节点不平衡力向量，按下式计算：
式中U是位移矩阵；B是几何函数矩阵，用于在已知位移的情况下计算应变。以平面问题为例，假定材料本构与设计变量xi无关，Mises应变εeq按下式计算：
由于混凝土开裂由其应变控制，而对于非一维杆系结构，单元应力与应变分布规律不尽相同，所以应力优化存在无法解决的基因缺陷。因此，本实施例提出一种考虑材料非线性、以结构峰值应变最小化为目标函数的应变拓扑优化方法。该法引入了p范数近似估计非线性应变，以伴随法推导灵敏度，能显著降低结构应变峰值，从而实现控制结构最大裂缝宽度目标下的混凝土结构优化设计。此外，在保证数值计算稳定性的前提下，p取值越大，最优拓扑结构越优。
式中εx、εy分别为x、y方向的应变，γxy为剪切应变；ε为单元应变向量，v是材料的泊松比；A是应变系数矩阵，按下式计算：
式中
由于应变具有和应力相似的局部性质，无法通过显式的数学表达式描述，因此，引入p范数函数近似计算结构峰值应变，将式(1)中的目标函数转化为：
最小化：
式中εPN为应变的凝聚函数；p为人为设定的范数值。当p＝1时，εPN为所有单元应变之和，p值越大，εPN越接近于但同时也可能导致优化结果的发散。因此，在实际操作中，p值的选取不仅需要考虑εPN和/>的近似程度，还需要权衡优化过程的稳定性。
1.2优化灵敏度
本实施例使用伴随方法推导灵敏度，忽略非线性分析状态转换引起的不可微性。首先，采用线性插值模型建立设计变量与单元材料属性的关系：
Di＝xiD0 (6)
式中Di表示第i号单元参与有限元计算时的实际弹性矩阵，D0表示实体材料的弹性矩阵。
再采用连续性材料Mises应力插值模型的Mises应力表达式：
σi＝DiBiui/xi＝D0Biui (7)
式中σi和ui分别为第i号单元的应力向量与节点位移向量。
然后，引入与荷载子步数nload相同数量的拉格朗日乘子λk，每一个λk与未知的位移矩阵U具有相同的维数，并且使任意荷载增量步(第k步)与其上一步(第k-1步)的不平衡力R的差值ΔRk为0：
式中Li是索引矩阵，用于从整体矩阵中索引局部矩阵，如：第i号单元的节点位移向量ui＝LiU，第i号单元的节点力向量fi＝LiF。于是，式(5)中的原目标函数等价于：
显然，式(5)与式(9)的这种等价关系，对于单元灵敏度同样成立，即满足 而且对于按力加载，/>所以，将式(9)等式两边对xj求偏导，可得：
式中是非线性系统在第k个荷载增量步中达到平衡时的结构全局有限元切线刚度矩阵。
在第k个荷载增量步的结构分析中，近似有：
式中ΔFk和ΔUk分别为第k个荷载增量步时的荷载增量矩阵和结构位移增量矩阵。
将式(11)两边同时对设计变量xj求偏导，可得：
于是，根据式(10)可转化得到：
对式(10)和式(13)取平均，可得：
取此时/>的系数为0，式(14)变换成：
式中为第i号单元在第k个荷载增量步时参与有限元计算的单元切线刚度矩阵，/>为第i号单元在第k个荷载增量步时单元切线刚度矩阵。为了与传统BESO算法的灵敏度保持相同的增减单调性，在式(15)上再增乘了一个反号因子“-xi”，得到每个单元的原始灵敏度表达式：
为了避免拓扑优化中的网格依赖性问题和棋盘格现象等常见数值问题，对单元灵敏度进行如下过滤操作：
式中为线性权重因子，计算如下：
式中rsen为灵敏度过滤半径，其取值不宜过大，以确保已被删除的单元存在被恢复的可能。通常取值为单元边长的2～3倍；Δ(i,j)为第i号单元和第j号单元中心点之间的距离。
BESO方法中设计变量的离散本质，以及应变问题较强的非线性特征，都易导致优化过程出现较大的波动，因此，本实施例还增加了将过滤后的灵敏度与前两代的历史灵敏度信息进行平均处理的操作以进一步过滤灵敏度：
式中(l)代表优化迭代循环次数。
1.3设计变量更新机制
得到每个单元灵敏度数值后，按下式执行第1次设计变量更新操作：
式中和/>分别为灵敏度上阈值和下阈值，首先，据下式计算当代的目标体积V(l)：/>
V(l)＝max{V*,(1-cer)V(l-1)} (21)
式中cer表示进化率，一般预设为1％～2％。再据V(l)计算当代优化执行完删除和恢复操作后将减少的单元总数目n0。然后，对全部存留单元按灵敏度值从小到大进行排序，将第n0位次上单元的灵敏度值作为阈值设置参考值再据事先预设的单元恢复数目最大比例/>(以目标体积V(l)包含单元总数目为基数)，计算当代优化中可执行恢复操作的最大单元数目n1，对全部已被删除单元按灵敏度值从大到小进行排序，取第n1位次上单元的灵敏度值/>与/>进行比较，以较小者作为/>的取值。若/>等于/>则/>也等于/>若等于/>则按n0+n1计算出当代需执行删除操作的单元数目n2，再取之前按灵敏度值从小到大的排序的全部存留单元中第n2位次上单元的灵敏度值，作为/>的取值。
为提高算法的稳定性和鲁棒性，本实施例在传统BESO算法仅对灵敏度执行过滤操作的基础上，增加了对设计变量的过滤操作。与灵敏度过滤类似，其过滤方程如下：
式中线性权重因子的大小计算：
式中rden为设计变量过滤半径，取值需大于灵敏度过滤半径rsen，以保证优化过程的收敛性和网格独立性。
在执行了以上过滤操作后，设计变量不再只有0和1两种单元状态，而是可能包含0到1之间的任意值。为了保证拓扑结构是一个明确的离散结构，需按式(20)再执行一次设计变量更新操作。
1.4优化流程
由于应变的高度非线性特征，传统的基于特定迭代次数的BESO方法启发式收敛准则通常不适用于基于应变的优化。因此，在优化至体积约束后，需要继续迭代循环直至满足下式时终止优化：
式中τ和Nitr分别为容许收敛误差和预设差值总数，在本实施例中分别设置为0.001和5；m为当前代数。最后，从优化至体积约束后迭代循环所得的拓扑结构中，按目标函数进行评价，找出最优拓扑结构。
整个优化流程如图1所示：
2.数值算例
2.1基本参数
实施例算例均基于matlab平台自主编程开展非线性有限元分析，采用钢筋混凝土整体模型，即将钢筋均匀弥散于混凝土中，视钢筋混凝土为一种复合材料，按下式构建单元刚度矩阵：
式中DC为混凝土的弹性矩阵，水平方向的配筋率ρx＝0.2％，竖直方向的配筋率ρy＝0.15％。混凝土的泊松比υ取0.3；考虑到算例中结构的整体应力水平不太高，混凝土受压和钢筋受拉均采用线弹性本构，杨氏模量Ec和Es分别取3.15×104N/mm2和2×105N/mm2；混凝土受拉采用两折线本构，线弹性段的杨氏模量Ec取3.15×104N/mm2，直到σs0＝2.2N/mm2时开裂，然后进入下降段，直到达到预设的断裂能Gf＝0.124N/mm2时失效。
优化中，均分为10个荷载步进行加载，cer和分别设为1％和0.5％，灵敏度和设计变量的过滤半径分别取3倍和5倍单元边长。考虑到应变问题的非线性特征，在优化至目标体积后，再根据式(24)进行至多20次的额外迭代，取其中/>最小的结构方案作为最优拓扑。
此外，由于本实施例优化的原始目标是最小化混凝土构件中的最大裂缝宽度，而优化基于的仿真采用弥散裂缝模型，无法直接读取裂缝宽度。因此，基于经典裂缝带理论，按下式定义裂缝应变作为描述裂缝宽度的等效指标：
式中εtol为带筋混凝土的总应变，εe为其中的弹性应变。
2.2开孔简支深梁
某钢筋混凝土开孔简支深梁的优化设计域、边界条件和荷载工况如图2所示，采用4节点平面应力平行四边形单元和边长为10mm的正方形网格进行有限元网格划分，结构目标体积率设为50％。为避免集中荷载作用点的应力集中，在加载点和两支座下均设置宽度为100mm的刚性垫块，约束垫块，使其不能与构件发生相对位移。
本实施例先以p＝6下开展的开孔简支深梁算例优化为例，展示优化过程。采用本实施例应变优化方法的拓扑结构随代数的演化过程及它们的应变分布如图3所示，从图3可知，随着优化整体上呈现下降的趋势，在第37次迭代后，应变水平基本维持在9.25×10-5左右，由此表明了优化目标的实现。对于该曲线上出现较大波动的原因，主要是由于结构删除单元造成了暂时的薄弱部位，一般可以通过提高复活率或者调整进化率来稳定优化过程。
不同p值下开孔简支深梁算例应变优化的最优拓扑结构及它们的应变基本情况列于图4中。由图4可知，随着p值的增大，对于演化出的最优拓扑结构，其包含的杆件愈发地规整和清晰，应变分布也越来越均匀，表明其能更准确地反映传力路径；裂缝应变与Mises应变峰值/>有着相同的变化规律，说明采用p范数函数近似原目标函数是一定程度上等效且合理的，可以实现控制结构应变峰值的优化目标，从而达到控制结构裂缝宽度的原始目的。但需要特别说明的是，p取值过大(如以上简支梁算例中，当p>6后)，可能造成拓扑结构较为混乱。
2.3悬臂深梁
某钢筋混凝土悬臂梁的优化设计域、边界条件和荷载工况如图5所示，结构目标体积率设为55％，其余未尽参数与上一小节的深梁算例相同。在相同的优化参数条件下，对该算例分别开展应变优化和应力优化，并将所得所有最优拓扑结构都基于相同的非线性本构进行有限元分析，得到它们的结果均列于图6中。
从图6可知，p取值对优化的影响与上一节的深梁算例所得结论基本一致；较之应力优化，采用应变优化可得到值较小的最优拓扑结构，表明了其在控制裂缝宽度发育上的优势。这是由于应力优化的灵敏度基于单元Mises应力建立，但由式(3)可知，由于对于由平面单元或实体单元组成的结构，结构中的单元Mises应力和Mises应变在分布规律上并不完全相同，而根据第二强度理论，混凝土开裂由应变控制，因此，应力优化方法存在基因缺陷，理论上不太可能得到最优解，同时也表明了在以混凝土裂缝宽度控制为目标时，应变优化的必要性。
3.结论
(1)考虑材料非线性的应变拓扑优化方法，以最小化结构应变峰值为目标，采用p范数函数近似估算结构应变峰值，再基于伴随法推导灵敏度以开展优化。该方法能显著降低结构应变峰值，从而达到控制结构最大裂缝宽度的目标。
(2)一般情况下，p的取值越大，应变优化得到的最优拓扑结构也会越符合优化目标，但是当p的取值过大时，可能引起数值计算不稳定，从而导致优化结果发散。因此，建议一般情况下p值取6或7。
(3)通过对灵敏度和设计变量的过滤和更新等操作，可以降低应变的高度非线性特性对优化过程的影响，提高优化过程的稳定性，演化出清晰的类杆件的最优拓扑结构。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干等同替代或明显变型，而且性能或用途相同，都应当视为属于本发明的保护范围之内。
Claims (1)
1.一种考虑材料非线性的应变BESO混凝土构件优化方法，其特征在于具体步骤如下：
S1：确定设计域，荷载工况和边界条件；
S2：离散设计域；
S3：预设优化参数；
S4：计算优化所需总迭代次数；
S5：计算当代目标体积率；
S6：进行有限元分析；
S7：求解不平衡力；
S8：判断不平衡力是否满足收敛条件，若否，则执行S9，若是，则执行S10；
S9：将荷载步一分为二，并重新加载，返回执行S6；
S10：计算灵敏度；
S11：过滤灵敏度；
S12：与历史灵敏度平均处理；
S13：进行第一次设计变量更新；
S14：过滤设计变量；
S15：再执行一次S13操作，接着执行S16操作；
S16：删减和增添单元；
S17：判断剩余体积是否小于目标体积，若否，则返回执行S5；若是，则执行S18；
S18：返回执行S6-S16，直至最小化混凝土构件中的最大裂缝宽度；
S19：寻找并确定最优拓扑结构；
S20：输出最优拓扑结构；
基于最大裂缝宽度最小化目标等价于结构峰值应变最小化目标，建立优化的数学模型转化为如下：
式中为结构中Mises应变峰值，εeq,i为第i个单元中心点处的Mises应变；V、vi、V*分别为结构材料总体积、第i个单元的体积和结构目标体积；x为设计变量向量，即结构拓扑的解空间；xi为第i个单元的设计变量值，取0或1，分别代表已被删除和存留状态；R为节点不平衡力向量，按下式计算：
式中U是位移矩阵；B是几何函数矩阵，用于在已知位移的情况下计算应变；假定材料本构与设计变量xi无关，Mises应变εeq按下式计算：
式中εx、εy分别为x、y方向的应变，γxy为剪切应变；ε为单元应变向量，v是材料的泊松比；A是应变系数矩阵，按下式计算：
式中
引入p范数函数近似计算结构峰值应变，将式(1)中的目标函数转化为：
最小化：
式中εPN为应变的凝聚函数；p为人为设定的范数值，当p＝1时，εPN为所有单元应变之和；
S10中，优化灵敏度的步骤如下：
使用伴随方法推导灵敏度，忽略非线性分析状态转换引起的不可微性；
首先，采用线性插值模型建立设计变量与单元材料属性的关系：
Di＝xiD0 (6)
式中Di表示第i号单元参与有限元计算时的实际弹性矩阵，D0表示实体材料的弹性矩阵；
再采用连续性材料Mises应力插值模型的Mises应力表达式：
σi＝DiBiui/xi＝D0Biui (7)
式中σi和ui分别为第i号单元的应力向量与节点位移向量；
然后，引入与荷载子步数nload相同数量的拉格朗日乘子λk，每一个λk与未知的位移矩阵U具有相同的维数，并且使任意荷载增量步，第k步与其上一步，第k-1步的不平衡力R的差值ΔRk为0：
式中Li是索引矩阵，用于从整体矩阵中索引局部矩阵，第i号单元的节点位移向量ui＝LiU，第i号单元的节点力向量fi＝LiF，于是，式(5)中的原目标函数等价于：
显然，式(5)与式(9)的这种等价关系，对于单元灵敏度同样成立，即满足 而且对于按力加载，/>所以，将式(9)等式两边对xj求偏导，可得：
式中是非线性系统在第k个荷载增量步中达到平衡时的结构全局有限元切线刚度矩阵；
在第k个荷载增量步的结构分析中，近似有：
式中ΔFk和ΔUk分别为第k个荷载增量步时的荷载增量矩阵和结构位移增量矩阵；
将式(11)两边同时对设计变量xj求偏导，可得：
于是，根据式(10)可转化得到：
对式(10)和式(13)取平均，可得：
取此时/>的系数为0，式(14)变换成：
式中为第i号单元在第k个荷载增量步时参与有限元计算的单元切线刚度矩阵，在式(15)上再增乘了一个反号因子“-xi”，得到每个单元的原始灵敏度表达式：
S11包括：
对单元灵敏度进行如下过滤操作：
式中为线性权重因子，计算如下：
式中rsen为灵敏度过滤半径；Δ(i,j)为第i号单元和第j号单元中心点之间的距离；
S12包括：
将过滤后的灵敏度与前两代的历史灵敏度信息进行平均处理以进一步过滤灵敏度：
式中(l)代表优化迭代循环次数；
S13中设计变量更新机制的步骤如下：
得到每个单元灵敏度数值后，按下式执行第1次设计变量更新操作：
式中和/>分别为灵敏度上阈值和下阈值，首先，据下式计算当代的目标体积V(l)：
V(l)＝max{V*,(1-cer)V(l-1)} (21)
式中cer表示进化率，根据V(l)计算当代优化执行完删除和恢复操作后将减少的单元总数目n0，然后，对全部存留单元按灵敏度值从小到大进行排序，将第n0位次上单元的灵敏度值作为阈值设置参考值再据事先预设的单元恢复数目最大比例/>以目标体积V(l)包含单元总数目为基数，计算当代优化中可执行恢复操作的最大单元数目n1，对全部已被删除单元按灵敏度值从大到小进行排序，取第n1位次上单元的灵敏度值/>与/>进行比较，以较小者作为/>的取值，若/>等于/>则/>也等于/>若/>等于/>则按n0+n1计算出当代需执行删除操作的单元数目n2，再取之前按灵敏度值从小到大的排序的全部存留单元中第n2位次上单元的灵敏度值，作为/>的取值；
S14包括：
增加对设计变量的过滤操作，其过滤方程如下：
式中线性权重因子的大小计算：
式中rden为设计变量过滤半径，取值需大于灵敏度过滤半径rsen，以保证优化过程的收敛性和网格独立性；
在执行了以上过滤操作后，设计变量不再只有0和1两种单元状态，而是包含0到1之间的任意值；
S17还包括：在优化至体积约束后，需要继续迭代循环直至满足下式时终止优化：
式中τ和Nitr分别为容许收敛误差和预设差值总数，m为当前代数；
S19包括：
最后，从优化至体积约束后迭代循环所得的拓扑结构中，按目标函数进行评价，取其中最小的结构方案作为最优拓扑，找出最优拓扑结构；
构建钢筋混凝土单元刚度矩阵：
式中DC为混凝土的弹性矩阵，ρx为水平方向的配筋率，ρy为竖直方向的配筋率，υ为混凝土的泊松比；Ec和Es分别为混凝土和钢筋的杨氏模量
按下式定义裂缝应变作为描述裂缝宽度的等效指标：
式中εtol为带筋混凝土的总应变，εe为其中的弹性应变。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311661221.9A CN117592338B (zh) | 2023-12-06 | 2023-12-06 | 一种考虑材料非线性的应变beso混凝土构件优化方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311661221.9A CN117592338B (zh) | 2023-12-06 | 2023-12-06 | 一种考虑材料非线性的应变beso混凝土构件优化方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117592338A CN117592338A (zh) | 2024-02-23 |
| CN117592338B true CN117592338B (zh) | 2024-06-21 |
Family
ID=89921822
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311661221.9A Active CN117592338B (zh) | 2023-12-06 | 2023-12-06 | 一种考虑材料非线性的应变beso混凝土构件优化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117592338B (zh) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102323976A (zh) * | 2011-06-24 | 2012-01-18 | 武汉理工大学 | 混凝土桥梁收缩徐变及预应力损失计算方法 |
| CN110348102A (zh) * | 2019-07-04 | 2019-10-18 | 广州大学 | 基于反正切的动态进化率beso拓扑优化方法及其应用 |
Family Cites Families (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9194855B2 (en) * | 2014-02-28 | 2015-11-24 | Quipip, Llc | Systems, methods and apparatus for providing to a driver of a vehicle carrying a mixture real-time information relating to a characteristic of the mixture |
| CN106372347B (zh) * | 2016-09-08 | 2019-05-07 | 厦门大学嘉庚学院 | 改进双向渐进法的等效静载荷法动态响应拓扑优化方法 |
| JP7058902B2 (ja) * | 2019-04-17 | 2022-04-25 | 大連理工大学 | ハイブリッド繊維複合材料の板巻きシェル構造に対する高速協調最適化方法 |
| CN110851904A (zh) * | 2019-11-12 | 2020-02-28 | 武汉理工大学 | 一种外框内筒电视塔结构参数灵敏度快速分析方法及系统 |
| CN111523162B (zh) * | 2020-03-13 | 2021-04-13 | 湖南科技大学 | 基于两种极限状态的深梁分离式配筋多目标拓扑设计方法 |
| CN113094945B (zh) * | 2021-03-22 | 2023-04-11 | 中山大学 | 一种sa-beso联合拓扑优化方法 |
| CN115464745A (zh) * | 2022-09-28 | 2022-12-13 | 西安建筑科技大学 | 一种路径宽度可变的混凝土3d打印路径优化方法 |
| CN115640719A (zh) * | 2022-10-13 | 2023-01-24 | 湖南科技大学 | 基于拓扑优化构建复杂应力构件的静定桁架模型的方法 |
| CN115758832A (zh) * | 2022-11-23 | 2023-03-07 | 湖南科技大学 | 一种钢筋混凝土剪力墙的拟静力渐进演化拓扑优化方法 |
| CN115906572A (zh) * | 2022-11-29 | 2023-04-04 | 中铁第四勘察设计院集团有限公司 | 一种基于apdl的钢管混凝土轨枕结构拓扑优化设计方法 |
| CN115828702B (zh) * | 2022-12-21 | 2023-07-11 | 中交二公局第一工程有限公司 | 一种改进的gbeso算法及其工程优化设计应用方法 |
| CN117172066A (zh) * | 2023-09-13 | 2023-12-05 | 中交路桥建设有限公司 | 混凝土梁桥抗火性能form有限元可靠度分析方法 |
-
2023
- 2023-12-06 CN CN202311661221.9A patent/CN117592338B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102323976A (zh) * | 2011-06-24 | 2012-01-18 | 武汉理工大学 | 混凝土桥梁收缩徐变及预应力损失计算方法 |
| CN110348102A (zh) * | 2019-07-04 | 2019-10-18 | 广州大学 | 基于反正切的动态进化率beso拓扑优化方法及其应用 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117592338A (zh) | 2024-02-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Kočvara | On the modelling and solving of the truss design problem with global stability constraints | |
| CN110008512B (zh) | 一种考虑承载特性的负泊松比点阵结构拓扑优化方法 | |
| CN111523162B (zh) | 基于两种极限状态的深梁分离式配筋多目标拓扑设计方法 | |
| CN110069864B (zh) | 一种结合变密度法的改进双向渐进结构拓扑优化方法 | |
| CN111353246B (zh) | 一种混凝土构件设计的静动力多目标拓扑演化方法 | |
| CN114970366B (zh) | 一种功能梯度超材料结构优化设计方法及系统 | |
| CN109255142B (zh) | 基于小生境遗传算法的环形张拉整体结构拓扑优化方法 | |
| CN117592338B (zh) | 一种考虑材料非线性的应变beso混凝土构件优化方法 | |
| Ma et al. | Stress relaxation and sensitivity weight for bi-directional evolutionary structural optimization to improve the computational efficiency and stabilization on stress-based topology optimization | |
| CN111597724B (zh) | 考虑频带约束的结构动力学拓扑优化方法、系统 | |
| CN116702391B (zh) | 基于正则化的保型拓扑优化设计方法 | |
| Chan et al. | Nonlinear stiffness design optimization of tall reinforced concrete buildings under service loads | |
| CN115828702A (zh) | 一种改进的gbeso算法及其工程优化设计应用方法 | |
| Åldstedt et al. | Nonlinear time-dependent concrete-frame analysis | |
| CN110704912B (zh) | 一种应力约束下的桥架托臂结构拓扑优化方法 | |
| CN114329702B (zh) | 基于改进的差分进化算法标定设计反应谱的高稳定性方法和设备 | |
| CN113887098A (zh) | 一种基于疲劳应力鲁棒性和可靠性的金属结构拓扑优化方法 | |
| CN113361176B (zh) | 考虑频率相关性材料的非线性特征值拓扑优化方法及系统 | |
| Finotto et al. | Discrete topology optimization of planar cable-truss structures based on genetic algorithms | |
| Talaslioglu | Design Optimization of Tubular Lattice Girders | |
| Jarallah | POST-PEAK RESPONSE AUTOMATIC SOLUTIONS IN STRUCTURAL ENGINEERING PROBLEMS–A REVIEW | |
| CN116595847A (zh) | 一种基于神经网络的复材固化应力计算方法 | |
| Eamon | Reliability-Based Design Optimization | |
| Xiong et al. | Modified proportional topology optimization algorithm for multiple optimization problems | |
| Cichon et al. | Large displacement analysis of elastic-plastic trusses with unstable bars |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant |