CN112380667A - 疏水多级微结构的纤维丛拓扑优化方法 - Google Patents

疏水多级微结构的纤维丛拓扑优化方法 Download PDF

Info

Publication number
CN112380667A
CN112380667A CN202011075569.6A CN202011075569A CN112380667A CN 112380667 A CN112380667 A CN 112380667A CN 202011075569 A CN202011075569 A CN 202011075569A CN 112380667 A CN112380667 A CN 112380667A
Authority
CN
China
Prior art keywords
equation
liquid
secondary structure
microstructure
vapor interface
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.)
Granted
Application number
CN202011075569.6A
Other languages
English (en)
Other versions
CN112380667B (zh
Inventor
张卫红
邓永波
刘震宇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202011075569.6A priority Critical patent/CN112380667B/zh
Publication of CN112380667A publication Critical patent/CN112380667A/zh
Application granted granted Critical
Publication of CN112380667B publication Critical patent/CN112380667B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/12—Cloth

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)
  • Nonwoven Fabrics (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及一种疏水多级微结构的纤维丛拓扑优化方法，通过对固体表面多级微结构支撑的液汽界面建模，计算与该界面凸起体积极小值对应的微结构几何构型，实现疏水多级微结构的纤维丛拓扑优化方法。该方法克服了目前表面微结构设计凭借研究人员直觉和仿生手段的局限性，提高了设计方法的适用性、灵活性和效率。

Description

疏水多级微结构的纤维丛拓扑优化方法
技术领域
本发明属于疏水多级微结构技术领域，本发明涉及疏水固体表面上周期性多级微结构的几何构型拓扑优化技术领域，特别涉及一种实现多级结构逆向设计的纤维丛拓扑优化方法。
背景技术
润湿是表面化学的一个重要研究方向。固体表面的润湿能力可以按疏水性，亲水性，疏油性，亲油性，双疏性和双亲性进行分类。在具有特殊润湿能力的人工表面微结构研究中，最新成果主要是针对疏水微结构。
对于一个粗糙平面，润湿取决于固体表面自由能和表面微结构构型。粗糙平面上的润湿现象具有两种不同的模态，分别是Wenzel模态和Cassie-Baxter模态。在Wenzel 模态中，液体完全填充粗糙平面上的微结构。在Cassie-Baxter模态中，液汽界面和微结构之间存在汽囊，且汽囊会被束缚在微结构中；此时，粗糙表面相当于一个固体和气体的复合表面，因而粗糙表面疏水性得到有效的增强。当增大施加于液体的压强时，固体和液体的接触模态可以从Cassie-Baxter模态转变为Wenzel模态。该转变过程中，液体会填满粗糙表面的微结构，进而导致疏水性的降低。因此，固体表面上微结构几何构型的合理性决定了微结构抵抗Cassie-Baxter模态向Wenzel模态转变的能力。该抵抗能力可通过Cassie-Baxter模态中液汽界面凸起的体积衡量；并且体积越小，微结构抵抗模态转变的能力越高。
固体表面微结构主要包括层状和多级微结构两类。层状微结构可视为对多级微结构的近似，且多级微结构能够支持更多处于亚稳态的Cassie-Baxter模态。因此，多级微结构的合理设计受到了更广泛的关注。多级微结构具有几何上的纤维丛构型。目前，多级微结构的设计普遍采用仿生方法，该方法通过模拟自然界植物表面上的纹理微结构获取相应的几何构型。该仿生学设计方法缺乏严格的数学模型，并且设计效率较低。因此，发展严格、高效的多级微结构逆向设计方法具有很强的必要性。
在已有的结构逆向设计方法中，拓扑优化被认为是最为有效的。目前，密度法和水平集法是两种普遍采用的拓扑优化实现方法。相比于水平集法，密度法在初值依赖性、收敛效率和多约束处理等方面具有显著的优势。因此，以下采用密度法给出并实现疏水多级微结构的纤维丛拓扑优化方法。
发明内容
要解决的技术问题
为了解决现有技术中表面微结构设计凭借研究人员直觉和仿生手段的局限性，本发明提出一种通过对固体表面多级微结构支撑的液汽界面建模，计算与该界面凸起体积极小值对应的微结构几何构型，实现疏水多级微结构的纤维丛拓扑优化方法。
技术方案
一种疏水多级微结构的纤维丛拓扑优化方法，其特征在于步骤如下：
步骤1：设计变量定义，包括主机构设计变量和二级结构设计变量：
所述的主结构设计变量：定义取值于[0,1]的设计变量zm，用以表达主结构表面Σ；为获得光滑的主结构表面，对zm进行以下滤波操作：
Figure RE-GDA0002896625040000021
对滤波处理后的主结构设计变量zf做进一步的滤波，以获得主结构上点的竖坐标zs并控制主结构高度：
Figure RE-GDA0002896625040000022
其中，Ω是zm的定义域，对应于归一化的固体表面周期性剖分单元；
Figure RE-GDA0002896625040000023
是二维平面内的梯度算子；i和j分别是x和y轴的方向向量；rm是滤波半径，其值为常数；bz是控制主结构高度的参数，该参数取非负值；主结构表面Σ上的法向量为
Figure RE-GDA0002896625040000031
主结构表面Σ称为主结构流形，该流形为二级结构的可变设计区域；
所述的二级结构设计变量：定义取值于[0,1]的变量γ为主结构上二级结构的设计变量；二级结构分布于主结构表面Σ上；为控制二级结构的特征尺寸和移除其中的灰度区域，将γ进行如下的滤波和投影处理：
Figure RE-GDA0002896625040000032
Figure RE-GDA0002896625040000033
其中，
Figure RE-GDA0002896625040000034
是Σ上的切向梯度算子；γp为投影后的二级结构设计变量，将其命名为材质密度；rf是滤波半径，其值为常数；ξ和β是投影参数，其值由数值实验获得；Σ上的切向梯度算子
Figure RE-GDA0002896625040000035
与Ω上的梯度算子
Figure RE-GDA0002896625040000036
之间的关系为
Figure RE-GDA0002896625040000037
因为
Figure RE-GDA0002896625040000038
依赖于zs，该算子对zs的一阶变分满足
Figure RE-GDA0002896625040000039
其中，
Figure RE-GDA00028966250400000310
是zs的试函数；
步骤2：基于以上的设计变量定义，所给出的拓扑优化方法将通过同时演化两套设计变量实现主结构和二级结构的最佳匹配；几何上，最佳匹配后的两套设计变量构成纤维丛(Σ×γp(Σ),Σ,proj1,γp(Σ))，其中Σ为该纤维丛的底流形，γp:Σ→[0,1]为纤维，proj1为满足以下条件的自然映射：
Figure RE-GDA0002896625040000041
因此，所给出的是一种纤维丛拓扑优化方法，并且该方法通过最佳匹配所定义的主结构和二级结构设计变量实现；
当疏水多级微结构上支撑的液汽界面施加的压强逐渐增大时，处于的 Cassie-Baxter模态的液汽界面曲率随之逐渐增大，Cassie-Baxter模态的液汽界面最终完全由二级结构支撑，此时对应的亚稳态称为Cassie-Baxter模态的终态；此后， Cassie-Baxter模态开始向Wenzel模态转变；因此，疏水多级微结构的拓扑优化目标是提高Cassie-Baxter模态的终态稳定性；
基于不相溶两相流体的界面自由能极小原理，多级微结构支撑的Cassie-Baxter终态液汽界面是常平均曲率曲面；该曲面由Young-Laplace方程描述，其物理意义在于液汽界面上毛细压和表面张力之间的平衡；在归一化尺度下，无量纲Young-Laplace 方程为
Figure RE-GDA0002896625040000042
其中，
Figure RE-GDA0002896625040000043
是液汽界面相对于Σ的无量纲化位移；d0是液汽界面原始位移的数量级；κ是对应于Σ平均曲率的静压强，将其命名为底流形压强；
Figure RE-GDA0002896625040000044
是无量纲化的表面张力；σ和P分别为表面张力和两相界面处的静压强；常量
Figure RE-GDA0002896625040000045
是无量纲化的液体表面张力；为保证无量纲Young-Laplace方程解的唯一性，液汽界面的边界设置为
Figure RE-GDA0002896625040000046
由滤波方程2可得
Figure RE-GDA0002896625040000051
进而底流形压强κ可变换为如下形式：
Figure RE-GDA0002896625040000052
在疏水多级微结构拓扑优化中，对无量纲化表面张力
Figure RE-GDA0002896625040000059
进行基于材质密度的插值，并通过材质密度对底流形压强κ进行惩罚，具体如下：
Figure RE-GDA0002896625040000053
Figure RE-GDA0002896625040000054
其中，
Figure RE-GDA0002896625040000055
是液固界面的无量纲化表面张力；q是用于调整材质插值和惩罚方程凹凸性的参数；pκ是底流形压强κ的惩罚因子，其最大和最小值分别为pκ,max和pκ,min；在表面张力的材质插值方程中，
Figure RE-GDA0002896625040000056
的理论值为正无穷大；在数值执行过程中，该参数取为足够大的正数，以同时保证数值计算的收敛性和液固界面的逼近精度；惩罚因子 pκ的作用是消除二级结构上的底流形压强，并在液汽界面上保持该压强；
基于上述底流形压强变换、表面张力材质插值和底流形压强惩罚，方程8中的无量纲Young-Laplace方程可变换为
Figure RE-GDA0002896625040000057
对Cassie-Baxter模态的终态稳定性采用如下最小二乘形式进行度量
Figure RE-GDA0002896625040000058
其中，|Σ|是底流形Σ的面积，具体表达式为
|Σ|＝∫Σ1dΣ (14)
方程13中的终态稳定性度量等价于单位底流形面积上液汽界面凸起体积的平方；因此，疏水多级微结构的拓扑优化目标设置为极小化方程13中的终态稳定性度量；
步骤3：建立如下拓扑优化问题，以实现主结构流形和二级结构版图的最佳匹配，进而实现纤维丛结构拓扑优化：
Figure RE-GDA0002896625040000065
for fiber bundle(Σ×γp(Σ),Σ,proj1,γp(Σ)),
to minimize
Figure RE-GDA0002896625040000061
with
Figure RE-GDA0002896625040000062
constrained by
Figure RE-GDA0002896625040000063
在方程15中，fd是二级结构于主结构流形上的占空比，其表达式为
Figure RE-GDA0002896625040000064
f0是设计者给定的占空比；
为求解方程15中的拓扑优化问题，采用基于梯度信息的迭代算法，其中的液汽界面终态稳定性度量和二级结构占空比的梯度信息可通过伴随分析法获得；基于伴随分析法，液汽界面终态稳定性度量的敏度为：
Figure RE-GDA0002896625040000071
其中δ为一阶变分算子，δ(J|Σ|2)为J|Σ|2的伴随敏度：
Figure RE-GDA0002896625040000072
其中，δzm和δγ分别是zm和γ的一阶变分；zfa和γfa分别是zf和γf的伴随变量；
Figure RE-GDA00028966250400000723
和
Figure RE-GDA00028966250400000722
分别是
Figure RE-GDA00028966250400000721
和
Figure RE-GDA0002896625040000075
的对偶空间；
Figure RE-GDA0002896625040000076
和
Figure RE-GDA0002896625040000077
分别是Σ和Ω上无穷光滑函数空间
Figure RE-GDA0002896625040000078
和
Figure RE-GDA0002896625040000079
的闭空间；方程18中的伴随变量由以下伴随方程求得：
find
Figure RE-GDA00028966250400000710
with
Figure RE-GDA00028966250400000711
on
Figure RE-GDA00028966250400000712
satisfying
Figure RE-GDA00028966250400000713
find
Figure RE-GDA00028966250400000714
satisfying
Figure RE-GDA00028966250400000715
find
Figure RE-GDA00028966250400000716
withγfa＝0 on
Figure RE-GDA00028966250400000717
satisfying
Figure RE-GDA00028966250400000718
find
Figure RE-GDA00028966250400000719
with zsa＝0 on
Figure RE-GDA00028966250400000720
satisfying
Figure RE-GDA0002896625040000081
find
Figure RE-GDA0002896625040000082
with zfa＝0 on
Figure RE-GDA0002896625040000083
satisfying
Figure RE-GDA0002896625040000084
其中，
Figure RE-GDA0002896625040000085
κa和zsa分别是
Figure RE-GDA0002896625040000086
κ和zs的伴随变量；
Figure RE-GDA0002896625040000087
和
Figure RE-GDA0002896625040000088
分别是
Figure RE-GDA0002896625040000089
κa，γa，zsa和zfa的试函数；方程17中，|Σ|的伴随敏度为：
Figure RE-GDA00028966250400000810
其中的伴随变量zfa由以下伴随方程求得：
find
Figure RE-GDA00028966250400000811
with zsa＝0 on
Figure RE-GDA00028966250400000812
satisfying
Figure RE-GDA00028966250400000813
find
Figure RE-GDA00028966250400000814
with zfa＝0 on
Figure RE-GDA00028966250400000815
satisfying
Figure RE-GDA00028966250400000816
二级结构占空比的敏度为：
Figure RE-GDA00028966250400000817
其中，δ|Σ|由方程24，25和26求得；fd|Σ|的伴随敏度为
Figure RE-GDA0002896625040000091
方程28中的伴随变量γfa和zfa由以下方程求得：
find
Figure RE-GDA0002896625040000092
withγfa＝0 on
Figure RE-GDA0002896625040000093
satisfying
Figure RE-GDA0002896625040000094
find
Figure RE-GDA0002896625040000095
with zsa＝0 on
Figure RE-GDA0002896625040000096
satisfying
Figure RE-GDA0002896625040000097
find
Figure RE-GDA0002896625040000098
with zfa＝0 on
Figure RE-GDA0002896625040000099
satisfying
Figure RE-GDA00028966250400000910
在进行伴随分析后，采用以下迭代步骤求解上述变分问题：
(a)求解滤波方程1，2和3，并用方程4投影滤波后的二级结构设计变量；
(b)由方程16计算当前设计变量对应的二级结构占空比；
(c)求解无量纲Young-Laplace方程12，并由方程13和14计算液汽界面终态稳定性度量；
(d)计算方程18中的伴随敏度，其中的伴随变量由方程19，20，21，22，23求得；
(e)计算方程24中的伴随敏度，其中的伴随变量由方程25，26求得；
(f)由方程17计算液汽界面终态稳定性度量的伴随敏度；
(g)计算方程28中的伴随敏度，其中的伴随变量由方程29，30，31求得；
(h)计算方程27中的二级结构占空比伴随敏度；
(i)演化主结构和二级结构的设计变量；
(j)判断是否满足收敛条件，若满足则终止迭代；若不满足则返回(a)。
本发明技术方案更进一步的说：步骤2中在底流形压强的惩罚因子中，pκ,max和pκ,min的取值分别为1和0。
本发明技术方案更进一步的说：步骤2中f0的取值范围为(0,1)。
本发明技术方案更进一步的说：步骤2中终止迭代的条件为：(1)迭代次数达到最大值315；(2)连续5次迭代的设计目标值变化量与5次迭代目标平均值的比值和占空比与f0的偏差都小于10-3。
本发明技术方案更进一步的说：步骤2中迭代过程中，β的初始值取为1，之后经每30次迭代其值翻倍一次；ξ取值0.5。
有益效果
本发明提出的一种疏水固体表面上的多级微结构几何构型的拓扑优化方法，该方法克服了目前表面微结构设计凭借研究人员直觉和仿生手段的局限性，提高了设计方法的适用性、灵活性和效率。由于该方法的主结构设计变量定义了二级结构的可变设计区域，本方法是一种具有可变设计区域的拓扑优化方法。由于多级微结构通过最佳匹配主结构流形和二级结构版图获得，且主结构流形和二级结构版图构成几何上的纤维丛，本方法是一种纤维丛拓扑优化方法。
附图说明
图1为疏水多级微结构图：1a，1b分别为疏水多级微结构的立体和截面示意图。
图2为归一化尺度上多级微结构底流形和二级结构版图，及对应的纤维丛沉降和缩放操作示意图。
图3为主结构设计变量的两次滤波操作示意图。
图4为二级结构设计变量的滤波和投影操作示意图。
图5为疏水多级微结构上Cassie-Baxter模态的受压演化过程示意图。
图6为正三角形相关图：图6a，6b1，6b2，6c1，6c2分别为固体表面的正三角形周期性划分，划分单元上获得的疏水多级微结构纤维丛构型的立体图、俯视图，沉降操作后的纤维丛构型立体图、俯视图。
图7为正四边形相关图：图7a，7b1，7b2，7c1，7c2分别为固体表面的正四边形周期性划分，划分单元上获得的疏水多级微结构纤维丛构型立体图、俯视图，沉降操作后的纤维丛构型立体图、俯视图。
图8为正六边形相关图：图8a1，8c1，8c2，8d1，8d2分别为固体表面的正六边形周期性划分，划分单元上获得的疏水多级微结构纤维丛构型立体图、俯视图，沉降操作后的纤维丛构型立体图、俯视图。
图9为30°顶角等腰三角形相关图：图9a，9b1，9b2，9c1，9c2分别为固体表面30°顶角等腰三角形的周向对称性划分，划分单元上获得的疏水多级微结构纤维丛构型立体图、俯视图，沉降操作后的纤维丛构型立体图、俯视图。
图10为15°顶角等腰三角形相关图：图10a，10b1，10b2，10c1，10c2分别为固体表面15°顶角等腰三角形的周向对称性划分，划分单元上获得的疏水多级微结构纤维丛构型立体图、俯视图，沉降操作后的纤维丛构型立体图、俯视图。
图11为等边四边形相关图：图11a，11b1，11b2，11c1，11c2分别为固体表面的等边四边形手性划分，划分单元上获得的疏水多级微结构纤维丛构型立体图、俯视图，沉降操作后的纤维丛构型立体图、俯视图。
图12为等边九边形相关图：图12a，12b1，12b2，12c1，12c2分别为固体表面的等边九边形手性划分，划分单元上获得的疏水多级微结构纤维丛构型立体图、俯视图，沉降操作后的纤维丛构型立体图、俯视图。
图13为菱形相关图：图13a，13b1，13b2，13c1，13c2分别为基于两个菱形单元的固体表面准周期性划分，划分单元上获得的疏水多级微结构纤维丛构型立体图、俯视图，沉降操作后的纤维丛构型立体图、俯视图。
具体实施方式
现结合实施例、附图对本发明作进一步描述：
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及具体实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，而不构成对本发明的限制。
本发明疏水多级微结构的纤维丛拓扑优化方法包括以下步骤：
固体表面微结构通常都是周期性分布的，并且通过周期性排布一个结构单元可获得疏水固体表面上的周期性微结构。因此，多级微结构的拓扑优化针对结构单元进行 (图1a)。多级微结构的结构单元由主结构和二级结构构成。几何上，这种结构对应着由底流形和纤维构成的纤维丛(图1b)。
多级微结构上的液汽界面采用Young-Laplace描述，其解具有比例缩放的性质，因而以下在归一化尺度上进行多级微结构的拓扑优化，以保证所得结构在一定尺度范围内的适用性和数值求解精度；进而，多级微结构可通过沉降和缩放拓扑优化所得的归一化纤维丛构型获得(图2)。
定义如下的多级微结构拓扑优化设计变量：
主结构设计变量：定义取值于[0,1]的设计变量zm，用以表达主结构表面Σ。为获得光滑的主结构表面，对zm进行以下滤波操作：
Figure RE-GDA0002896625040000121
对滤波处理后的主结构设计变量zf做进一步的滤波，以获得主结构上点的竖坐标zs并控制主结构高度：
Figure RE-GDA0002896625040000122
其中，Ω是zm的定义域，对应于归一化的固体表面周期性剖分单元；
Figure RE-GDA0002896625040000123
是二维平面内的梯度算子；i和j分别是x和y轴的方向向量；rm是滤波半径，其值为常数；bz是控制主结构高度的参数，该参数取非负值；主结构表面Σ上的法向量为
Figure RE-GDA0002896625040000131
主结构表面Σ称为主结构流形，该流形为二级结构的可变设计区域。上述主结构设计变量经两次滤波后获得主结构形貌的过程如图3所示。
二级结构设计变量：定义取值于[0,1]的变量γ为主结构上二级结构的设计变量。二级结构分布于主结构表面Σ上。为控制二级结构的特征尺寸和移除其中的灰度区域，将γ进行如下的滤波和投影处理：
Figure RE-GDA0002896625040000132
Figure RE-GDA0002896625040000133
其中，
Figure RE-GDA0002896625040000134
是Σ上的切向梯度算子；γp为投影后的二级结构设计变量，将其命名为材质密度；rf是滤波半径，其值为常数；ξ和β是投影参数，其值由数值实验获得。Σ上的切向梯度算子
Figure RE-GDA0002896625040000135
与Ω上的梯度算子
Figure RE-GDA0002896625040000136
之间的关系为
Figure RE-GDA0002896625040000137
因为
Figure RE-GDA0002896625040000138
依赖于zs，该算子对zs的一阶变分满足
Figure RE-GDA0002896625040000139
其中，
Figure RE-GDA00028966250400001310
是zs的试函数。上述二级结构设计变量经滤波和投影后获得二级结构版图的过程如图4所示。
基于以上的设计变量定义，所给出的拓扑优化方法将通过同时演化两套设计变量实现主结构和二级结构的最佳匹配。几何上，最佳匹配后的两套设计变量构成纤维丛 (Σ×γp(Σ),Σ,proj1,γp(Σ))，其中Σ为该纤维丛的底流形，γp:Σ→[0,1]为纤维，proj1为满足以下条件的自然映射：
Figure RE-GDA0002896625040000141
因此，所给出的是一种纤维丛拓扑优化方法，并且该方法通过最佳匹配所定义的主结构和二级结构设计变量实现。
当如图1所示的疏水多级微结构上支撑的液汽界面施加的压强逐渐增大时，处于的Cassie-Baxter模态的液汽界面曲率随之逐渐增大。如图5所示，Cassie-Baxter模态的液汽界面最终完全由二级结构支撑，此时对应的亚稳态称为Cassie-Baxter模态的终态。此后，Cassie-Baxter模态开始向Wenzel模态转变。因此，疏水多级微结构的拓扑优化目标是提高Cassie-Baxter模态的终态稳定性。
基于不相溶两相流体的界面自由能极小原理，多级微结构支撑的Cassie-Baxter终态液汽界面是常平均曲率曲面。该曲面由Young-Laplace方程描述，其物理意义在于液汽界面上毛细压和表面张力之间的平衡。在归一化尺度下，无量纲Young-Laplace 方程为
Figure RE-GDA0002896625040000142
其中，
Figure RE-GDA0002896625040000143
是液汽界面相对于Σ的无量纲化位移；d0是液汽界面原始位移的数量级；κ是对应于Σ平均曲率的静压强，将其命名为底流形压强；
Figure RE-GDA0002896625040000144
是无量纲化的表面张力；σ和P分别为表面张力和两相界面处的静压强；常量
Figure RE-GDA0002896625040000145
是无量纲化的液体表面张力；为保证无量纲Young-Laplace方程解的唯一性，液汽界面的边界设置为
Figure RE-GDA0002896625040000146
由滤波方程2可得
Figure RE-GDA0002896625040000151
进而底流形压强κ可变换为如下形式：
Figure RE-GDA0002896625040000152
在疏水多级微结构拓扑优化中，对无量纲化表面张力
Figure RE-GDA0002896625040000153
进行基于材质密度的插值，并通过材质密度对底流形压强κ进行惩罚，具体如下：
Figure RE-GDA0002896625040000154
Figure RE-GDA0002896625040000155
其中，
Figure RE-GDA0002896625040000156
是液固界面的无量纲化表面张力；q是用于调整材质插值和惩罚方程凹凸性的参数；pκ是底流形压强κ的惩罚因子，其最大和最小值分别为pκ,max和pκ,min。在表面张力的材质插值方程中，
Figure RE-GDA0002896625040000159
的理论值为正无穷大；在数值执行过程中，该参数取为足够大的正数，以同时保证数值计算的收敛性和液固界面的逼近精度。在底流形压强的惩罚因子中，pκ,max和pκ,min的取值分别为1和0。惩罚因子pκ的作用是消除二级结构上的底流形压强，并在液汽界面上保持该压强。
基于上述底流形压强变换、表面张力材质插值和底流形压强惩罚，方程8中的无量纲Young-Laplace方程可变换为
Figure RE-GDA0002896625040000157
对图5所示的Cassie-Baxter模态的终态稳定性采用如下最小二乘形式进行度量
Figure RE-GDA0002896625040000158
其中，|Σ|是底流形Σ的面积，具体表达式为
|Σ|＝∫Σ1dΣ (14) 方程13中的终态稳定性度量等价于单位底流形面积上液汽界面凸起体积的平方。因此，疏水多级微结构的拓扑优化目标设置为极小化方程13中的终态稳定性度量。
综上所述，可建立如下拓扑优化问题，以实现主结构流形和二级结构版图的最佳匹配，进而实现纤维丛结构拓扑优化：
Figure RE-GDA0002896625040000165
for fiber bundle(Σ×γp(Σ),Σ,proj1,γp(Σ)),
to minimize
Figure RE-GDA0002896625040000161
with
Figure RE-GDA0002896625040000162
constrained by
Figure RE-GDA0002896625040000163
在方程15中，fd是二级结构于主结构流形上的占空比，其表达式为
Figure RE-GDA0002896625040000164
f0是设计者给定的占空比，其取值范围为(0,1)。
为求解方程15中的拓扑优化问题，采用基于梯度信息的迭代算法，其中的液汽界面终态稳定性度量和二级结构占空比的梯度信息可通过伴随分析法获得。基于伴随分析法，液汽界面终态稳定性度量的敏度为：
Figure RE-GDA0002896625040000171
其中δ为一阶变分算子。δ(J|Σ|2)为J|Σ|2的伴随敏度：
Figure RE-GDA0002896625040000172
其中，δzm和δγ分别是zm和γ的一阶变分；zfa和γfa分别是zf和γf的伴随变量；
Figure RE-GDA0002896625040000173
和
Figure RE-GDA0002896625040000174
分别是
Figure RE-GDA0002896625040000175
和
Figure RE-GDA0002896625040000176
的对偶空间；
Figure RE-GDA0002896625040000177
和
Figure RE-GDA0002896625040000178
分别是Σ和Ω上无穷光滑函数空间
Figure RE-GDA0002896625040000179
和
Figure RE-GDA00028966250400001710
的闭空间。方程18中的伴随变量由以下伴随方程求得：
find
Figure RE-GDA00028966250400001711
with
Figure RE-GDA00028966250400001712
on
Figure RE-GDA00028966250400001713
satisfying
Figure RE-GDA00028966250400001714
find
Figure RE-GDA00028966250400001715
satisfying
Figure RE-GDA00028966250400001716
find
Figure RE-GDA00028966250400001717
withγfa＝0 on
Figure RE-GDA00028966250400001718
satisfying
Figure RE-GDA00028966250400001719
find
Figure RE-GDA00028966250400001720
with zsa＝0 on
Figure RE-GDA00028966250400001721
satisfying
Figure RE-GDA0002896625040000181
find
Figure RE-GDA0002896625040000182
with zfa＝0 on
Figure RE-GDA0002896625040000183
satisfying
Figure RE-GDA0002896625040000184
其中，
Figure RE-GDA0002896625040000185
κa和zsa分别是
Figure RE-GDA0002896625040000186
κ和zs的伴随变量；
Figure RE-GDA0002896625040000187
和
Figure RE-GDA0002896625040000188
分别是
Figure RE-GDA0002896625040000189
κa，γa，zsa和zfa的试函数。方程17中，|Σ|的伴随敏度为：
Figure RE-GDA00028966250400001810
其中的伴随变量zfa由以下伴随方程求得：
find
Figure RE-GDA00028966250400001811
with zsa＝0 on
Figure RE-GDA00028966250400001812
satisfying
Figure RE-GDA00028966250400001813
find
Figure RE-GDA00028966250400001814
with zfa＝0 on
Figure RE-GDA00028966250400001815
satisfying
Figure RE-GDA00028966250400001816
二级结构占空比的敏度为：
Figure RE-GDA00028966250400001817
其中，δ|Σ|由方程24，25和26求得；fd|Σ|的伴随敏度为
Figure RE-GDA0002896625040000191
方程28中的伴随变量γfa和zfa由以下方程求得：
find
Figure RE-GDA0002896625040000192
withγfa＝0 on
Figure RE-GDA0002896625040000193
satisfying
Figure RE-GDA0002896625040000194
find
Figure RE-GDA0002896625040000195
with zsa＝0 on
Figure RE-GDA0002896625040000196
satisfying
Figure RE-GDA0002896625040000197
find
Figure RE-GDA0002896625040000198
with zfa＝0 on
Figure RE-GDA0002896625040000199
satisfying
Figure RE-GDA00028966250400001910
在进行伴随分析后，采用以下迭代步骤求解上述变分问题：
(a)求解滤波方程1，2和3，并用方程4投影滤波后的二级结构设计变量；
(b)由方程16计算当前设计变量对应的二级结构占空比；
(c)求解无量纲Young-Laplace方程12，并由方程13和14计算液汽界面终态稳定性度量；
(d)计算方程18中的伴随敏度，其中的伴随变量由方程19，20，21，22，23 求得；
(e)计算方程24中的伴随敏度，其中的伴随变量由方程25，26求得；
(f)由方程17计算液汽界面终态稳定性度量的伴随敏度；
(g)计算方程28中的伴随敏度，其中的伴随变量由方程29，30，31求得；
(h)计算方程27中的二级结构占空比伴随敏度；
(i)演化主结构和二级结构的设计变量；
(j)判断是否满足收敛条件，若满足则终止迭代；若不满足则返回(a)，其中，终止迭代的条件为：(1)迭代次数达到最大值315；(2)连续5次迭代的设计目标值变化量与5次迭代目标平均值的比值和占空比与f0的偏差都小于10-3。
在上述迭代过程中，β的初始值取为1，之后经每30次迭代其值翻倍一次，ξ取值0.5。
在一个实施例中，采用上述疏水多级微结构的纤维丛拓扑优化方法，将相关参数设置为表1中的值，可分别得到如图6，7，8所示基于正三角形、正四边形和正六边形的三种固体表面轴对称周期性划分单元上的微结构纤维丛构型，如图9，10所示基于顶角分别为30°和15°等腰三角形的两种固体表面周向对称周期性划分单元上的微结构纤维丛构型，如图11，12所示基于等边四边形和等边九边形的两种固体表面手性划分单元上的微结构纤维丛构型，如图13所示基于两个菱形的固体表面准周期性划分单元上的微结构纤维丛构型。
Figure RE-GDA0002896625040000201
表1
基于所介绍的纤维丛拓扑优化方法，通过表1中的参数，求解方程15中的拓扑优化问题，可得如图6，7，8，9，10，11，12，13所示分别具有轴对称性、周向对称性、手性和准周期性的固体表面划分的单元多级微结构几何构型，而且它们可以使用双光子光刻等微纳三维打印类加工工艺制造出来。在缩放操作中，周期性单元的特征尺寸是其缩放系数；在沉降操作中，拉伸距离应该比液汽界面的深度要大，以避免液汽界面接触到微结构的底部而导致的Cassie-Baxter模态崩溃。如图6，7，8，9，10，11， 12，13所示，多级微结构单元的纤维丛构型具有峰状的主结构和密布于主结构顶端的二级结构。
以上所述本发明的具体实施方式，并不构成对本发明保护范围的限定。任何根据本发明的技术构思所作出的各种其他相应的改变与变形，均应包含在本发明权利要求的保护范围内。

Claims (5)

1.一种疏水多级微结构的纤维丛拓扑优化方法，其特征在于步骤如下：
步骤1：设计变量定义，包括主机构设计变量和二级结构设计变量：
所述的主结构设计变量：定义取值于[0,1]的设计变量zm，用以表达主结构表面Σ；为获得光滑的主结构表面，对zm进行以下滤波操作：
Figure FDA0002716538780000011
对滤波处理后的主结构设计变量zf做进一步的滤波，以获得主结构上点的竖坐标zs并控制主结构高度：
Figure FDA0002716538780000012
其中，Ω是zm的定义域，对应于归一化的固体表面周期性剖分单元；
Figure FDA0002716538780000013
是二维平面内的梯度算子；i和j分别是x和y轴的方向向量；rm是滤波半径，其值为常数；bz是控制主结构高度的参数，该参数取非负值；主结构表面Σ上的法向量为
Figure FDA0002716538780000014
主结构表面Σ称为主结构流形，该流形为二级结构的可变设计区域；
所述的二级结构设计变量：定义取值于[0,1]的变量γ为主结构上二级结构的设计变量；二级结构分布于主结构表面Σ上；为控制二级结构的特征尺寸和移除其中的灰度区域，将γ进行如下的滤波和投影处理：
Figure FDA0002716538780000015
Figure FDA0002716538780000016
其中，
Figure FDA0002716538780000021
是Σ上的切向梯度算子；γp为投影后的二级结构设计变量，将其命名为材质密度；rf是滤波半径，其值为常数；ξ和β是投影参数，其值由数值实验获得；Σ上的切向梯度算子
Figure FDA0002716538780000022
与Ω上的梯度算子
Figure FDA0002716538780000023
之间的关系为
Figure FDA0002716538780000024
因为
Figure FDA0002716538780000025
依赖于zs，该算子对zs的一阶变分满足
Figure FDA0002716538780000026
其中，
Figure FDA0002716538780000027
是zs的试函数；
步骤2：基于以上的设计变量定义，所给出的拓扑优化方法将通过同时演化两套设计变量实现主结构和二级结构的最佳匹配；几何上，最佳匹配后的两套设计变量构成纤维丛(Σ×γp(Σ),Σ,proj1,γp(Σ))，其中Σ为该纤维丛的底流形，γp:Σ→[0,1]为纤维，proj1为满足以下条件的自然映射：
Figure FDA0002716538780000028
因此，所给出的是一种纤维丛拓扑优化方法，并且该方法通过最佳匹配所定义的主结构和二级结构设计变量实现；
当疏水多级微结构上支撑的液汽界面施加的压强逐渐增大时，处于的Cassie-Baxter模态的液汽界面曲率随之逐渐增大，Cassie-Baxter模态的液汽界面最终完全由二级结构支撑，此时对应的亚稳态称为Cassie-Baxter模态的终态；此后，Cassie-Baxter模态开始向Wenzel模态转变；因此，疏水多级微结构的拓扑优化目标是提高Cassie-Baxter模态的终态稳定性；
基于不相溶两相流体的界面自由能极小原理，多级微结构支撑的Cassie-Baxter终态液汽界面是常平均曲率曲面；该曲面由Young-Laplace方程描述，其物理意义在于液汽界面上毛细压和表面张力之间的平衡；在归一化尺度下，无量纲Young-Laplace方程为
Figure FDA0002716538780000031
其中，
Figure FDA0002716538780000032
是液汽界面相对于Σ的无量纲化位移；d0是液汽界面原始位移的数量级；κ是对应于Σ平均曲率的静压强，将其命名为底流形压强；
Figure FDA0002716538780000033
是无量纲化的表面张力；σ和P分别为表面张力和两相界面处的静压强；常量
Figure FDA0002716538780000034
是无量纲化的液体表面张力；为保证无量纲Young-Laplace方程解的唯一性，液汽界面的边界设置为
Figure FDA0002716538780000035
由滤波方程2可得
Figure FDA0002716538780000036
进而底流形压强κ可变换为如下形式：
Figure FDA0002716538780000037
在疏水多级微结构拓扑优化中，对无量纲化表面张力
Figure FDA0002716538780000038
进行基于材质密度的插值，并通过材质密度对底流形压强κ进行惩罚，具体如下：
Figure FDA0002716538780000039
Figure FDA00027165387800000310
其中，
Figure FDA00027165387800000311
是液固界面的无量纲化表面张力；q是用于调整材质插值和惩罚方程凹凸性的参数；pκ是底流形压强κ的惩罚因子，其最大和最小值分别为pκ,max和pκ,min；在表面张力的材质插值方程中，
Figure FDA00027165387800000312
的理论值为正无穷大；在数值执行过程中，该参数取为足够大的正数，以同时保证数值计算的收敛性和液固界面的逼近精度；惩罚因子pκ的作用是消除二级结构上的底流形压强，并在液汽界面上保持该压强；
基于上述底流形压强变换、表面张力材质插值和底流形压强惩罚，方程8中的无量纲Young-Laplace方程可变换为
Figure FDA0002716538780000041
对Cassie-Baxter模态的终态稳定性采用如下最小二乘形式进行度量
Figure FDA0002716538780000042
其中，|Σ|是底流形Σ的面积，具体表达式为
|Σ|＝∫Σ1dΣ (14)
方程13中的终态稳定性度量等价于单位底流形面积上液汽界面凸起体积的平方；因此，疏水多级微结构的拓扑优化目标设置为极小化方程13中的终态稳定性度量；
步骤3：建立如下拓扑优化问题，以实现主结构流形和二级结构版图的最佳匹配，进而实现纤维丛结构拓扑优化：
Figure FDA0002716538780000051
for fiber bundle(Σ×γp(Σ),Σ,proj1,γp(Σ)),
Figure FDA0002716538780000052
Figure FDA0002716538780000053
在方程15中，fd是二级结构于主结构流形上的占空比，其表达式为
Figure FDA0002716538780000054
f0是设计者给定的占空比；
为求解方程15中的拓扑优化问题，采用基于梯度信息的迭代算法，其中的液汽界面终态稳定性度量和二级结构占空比的梯度信息可通过伴随分析法获得；基于伴随分析法，液汽界面终态稳定性度量的敏度为：
Figure FDA0002716538780000055
其中δ为一阶变分算子，δ(J|Σ|2)为J|Σ|2的伴随敏度：
Figure FDA0002716538780000061
其中，δzm和δγ分别是zm和γ的一阶变分；zfa和γfa分别是zf和γf的伴随变量；
Figure FDA0002716538780000062
和
Figure FDA0002716538780000063
分别是
Figure FDA0002716538780000064
和
Figure FDA0002716538780000065
的对偶空间；
Figure FDA0002716538780000066
和
Figure FDA0002716538780000067
分别是Σ和Ω上无穷光滑函数空间
Figure FDA0002716538780000068
和
Figure FDA0002716538780000069
的闭空间；方程18中的伴随变量由以下伴随方程求得：
Figure FDA00027165387800000610
Figure FDA00027165387800000611
Figure FDA00027165387800000612
Figure FDA00027165387800000613
Figure FDA00027165387800000614
Figure FDA00027165387800000615
Figure FDA00027165387800000616
Figure FDA0002716538780000071
Figure FDA0002716538780000072
Figure FDA0002716538780000073
其中，
Figure FDA0002716538780000074
κa和zsa分别是
Figure FDA0002716538780000075
κ和zs的伴随变量；
Figure FDA0002716538780000076
和
Figure FDA0002716538780000077
分别是
Figure FDA0002716538780000078
κa，γa，zsa和zfa的试函数；方程17中，|Σ|的伴随敏度为：
Figure FDA0002716538780000079
其中的伴随变量zfa由以下伴随方程求得：
Figure FDA00027165387800000710
Figure FDA00027165387800000711
Figure FDA00027165387800000712
Figure FDA00027165387800000713
二级结构占空比的敏度为：
Figure FDA00027165387800000714
其中，δ|Σ|由方程24，25和26求得；fd|Σ|的伴随敏度为
Figure FDA0002716538780000081
方程28中的伴随变量γfa和zfa由以下方程求得：
Figure FDA0002716538780000082
Figure FDA0002716538780000083
Figure FDA0002716538780000084
Figure FDA0002716538780000085
Figure FDA0002716538780000086
Figure FDA0002716538780000087
在进行伴随分析后，采用以下迭代步骤求解上述变分问题：
(a)求解滤波方程1，2和3，并用方程4投影滤波后的二级结构设计变量；
(b)由方程16计算当前设计变量对应的二级结构占空比；
(c)求解无量纲Young-Laplace方程12，并由方程13和14计算液汽界面终态稳定性度量；
(d)计算方程18中的伴随敏度，其中的伴随变量由方程19，20，21，22，23求得；
(e)计算方程24中的伴随敏度，其中的伴随变量由方程25，26求得；
(f)由方程17计算液汽界面终态稳定性度量的伴随敏度；
(g)计算方程28中的伴随敏度，其中的伴随变量由方程29，30，31求得；
(h)计算方程27中的二级结构占空比伴随敏度；
(i)演化主结构和二级结构的设计变量；
(j)判断是否满足收敛条件，若满足则终止迭代；若不满足则返回(a)。
2.根据权利要求1所述的一种疏水多级微结构的纤维丛拓扑优化方法，其特征在于步骤2中在底流形压强的惩罚因子中，pκ,max和pκ,min的取值分别为1和0。
3.根据权利要求1所述的一种疏水多级微结构的纤维丛拓扑优化方法，其特征在于步骤2中f0的取值范围为(0,1)。
4.根据权利要求1所述的一种疏水多级微结构的纤维丛拓扑优化方法，其特征在于步骤2中终止迭代的条件为：(1)迭代次数达到最大值315；(2)连续5次迭代的设计目标值变化量与5次迭代目标平均值的比值和占空比与f0的偏差都小于10-3。
5.根据权利要求1所述的一种疏水多级微结构的纤维丛拓扑优化方法，其特征在于步骤2中迭代过程中，β的初始值取为1，之后经每30次迭代其值翻倍一次；ξ取值0.5。
CN202011075569.6A 2020-10-09 2020-10-09 疏水多级微结构的纤维丛拓扑优化方法 Active CN112380667B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011075569.6A CN112380667B (zh) 2020-10-09 2020-10-09 疏水多级微结构的纤维丛拓扑优化方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011075569.6A CN112380667B (zh) 2020-10-09 2020-10-09 疏水多级微结构的纤维丛拓扑优化方法

Publications (2)

Publication Number Publication Date
CN112380667A true CN112380667A (zh) 2021-02-19
CN112380667B CN112380667B (zh) 2022-09-09

Family

ID=74581094

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011075569.6A Active CN112380667B (zh) 2020-10-09 2020-10-09 疏水多级微结构的纤维丛拓扑优化方法

Country Status (1)

Country Link
CN (1) CN112380667B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108052729A (zh) * 2017-12-11 2018-05-18 中国科学院长春光学精密机械与物理研究所 超疏水表面微结构的逆向设计方法
CN110008512A (zh) * 2019-03-04 2019-07-12 三峡大学 一种考虑承载特性的负泊松比点阵结构拓扑优化方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108052729A (zh) * 2017-12-11 2018-05-18 中国科学院长春光学精密机械与物理研究所 超疏水表面微结构的逆向设计方法
CN110008512A (zh) * 2019-03-04 2019-07-12 三峡大学 一种考虑承载特性的负泊松比点阵结构拓扑优化方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
庄春刚等: "三维线弹性连续体的结构拓扑优化设计", 《中国机械工程》 *

Also Published As

Publication number Publication date
CN112380667B (zh) 2022-09-09

Similar Documents

Publication Publication Date Title
Seifi et al. Design optimization and additive manufacturing of nodes in gridshell structures
Yang et al. A parameter-free generalized moment limiter for high-order methods on unstructured grids
CN102720295B (zh) 一种基于索穹顶张拉和承载全过程分析的预应力确定方法
CN112417692B (zh) 基于载荷不确定性的材料结构多尺度拓扑优化设计方法
CN104866673A (zh) 一种轴压加筋柱壳的开口补强方法
Feng et al. Stiffness optimization design for TPMS architected cellular materials
CN108763658A (zh) 基于等几何方法的组合薄壁结构固有频率设计方法
CN110210130B (zh) 针对工字梁二维模型的形状优化方法
CN113326582B (zh) 一种基于应力分布的变密度点阵结构及其设计方法
CN111625910B (zh) 一种基于静强度和疲劳约束的柔顺机构拓扑优化方法
CN110555263A (zh) 一种用于曲壳结构优化设计的水平集拓扑优化方法
CN113887095B (zh) 一种基于等几何分析的渐进式结构拓扑优化方法
CN102493569B (zh) 一种建筑结构基于抗震性能的优化方法和系统
CN108829914A (zh) 一种frp结构件的结构与工艺一体化设计方法
CN110321611A (zh) 一种多材料结构拓扑优化方法
Mulani et al. Optimal design of unitized structures using response surface approaches
CN112380667B (zh) 疏水多级微结构的纤维丛拓扑优化方法
CN103353916A (zh) 基于工程的复合材料层合板铺层优化后处理方法
Huynh et al. A polytree-based adaptive scheme for modeling linear fracture mechanics using a coupled XFEM–SBFEM approach
CN114357546B (zh) 梯度随机点阵/多孔结构多尺度优化设计方法
CN109502017B (zh) 一种拓扑优化仿生无人机及其设计方法
CN107515960B (zh) 基于特征建模的循环对称筒形支撑结构拓扑优化设计方法
Xiao et al. Infill topology and shape optimization of lattice‐skin structures
CN108052729B (zh) 超疏水表面微结构的逆向设计方法
CN108446507A (zh) 基于网格质量反馈优化的弹性体网格变形方法

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
GR01 Patent grant