COMSOL-Multiphysics自定义PDE

仿真智领创新Simulatinginspiresinnovation中仿科技2014年1月COMSOL自定义方程方法技术部仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation主要内容自定义PDE系数型边界约束通式型边界失效弱解型弱边界ODEs和DAEs仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation系数型方程注意：时刻记住什么是边界条件，什么是初始条件！1.线性：系数只是空间或时间的函数或常量2.非线性：系数是因变量的函数仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation方程中各项含义质量阻尼质量扩散对流源对流吸收源fauuuuctudtueaa)(22守恒通量仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation质量阻尼质量弹性力初始/热应力惯性力(重力)aaedccu22()aaedcafttuuuuuu系数形式——波动方程密度阻尼系数应力刚性，“弹簧常数”仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovationfauuuuctudtueaa)(22堆积/储存扩散对流源对流吸收源系数形式——输送扩散方程仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation扩散Helmholtz项源22()2cukufakkHelmholtz方程：波数波长系数形式——频率响应波动方程fauuuuctudtueaa)(22仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation案例1——Telegraph方程式•电报线按照Telegraph方程传递电压脉冲。仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovationα+β=0.5时结果仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovationα+β=0.5时波形传播仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation有限元方法•将偏微分方程转化为离散形式的线性方程组•自动检测–非线性(线性求解或者非线性求解)–对称性(节约内存)•非线性问题：未知场或者场的导数为•材料系数•约束条件•矩阵K称为刚度矩阵仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation线性求解•通过线性变换转换矩阵的形式，即u=K-1F–将Ku=F转换为LUu=F，简便，稳定–u=U-1L-1F•优点:非常稳定•缺点:占用内存大仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation刚度矩阵的对称性•对于有限元问题总可得到对称矩阵K•刚度矩阵沿对角线对称，即：K=KT•降低内存占用量仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation刚度矩阵的正定性、对称性•通常的有限元法问题总有正定的、对称刚度矩阵K(所有特征值为正)•下列问题的解等价于最小势能：–线性结构(弹性势能)–热传导(热能)–静电学(电学)–扩散现象(质量势能)•这些问题均适合用有限元方法求解仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation迭代求解•对于大型问题•矩阵的变换占用大量内存–L和U比刚度矩阵有更多的非零项•如何避免矩阵变换？•迭代求解：–没有矩阵变换过程–在每次迭代计算后判定r=Ku-F~0–r为余量•优点:节省内存空间•缺点:稳定性较差，计算速度慢仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation前处理器•ConjugateGradient,F/GMRES,BiCGStab,GeometricMultigrid–本次课程不详细讲解这些算法•这些方法均使用前处理器•平衡方程Ku=F的两端乘以预调节器矩阵M，以增加条件数Ku=FMKu=MF仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation求解器•直接法–数值稳定性与矩阵条件数K相关性不大–COMSOL中的直接求解方法有:•PARDISO:速度快，适合于共享内存计算•SPOOLES:速度慢，占用内存最小(离散式计算)•MUMPS:速度慢，稳定，占用内存最多，适合于集群计算•迭代法：–比直接法计算速度慢–收敛性决定于矩阵的条件数K–不同的物理场需要不同的求解器配置(缺点)•发展中的研究课题•计算方法仍然在发展完善仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation线性问题•线性问题的特点1.载荷为0则解亦为02.解与载荷成比例3.解的正负随载荷的正负变化载荷解载荷解线性非线性仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation非线性问题•材料属性是u的函数，k(u)•载荷是u的函数，b(u)uk,b(u)线性uk,b(u)弱非线性uk,b(u)强非线性uk,b(u)高度非线性三种情形之间没有明显的界限19仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation余量无解余量多个解解解并非所有非线性问题都有唯一解仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation有限元方法示例k=4N/muF=2N作用点外力平衡:f(u)usolution=0.5uk=euN/muF=2N作用点外力平衡:uf(u)usolution≈0.853仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation求解线性问题f(u)=2-4u=0f(u)usolution=0.5方法：1.初始值：u0=02.斜率：f’(u0)3.求解：(1)(2)(3)对于线性问题，初始值可任意选取牛顿拉弗森迭代法：仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation这里，有限元法与牛顿拉弗森迭代法等价仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation用同样的方法求解非线性问题计算方法：1.初始值：u0=02.线性点u0=03.牛顿拉弗森迭代法4.精度设置：f(u)=2-exp(u)uuf(u)24仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovationf(u)=b(u)-K(u)uf′(u)=S(u)迭代精度:p=2k2=exp(u2-u1)k3=4u1u2u3k1=6+(u3-u1)2S(u)为雅可比矩阵用同样的方法得到非线性问题的矩阵仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation非线性问题都具有收敛半径f(u)=2-exp(u)uuf(u)若初始值在线左侧，则解不收敛u0=-1非线性问题的收敛性决定于初始条件仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovationf(u)=2-exp(u)uuf(u)缓慢收敛快速收敛收敛速度的快慢没有严格的界限。应根据实际经验选择合适的初始值。选择合适的初始值可提高收敛速度27仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovationuf(u)choose(α1)&(ααprevious)while(|fi+1||fi|)recomputefi+1(ui+1)end|fi|无阻尼速度|fi+1|中间值|fi+1|最终值|fi+1|f(u)=2-exp(u)u引入阻尼系数α提高收敛速度28仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation绝对容差与相对容差相对容差内的解绝对容差内的解tutu精确解tu0附近的解相对容差不可用，因为这会导致非常小的时间步长tu限制绝对容差与相对容差，直到结果满意为止仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation向后差分和广义α•收敛判断准则•BDF•广义alpha仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation稳态问题——系数c0022tutu稳态：仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation各项异性时的c仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation系数α和γ仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation系数β和a，f三个量都是矢量,,仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation瞬态问题——系数ea和da仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation系数的时间微分设为γ项设为f项有空间分布的时间微分，可以利用γ，f项添加：仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation特征值问题1.对比瞬态求解，将时间微分项替换为lambd2.处理线性特征值，求解非线性特征值要结合ODE约束处理仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatinginspiresinnovation边界条件1通量/源边界零通量少了一项？约束狄氏边界周期边界仿真智领创新Simulatinginspiresinnovation仿真智领创新Simulatingin