基于模糊控制理论的多物理场耦合计算自适应步长方法
技术领域
本发明涉及多物理场耦合条件下数值仿真技术领域，特别涉及一种基于模糊控制理论的多物理场耦合计算自适应步长方法，为指导工程技术人员利用知识与经验实现复杂多物理场耦合问题的高效求解提供了可行的处理办法。
背景技术
在航空航天、桥梁、船舶、化工、能源等领域往往涉及多物理场耦合问题。多物理场耦合仿真是高度非线性、时变性、强耦合的过程，计算量大、周期长、资源负担重，导致研发成本较高。考虑到耦合求解复杂性，人们往往采取单物理场各自求解、耦合界面数据传递的松耦合方法。松耦合方法由于物理时间与仿真时间的不统一，一般只有一阶的时间精度，为了提高求解精度必须将时间步长取很小，导致松耦合方法计算效率显著下降。此外，传统的多场耦合分析与结构设计中，工程技术人员积累下来的知识与经验难以在耦合求解中得到充分利用。因此，充分利用积累的知识与经验，实现精度与效率兼顾的多场耦合分析，具有显著的工程实用价值。
在数值仿真领域，自适应步长技术往往是提高计算效率的有效方法。多物理场耦合求解中常用的基于经验的控制方法、经典控制方法等均有不足。前者设计方便，但步长变化不平滑，可能影响数值解的平滑性及计算稳定性；后者步长变化平滑，但为获得最佳控制性能，需要建立被控对象的精确模型，且难以利用以往知识经验。模糊控制是一种智能控制方法，它利用多值模糊逻辑和人工智能要素(简化推理原则)来模仿人的思维，设计方便、不需要建立被控对象精确模型，可解决非线性、时变性、强耦合系统的控制问题，易于实现基于知识甚至语义描述的控制规律。因此模糊控制方法适用于复杂多物理场耦合分析的变步长控制。
综上，将以往知识和经验转化为控制规律、将模糊控制理论用于多物理场耦合分析的自适应步长控制，可以有效提高仿真效率，再结合各类高精度耦合求解格式，即可实现精度与效率兼顾的多场耦合分析。但是目前多场耦合计算领域，对基于模糊控制理论的变步长方法研究相对较少，基础薄弱。
发明内容
为了解决上述技术问题，克服现有技术的不足，本发明提供一种基于模糊控制理论的多物理场耦合计算自适应步长方法，充分考虑实际工程中多物理场耦合计算精度与效率难以兼顾、积累的知识经验难以重用的现状，以模糊控制理论为基础，建立一种多物理场耦合分析的自适应步长方法。所提方法相比传统定步长方法，可在保证多场精度的同时实现耦合问题的高效计算。本发明将重点探讨模糊控制理论框架下多物理场耦合求解的自适应步长方法，根据多场耦合界面数据获得控制系统输入，采用专家经验法建立模糊规则库，通过建立通用的计算流程，为多物理场耦合的高效计算提供技术保障。
本发明采用的技术方案为：一种基于模糊控制理论的多物理场耦合计算自适应步长方法，该方法实现步骤如下：
第一步：进行模糊控制系统输入、输出变量的选取。对于多物理场求解的场间数据传递过程，以S＝{Ii|i＝1,2,…,α}表征交换的数据组成的集合，α为数据类型总数，Ii为第i类数据矢量，Ii中元素按照结点编号或自由度进行排序，即：
其中ji表示第i类数据的总结点数或自由度数：若第i类数据为矢量，则ji表示自由度数；若第i类数据为标量，则ji表示结点数。
利用
其中Ii,n表示第n步第i类数据矢量，||·||为矢量范数。利用Δen＝|en-en-1|表示第n步计算后en的变化量。利用e0、Δe0分别表示en与Δen的容差，后续计算过程中需保证下述关系式成立：
en≤e0，Δen≤Δe0
模糊控制系统输入变量定义为E、ΔE，第n步计算输入变量具体表达式如下：
模糊控制系统输出变量定义为步长调节因子k，与步长满足下述关系：
hn+1＝khn
其中hn表示第n步计算采用的步长。
第二步：利用第一步提出的输入变量E、ΔE，以及输出变量k，进行模糊化设计、模糊推理和精确化设计。将输入变量E、ΔE的论域设置为U＝[0,1]，输出变量k的论域设置为V＝[0,k0]，k0为步长调节因子上界。对E、ΔE、k进行模糊集划分，隶属函数分别为
根据以往多场耦合计算的知识与经验，制定模糊控制规则。根据仿真经验，如果某个多场耦合问题选取的模糊控制变量对耦合计算结果影响较大，则该变量可以划分为更多模糊集，否则划分为较少模糊集；根据仿真经验，如果某个多场耦合问题各场变化幅度不大，且按常规定步长方法精度足够，则在所建立模糊控制规则中，输入变量的各个模糊集可以对应较“大”的输出变量模糊集，否则对应较“小”的输出变量模糊集。利用a、b分别表示某个输入变量或输出变量，
μR,Zadeh(a,b)＝[μA(a)∧μB(b)]∨[1-μA(a)]
μR,Mandani(a,b)＝[μA(a)∧μB(b)]
其中∧表示取最小值，∨表示取最大值。采用直积建立模糊关系矩阵
其中U1×U2＝{(E,ΔE)|E∈U1,ΔE∈U2}。针对每条模糊规则进行模糊推理，求出推理结果的并集，得到总输出，表达式为：
其中
用
从而，本发明完成了基于知识与经验的多物理场耦合模糊控制器的设计。
第三步：以第一步至第二步定义的模糊控制器作为步长控制器，设定基本参数：仿真总时间ttot，仿真总步数ntot和步长调节因子上界的最大容许值k′，且k′满足：
k′＜0.1ntot
第四步：以第一步至第二步定义的模糊控制器作为步长控制器，取仿真时间ttot0和仿真步数ntot0进行试算：分别进行定步长计算和基于模糊控制的变步长计算，将两种方法计算结果的差异定义为η，将容许的差异的上界定义为η0，在η≤η0的前提下，确定初始步长h0、步长调节因子上界k0(k0≤k′)和容差e0、Δe0。
第五步：以第一步至第二步定义的模糊控制器作为步长控制器，制定特例处理措施。
将步长取值区间[hmin,hmax]定义为：
在计算中，若khn-1∈[hmin,hmax]，则取hn＝khn-1；若khn-1≤hmin，则取hn＝hmin；若khn-1≥hmax，则取hn＝hmax。
若当前时间步容差超出许可值，如果当前步长hn大于hmin，则取二者的平均值作为新步长，重新计算当前时间步；若当前步长hn接近hmin，则扩大步长区间，并取hn与新的最小步长hmin,new的平均值代替hn，重新计算当前时间步。一旦满足容差要求，继续进行下一个时间步计算，并将步长取值区间重置为初始区间。
第六步：以第一步至第二步定义的模糊控制器作为步长控制器，确定计算流程。以tn为第n步计算后的总计算时间。以h1＝h0为初始步长，进行第n＝1步计算，根据第1步计算中传递的数据S1＝{Ii,1|i＝1,2,…,α}及初始条件，得到e1，完成第1步计算后令t1＝h1；以h2＝h1为步长进行第n＝2步计算，根据第1、2步计算中传递的数据S1＝{Ii,1|i＝1,2,…,α}、S2＝{Ii,2|i＝1,2,…,α}及e1，得到e2与Δe2，完成第2步计算后令t2＝t1+h2；根据第一步至第二步定义的模糊控制器计算第n(n≥3)步步长hn，即：
hn＝khn-1
进行第n步计算，得到en与Δen，通过第五步措施处理异常情况，完成第n步计算后令tn＝tn-1+hn；每完成一次计算，判断n≥ntot或t≥ttot是否成立，若成立则终止计算。
有益效果：
本发明提供了将积累的知识与经验用于多物理场耦合计算自适应步长的新方法，弥补与完善了传统定步长方法以及基于经典控制理论的自适应步长方法的局限性。所构建的基于模糊控制的变步长方法，一方面使用的模糊控制器设计方便，无需建立复杂系统的模型，另一方面能够将工程人员的仿真经验转化为控制规则，实现了知识经验的重用，为实现复杂多物理场耦合的高效计算奠定了一定的基础。
附图说明
图1是本发明建立的基于模糊控制理论的多物理场耦合计算自适应步长方法流程图；
图2是本发明定义的模糊控制系统结构图；
图3是本发明定义的模糊控制器输入数据来源及步长计算方法示意图；
图4是本发明实施例中二元机翼气动弹性模型示意图；
图5(a)是本发明实施例中输入变量与输出变量E的隶属函数图；
图5(b)是本发明实施例中输入变量与输出变量ΔE的隶属函数图；
图5(c)是本发明实施例中输入变量与输出变量k的隶属函数图；
图6(a)是本发明实施例中定步长方法及自适应步长方法所得刚心沉浮位移与绕刚心俯仰角变化曲线图：刚心沉浮位移曲线图；
图6(b)是本发明实施例中定步长方法及自适应步长方法所得刚心沉浮位移与绕刚心俯仰角变化曲线图：绕刚心俯仰角变化曲线图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅为本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域的普通技术人员在不付出创造性劳动的前提下所获得的所有其他实施例，都属于本发明的保护范围。
如图1所示，本发明提出了一种基于模糊控制理论的多物理场耦合计算自适应步长方法，包括以下步骤：
(1)选取模糊控制系统输入、输出变量。对于多物理场求解的场间数据传递过程，以S＝{Ii|i＝1,2,…,α}表征交换的数据组成的集合，α为数据类型总数，Ii为第i类数据矢量，Ii中元素按照结点编号或自由度进行排序，即：
其中ji表示第i类数据的总结点数或自由度数：若第i类数据为矢量，则ji表示自由度数；若第i类数据为标量，则ji表示结点数。
利用
其中Ii,n表示第n步第i类数据矢量，||·||为矢量范数。若||·||取欧几里得范数，
利用Δen＝|en-en-1|表示第n步计算后en的变化量。利用e0、Δe0分别表示en与Δen的容差，后续计算过程中需满足下述关系式：
en≤e0，Δen≤Δe0
模糊控制系统输入变量定义为E、ΔE，第n步计算输入变量具体表达式如下：
模糊控制系统输出变量定义为步长调节因子k，与步长满足下述关系：
hn+1＝khn
其中hn表示第n步计算采用的步长。
(2)利用第一步提出的输入变量E、ΔE，以及输出变量k，进行模糊化设计。将输入变量E、ΔE的论域设置为U＝[0,1]，输出变量k的论域设置为V＝[0,k0]，k0为步长调节因子上界。对E、ΔE、k进行模糊集划分，隶属函数分别为
根据以往多场耦合计算的知识与经验，制定模糊控制规则。根据仿真经验，如果某个多场耦合问题选取的模糊控制变量对耦合计算结果影响较大，则该变量可以划分为更多模糊集，否则划分为较少模糊集；根据仿真经验，如果某个多场耦合问题各场变化幅度不大，且按常规定步长方法精度足够，则在所建立模糊控制规则中，输入变量的各个模糊集可以对应较“大”的输出变量模糊集，否则对应较“小”的输出变量模糊集。利用a、b分别表示某个输入变量或输出变量，
μR,Zadeh(a,b)＝[μA(a)∧μB(b)]∨[1-μA(a)]
μR,Mandani(a,b)＝[μA(a)∧μB(b)]
其中∧表示取最小值，∨表示取最大值。采用直积建立模糊关系矩阵
其中U1×U2＝{(E,ΔE)|E∈U1,ΔE∈U2}。针对每条模糊规则进行模糊推理，以Mandani方法为例，第s条模糊规则的推理结果为：
其中
其中ns为模糊规则总数。
用
从而，本发明完成了基于知识与经验的多物理场耦合模糊控制器的设计，为后续多场耦合的自适应步长计算提供了基础。该模糊控制系统的结构图见图2。
(3)以第一步至第二步定义的模糊控制器作为步长控制器，设定基本参数：仿真总时间ttot，仿真总步数ntot和步长调节因子上界的最大容许值k′，且满足：
k′＜0.1ntot
(4)以第一步至第二步定义的模糊控制器作为步长控制器，取仿真时间ttot0和仿真步数ntot0进行试算，确定初始步长h0、步长调节因子上界k0(k0≤k′)和容差e0、Δe0。其中ttot0和ntot0应满足如下关系：
ttot0≤0.01ttot,ntot0≤0.01ntot
令初始步长h0为：
初始步长调节因子上界为k0＝2。分别进行定步长计算和基于模糊控制的变步长计算，将两种方法计算结果的差异定义为η。
进行试算，并比较定步长与自适应步长结果的偏差η，将容许的差异上界定义为η0。以η0＝0.5％为例，若k0＝2且η≥0.5％，则减小初始步长h0，并重新计算；若η＜0.5％，则增大步长调节因子上界，例如：
k0+1→k0
并重新计算；当出现下述情况时终止计算：
1.k0＝k′；
2.取k0时η＜0.5％，取k0+1时η≥0.5％。
计算终止后，得到最优步长调节因子k0。
根据第一步的定义，利用最后一次试算中每一时间步的物理场变化量ei,Δei定义容差e0、Δe0：
(5)以第一步至第二步定义的模糊控制器作为步长控制器，制定异常处理措施。
将步长取值区间[hmin,hmax]定义为：
在后续计算中，若khn-1∈[hmin,hmax]，则取hn＝khn-1；若khn-1≤hmin，则取hn＝hmin；若khn-1≥hmax，则取hn＝hmax。
若当前时间步容差超出许可值(e＞e0,Δe＞Δe0)，如果步长hn未达到最小步长hmin，则取二者的平均值为新步长hn,new，即：
并重新计算该时间步；若步长hn十分接近hmin，则扩大步长取值区间，定义小参数m0，0＜m0＜1，将新的步长上下界hmax,new、hmin,new分别定义为：
hmin,new＝hmin/(1+m0),hmax,new＝hmax(1+m0)
再取步长hn与新的最小步长hmin,new的平均值为新步长hn,new，即：
并重新计算该时间步；一旦满足容差要求(e≤e0,Δe≤Δe0)，则将步长取值区间重置为初始区间[hmin,hmax]，继续下一个时间步计算。
(6)以第一步至第二步定义的模糊控制器作为步长控制器，确定计算流程，如图3所示。以tn表示第n步计算后的总计算时间。
1.以h1＝h0为初始步长，进行第n＝1步计算，根据第1步计算中传递的数据S1＝{Ii,1|i＝1,2,…,α}及初始条件计算e1，完成第一步计算后令t1＝h1；
2.以h2＝h1为步长进行第n＝2步计算，根据第1、2步计算中传递的数据S1＝{Ii,1|i＝1,2,…,α}、S2＝{Ii,2|i＝1,2,…,α}及e1，计算e2与Δe2，完成第2步计算后令t2＝t1+h2；
3.第n(n≥3)步步长根据第一步至第二步定义的模糊控制器计算，即
hn＝khn-1
计算en与Δen，通过第五步措施处理异常情况，完成第n步计算后令tn＝tn-1+hn；
4.每完成一次计算，判断n≥ntot或t≥ttot是否成立，若成立则终止计算。
实施例：
为了更充分地了解该发明的特点及其对工程实际的适用性，本发明针对如图4所示的二元机翼气动弹性模型进行基于模糊控制理论的多物理场耦合自适应步长计算。该机翼刚心为E，刚心处固定一个刚度为Kh的线弹簧及一个刚度为Kα的扭转弹簧，气动中心到刚心距离为e。机翼弦长2b，机翼面积S，来流速度为v，来流密度为ρ。单位展长的机翼质量为m，单位展长机翼对刚心的质量静矩为Sα，单位展长机翼对刚心的质量惯性矩为Iα，机翼沉浮运动与俯仰运动的非耦合固有频率分别为ωh与ωα。取刚心的沉浮位移h以及绕刚心俯仰角α为广义坐标，
本实施例中，上述参数取值列于表1，计算初始条件列于表2。取传递数据集合
表1
表2
表3
分别采用四阶Runge-Kutta算法求解该问题，刚心沉浮位移曲线与绕刚心俯仰角变化曲线如图6，定步长方法及自适应步长方法所得刚心沉浮位移与绕刚心俯仰角变化曲线图：图中三角形为自适应步长计算结果，实线为定步长计算结果，(a)刚心沉浮位移曲线图；(b)绕刚心俯仰角变化曲线图。自适应步长求解时，给定的初始步长小于定步长；由于变量变化缓慢，若干次迭代后自适应步长大于定步长，且计算精度与定步长接近，但总积分步数减少，因此计算效率提高。
综上所述，本发明提出了一种基于模糊控制理论的多物理场耦合计算自适应步长方法。该方法将场间传递数据处理后作为模糊控制系统的输入变量，设定输入、输出变量的论域，划分模糊集并设置隶属函数；进而，基于耦合计算的知识与经验制定模糊控制规则，通过模糊推理得到模糊控制系统输出；之后，设置步长取值区间与超出容差处理措施；最后，以模糊控制器为步长控制器，确定计算流程，并设置计算总步数或总时间作为终止条件，完成了多物理场耦合计算步长的自适应模糊控制过程。
本发明未详细阐述部分属于本领域技术人员的公知技术。
以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制；其可扩展应用于步长模糊控制的多物理场耦合计算领域，凡采用等同变换或者等效替换而形成的技术方案，均落在本发明权利保护范围之内。
