我们知道一般仿真分析是将物理模型转化为数学模型再进行求解。通常将物理模型转化为数学模型这一步骤会由仿真软件完成,用户只需要选择相应的模块,并对数学方程的某些参数进行完善即可。目前COMSOL内置了电磁模块、流体流动&传热模块,结构力学&声学模块以及化工等模块,能够满足绝大多数的仿真需求。但是目前COMSOL没有预置所有能想象到的物理场方程,或者现有的模块我们需要太多的修正才能够使用,这时候,选择基于方程的建模方式是一种可行的方案。COMSOL的自定义方程可以与任何其他内置的物理场接口耦合,以获得强大且灵活的建模能力,又能避免用户编写繁琐的程序代码。 我们在使用COMSOL时,有四种建模方式,最方便的是用软件预置的物理场接口进行计算,其二是以预置物理场接口建模并添加修正项,第三是使用预置物理场接口和自定义方程配合使用,第四种也是比较困难的是直接根据原理以方程接口建模。这四种难度依次增大,灵活度逐渐增加。作为COMSOL仿真最后的大招,基于方程的建模形式可以帮助我们建立任意方程,帮助我们实现用内置物理场接口难以实现的数值模拟。 在COMSOL中,我们常用的自定义方程接口包括偏微分方程、经典偏微分方程以及常微分和微分代数方程,下表列出了各方程形式和方程描述。
系数形式偏微分方程是最基本的方程形式,通过控制方程各项的系数,一般形式和波形式微分方程、经典偏微分方程等均可以由系数型偏微分方程衍生得到。弱形式偏微分方程没有指定具体的微分方程形式,其原因是弱形式能够适应所有的微分方程,我们只需要把方程及其边界条件通过一定的变换变为弱形式即可,我们之前通过案例展示过微分方程转化为弱形式的方法。 在结束不久的正月闹新春活动中,敲锣打鼓是最常见的一种活动,那么有没有想过一锤子敲下去,鼓面是怎么振动的呢?我们今天就可以用COMSOL基于方程建模来一探究竟。鼓面可以看做薄膜模型,我们知道薄膜自由振动的方程可以写为: 其中W表示薄膜的振动挠度。这里我们需要的振动模态,也就是在频域对上述方程进行求解,通过傅里叶变换,上述方程变换为: 上式也即COMSOL中波动方程在频域下的表现形式。 我们首先建立几何模型,在COMSOL中绘制一个半径为1.5m的圆面。 对比我们推倒的方程形式,在COMSOL中完善模型,源项f为0,质量系数为1。我们推倒方程中的a2为张力T和材料密度ρ的比值,取T = 23000N/m,密度ρ = 7.805kg/m2,计算得到c = -2946.8,由于采用了统一的SI单位制,我们可以不在COMSOL中设置变量单位。 对于边界条件,鼓面边界上挠度为0,也即狄利克雷边界条件,指定挠度值为0。 求解器选用特征值求解器对鼓面的前六阶特征值进行求解,下面的图中展示出了前六阶特征值及其振动模态,计算结果和理论计算结果基本一致(通过点击文末阅读原文了解理论计算过程)。仿真结果表明,由于人敲打鼓面的频率远小于一阶频率,鼓面主要的振动形式是一阶振动构型也即整体的上下振动,几乎不存在高阶的振动模态(难以激起且高阶阻尼更大)。 系数形式偏微分方程是包含项最全面的微分方程,通过它我们可以建立各种形式的微分方程,热方程当然也包含在内。我们也可以直接用经典偏微分方程->热方程来建模,此处为了演示系数形式偏微分方程的用法,我们暂不采用热方程。 我们知道固体中的热方程具有以下的形式: 其中α为热扩散系数。对比系数形式偏微分方程,留下的仅有堆积储存项、扩散项和源项,因此我们只需要对da,c,f赋以恰当的值,其他系数均设置为0即可。 几何模型如图所示,中间为热源,边界条件为两端温度恒定,其他边界热通量为0。 由于不同域的方程参数不同,模型中建立了两个系数形式偏微分方程,分别指定热源区域和其他两部分区域的方程参数。取c=1,da=1,热源域f=1,其他域f=0。 下面的图形展示了随时间变化求解域中的温度变化情况。 接下来我们考虑,随着时间的增长,热能累积是多少?这就需要偏微分方程和常微分方程联合求解。 首先,我们需要对u在空间上进行积分,然后再在时间上积分。 空间积分在软件中可以很简单的操作,首先定义一个积分算子intop1,积分域为所有域。然后定义一个变量U,使其等于intop1(u)即可。 时间积分不能通过该方法实现,我们需要借助于常微分和微分代数方程下的全局常微分方程和微分代数方程来实现,设置方法如下图所示,实际上将积分方程转化成了微分方程来处理。定义了一个变量w,设置其导数(wt)等于U。 点击计算,下图绘制了热能累积随时间变化的情况。 在上述算例中,全局方程在整个计算域内都是有效的,但是有些情况下,我们需要变量考虑某些点上的行为。比如假设我们的模型是一根木条,在不停加热的情况下,局部区域肯定先会被烧焦,这种情况下我们就需要用到域常微分和微分代数方程了。判断是否被烧焦的依据就是该点的能量累积超过某个阈值,我们定义一个参数P用来表示域内各点的能量累积,并假定P超过20时,该点被烤焦。 积分方程同样被转化为微分方程求解,下图中的二阶导系数置0,一阶导系数为1,求导结果等于u。 计算即可得到P随时间的变化情况,我们区别绘制了P>20和P<20的区域,可见,随着时间的推移,木条很快被烤焦。 本次我们列举了COMSOL中常见的偏微分方程和常微分方程,并通过薄膜振动和固体传热两个案例说明了部分方程的用法。基于方程的建模形式不仅为我们提供了物理场接口之外的建模方案,对于一些简单的微分方程,我们也可以方便地通过COMSOL对其求解,进而研究其行为,促进我们对微分方程的认识。 |
|