有朋友反映咱们杂货店里推送的高水平SCI论文太过物理，对于一些工科或者其他专业人来说可能有些障碍。模数哥以后也会经常写一写具体的Comsol操作教程，当然也不会像案例库里那样具体“手把手”，而会强调对软件算法或者底层数学规律的介绍。

今天的话题是如何对变量(variable)进行空间微分或者时间积分。可能朋友们会比较诧异：这有什么难的呢？确实在Comsol软件里有具体的函数或者算法，可以帮助我们进行这些操作，如d(f,x)即是对变量f沿左边x方向偏微分，timeint (t1, t2, expr)是对表达式expr在时间t1到t2内进行积分，但是有些情况下，这些操作符并不能直接解决我们的问题。福利1：百度“常用comsol的操作符和数学函数”，即可很容易找到相关更多资料。

问题1：变量的空间微分

问题描述：我们拿一个案例库中ACDC模块的helmholtz coil为例，请问各位能否画出mf.Hx沿x方向的微分(即Hx在x方向的梯度)呢？当我们在后处理中尝试画出d(mf.Hx, x)时，Opps，结果为零，这明显是有背于真实数学和物理的！福利2：可以百度一下“Plotting Spatial Derivatives of the Magnetic Field”找到Comsol官网的相关教程。该教程介绍了如何通过新建耦合偏微分方程(Partial differential equation, PDE)来对mf.Hx变量进行空间求导，但是里面语焉不详，并未介绍一些底层的数学原因。

要搞清楚这里的原因，我们首先要知道有限元算法中形函数(shape funcion)的概念。在用网格格点对整个仿真系统进行离散化时discretization，每个格点都是一个未知的方程，而这个方程的表达形式，就可以简单理解为形函数(其实并不严谨)。比如一般在力学模块里，常用的形函数是朗格朗日二次形函数(Quadratic Lagrange)，即每一个格点的数学展开形式是二次多项式。如下图所示，对于空间任意一点x的数值f，我们均可以通过其周围几个点的数据进行插值interpolation得到。

如果我们要求解空间任意点x处变量f的偏微分，数学上其实是先对周围各点的形函数进行偏微分，再做插值处理得到x点的空间导数。以x1点为例，它的一阶微分是f1x=2a1*x1+b1，f1x中的x表示f1的x偏微分。它的二阶微分是f1xx = 2a1，很明显，它的三阶微分就是0了。所以这种有限元求微分的方法决定了朗格朗日二阶形函数下x点的三阶微分就是0——物理和数学上都是错误的。

那么怎么样才能求解变量f的三阶微分呢？方法就是我们可以新建一个PDE方程并将f1xx耦合为该PDE的主变量u，一般地，该PDE方程也默认是用二阶朗格朗日形函数表示，如下图，所以我们便可以再次对u(即f1xx)进行空间求导了。所以理论上，利用这种不断耦合到新PDE方程的形式，我们可以实现对f变量的任意阶空间求导或微分。

现在让我们再回到前文所述helmholtz coil案例中，该模型用的接口是3D的磁场Magnetic Field (mf)，该物理场的变量是磁矢势Magnetic Vector Potential A，3D情况下A是一个矢量，而磁感应强度便是该矢量的旋度：▽×A = B，对矢量求旋度的数学处理就等效为二次偏微分了，所以我们无法再对磁感应强度B进行空间微分。这也是我们发现d(mf.Hx, x) =0的原因。

具体如何新建并耦合系数型PDE方程来对磁场H求其空间梯度，朋友们可以参考上文所述Comsol官方blog：

结果如下，磁场的梯度u1x(即d(mf.Hx, x))主要出现在线圈的周围：

我这边重点介绍一下这个系数型偏微分方程。如下图所示，该方程其实是一个万金油式的功能型方程，几乎可以包含宏观尺寸各类物理问题：

我们所享用的Comsol中现成的各个物理场接口，其实也大都是该方程的变种(各种情况的简化)罢了，比如稀物质扩散(tds)或者热传输场(heat transfer)，都是用的下图所示扩散方程：

所以Comsol软件的前身是Matlab求解方程的一个工具箱FEMLab，后来独立更名为Comsol并一直沿用至今，所以Comsol的二次开发就是采用Matlab编程，这也决定了Comsol仿真的一些劣势(因为底层程序语言效率的不足，在非线性、流场、结构场粘弹性仿真等方面被Fluent、Abaqus碾压)。学习并掌握comsol中PDE方程的建模，设置各类边界条件和初始值，对于我们把握Comsol多物理场耦合有限元的数学思想是非常有帮助的。

问题2：变量的时间积分

问题描述：以案例库Laser Heating wafer为例，当一个固定加热光源打在一个旋转的硅片上，即等效为一个固定的硅片被一个旋转的光源所加热ht.q0。以其中某点为例，如何求出热通量的时间积分并画出随时间的曲线呢？

在后处理里画出Qtot随时间的曲线是不方便的，其实最简单的办法是新建并耦合一个常微分方程ODE，将热通量q0耦合进来进而求出其时间积分u：

若我们只想考察某点的热通量总和，我们只需要建立一个分布式的(distributed)点的常微分方程即可，其方程形式如上图所示，设源项Source term是热通量ht.q0，损耗系Damping coefficient数为1，质量系数为0，则ODE的变量u即为热通量的时间积分。所以我们便可以画出该点的温度T和热通量总和u(即Qtot)随时间的变化：

该点的温度随时间变化呈现一个个脉冲的形式，即该点在旋转过程中一次次经过热源并被加热。而总热通量随时间呈现出阶跃的变化，即是热通量的一次次累积。这样一来，物理意义十分清晰。

最后总结一下，有些朋友找到模数哥，说想在comsol里面建PDE、ODE的耦合方程，一下子就晕了，不知道各个系数怎么设置。我觉得不要太从数学的眼光来看待这些方程，最好是从物理的角度来看待他们，每一个系数都是有明确的物理意义的，这样理解起来也是最方便的。当然如果想尽快解决你手上的问题，最快捷的方法是找模数哥付费咨询了。