[转帖]怎么在femlab里写入自己得方程

cora

femlab的自定义方程主要有三种形式：
1、参数形式（coefficient form）
2、普通形式（general form）
3、弱解形式（weak form）
其中参数形式和普通形式十分类似，参数形式主要解决线性问题，而普通形式可以解决线性和弱的线性问题。而弱解形式主要用于解决非线性问题以及一些无法用前两者表达的线性问题。弱解形式是功能最为强大的一种求解方法，前两者可以解决的问题都可以用弱解形式解决，只是要费一些脑筋。
既然是初级讲座就只讲前两者吧，因为弱解形式较为复杂，需要有有限元解法的一些基础知识，而且使用上需要用到分步积分法。

一、参数形式

1、稳态问题
稳态问题是指求解量不随时间变化。对于静力学，求解量确实是与时间无关的静态量，例如位移、应力、应变等。而对于动力学、波动力学、电磁学等问题，求解量不随时间变化只是指其幅值与时间无关。
femlab指定的参数形式的方程如下（图1）：


其中第一个式子是求解域上的方程。第二、三个式子是边界条件，分别为第三类边界条件和第一类边界条件，即罗宾边界和狄利克边界条件，两者只能选取一个。
注：罗宾边界条件在femlab里常常被成为（广义的）牛曼边界条件。
方程里各参量的意义，在不同的问题里分别有不同的物理意义，但由于内在数学形式的一致性，因此其意义也十分类似，现以连续介质力学即扩散力学为例进行说明。
c－扩散系数
a－吸收系数|
f－源项
α－保守通量对流系数
β－对流系数
γ－保守通量源项
q－边界上的吸收系数
g－边界上的源项
这里需要注意的是源项f，对于弹性力学，对应于加在弹性体上的荷载。对于许多的物理场(声场、电磁场、温度场)问题，对应于场源（声源、点电荷、热源）。
而其他系数大都用于定义材料和媒质的特性。其中α、β可以是矢量，c可以是矩阵，用于表示各项异性的材料特性。
现在给出一个简单地判断问题是线性还是非线性的方法，若这些参量是求解变量u的函数，则问题是非线性的，反之为线性。当然还有个别的问题不符合这条准则，但要准确判断一个问题是否线性是比较困难的，详细说明可能要花上一章的功夫。呵呵，不说了，大家就用这条简单的准则判断吧。
femlab指定的PDE的参数形式已经包含了对空间坐标的二阶、一阶、0阶导数，因此大部分常见的线性PDE都可以转换成femlab指定的参数形式。其方法是将求解变量对空间坐标的导数按照阶数的高低依次排列，然后对照写出相应的参量。现举一例如下：
例：稳态下的声学波动方程如下（图2）：

cora

其中ρ为媒质密度，cs为声速，ω为声波角频率，p为声压。因此对比femlab的参数形式方程可以得出相应的参量的值为：c=-1/ρ,a=ω²/ρ/cs²,其余的参量为0。
具体操作步骤如下：
（1）在Model Navigator 窗口里选择PDE modes-PDE coefficient form-stationary analysis
设置求解变量的名称为p，如下图3：

（2）画出求解域图形，这里我用的是长为2，宽为1.2的矩形，中心在原点。选择physics－subdomain settings，在c一栏出填入-1/rho；在a一栏处填入omega^2/rho/cs^2，如下图4：

（3）在options－constants处根据实际情况分别填入rho、omega、cs的值，如图5：

（4）最后可以设置边界条件求解了。
假设边界条件是声学硬边界，即∂p/∂n=0。此为牛曼边界，且g=q=0。在physics－boundary setting中选择要设置的边界，并选择Neumann boundary condition，设置q、g都为0。如图6：

下表列出了各边界的边界条件，大家可以根据步骤（4）的方法依次设置各边界。*
边界 边界条件 类型 设置项仿
1 ∂p/∂n＝0.05 neumann q=0,g=0.05?
2,3 ∂p/∂n＝0 neumann q=0,g=0
4 ∂p/∂n＝i*ω*p/100 neumann q=-i*omega/100,g=0
设置好后就可以网格化，求解了。p的幅值求解结果如图7：

cora

（5）现在给一个复杂一点的情况。
其实方程和边界条件里的各参量都可以是空间坐标的函数，对于本例，假设媒质的密度是x的函数，即ρ＝ρ（x），原方程变为（图8）：

这里假设ρ（x）＝1.2*abs（x）
我们增加一个名为rhox的变量用以代替原来的rho常量。
在option-expressions-scalar expressions中加入一个名为rhox、表达式为rho*abs（x）的变量。如图9：


然后在subdomain settings的页面中，出现rho的地方全部用rhox代替。
边界条件中各参量也可以是空间坐标的函数，假设边界1的条件变为：
∂p/∂n＝0.05*y，则只要在boundary setting中将 g一项设置为0.05*y即可。更改后的求解结果如图10，注意图中显示的是幅值，即abs（p）。

cora

2、参数形式的动态问题。
所谓动态问题是指求解量与时间相关，反映在方程中就是含有时间t的项，通常是关于时间的一阶或二阶的导数。如图11：

比较稳态问题和动态问题的方程，可以发现后者只是比前者多了一个d_a*∂u/∂t 一项，其中的d_a被称为质量系数，这是因为这个系数通常都要被研究物体的质量或密度有关。
当然，很多的动态问题都含有对时间t的二阶导数的项（例如波动方程），这类动态问题的方程可表示为图12所示：

其实含有时间二阶导数项的方程，可以适当转换为方程组，使每个方程都只含有时间一阶导数项。例如图12所示的方程就可以通过引入一个中间变量v，从而建立一个只含时间一阶导数项的方程组，如图13：

建立方程组可在model navigator窗口的dependent viabless一栏中分别填入方程组的所有变量名，如图14。

注意填写方程组的参量时，各参量都应该是矩阵或数列的形式，会比较复杂。当然femlab已经内置了含时间二阶导数项方程的模式，称为"time-dependent analysis,wave extension"，大家可以直接使用和以之为参考。
当然，我们也可以用弱解形式来求解这类问题，在这里我就不谈了。

cora

3、特征值问题
所谓特征值问题，其方程的形式如图15。

与稳态问题的方程相比较，多了 一项λ*d_a*u。其中λ就是特征值，由特征值我们还可以算出特征频率。
求解特征值问题可以得到各特征值的值，以及不同特征值所对应的变量u的值。
二、普通形式
普通形式和参数形式十分类似，也同样分为稳态、动态、特征值问题。为简练起见，以下仅以稳态问题为例。
普通形式的稳态问题的方程如下：
图16：

比较参数形式和普通形式的方程，可以发现存在以下关系：
Γ=-c▽u-α*u+γ
F=f-β▽u-a*u
G=g-q*u
R=r-h*u
因此，参数形式的方程都可以通过以上公式转换为普通形式的方程。
例如对于前面的第一个例子，因为c=-1/ρ,a=ω²/ρ/cs²，所以Γ=▽u/ρ，F=-ω²*u/ρ/cs²在femlab中的设置如图17。

边界条件同样可以根据以上公式进行转换，边界1～4的设置见下表：
边界 类型 g
1 Neumann 0.05w
2 Neumann 0
3 Neumann 0
4 Neumann i*omega/100*u
图18为边界4的设置，其他边界类似，略。
图18：

然后是划分网格和求解，求解结果与图7完全一致。
最后可能有人问，既然参数形式可以转化为普通形式，那么普通形式的存在意义是什么呢？其原因是，参数形式适用于求解线性问题，而普通形式可以求解线性和一些非线性的问题。还有更重要的是，很多问题容易用普通形式表达，却难以使用参数形式表达。
例如对于非线性方程：∂²u/∂x²=u*∂u/∂x，我们可使用普通形式令Γ=-∂u/∂x，F=-u*∂u/∂x即可。但由于方程是非线性的，因此难以使用参数形式表示。

abror

　　

谢谢，受益匪浅!!

abror

？


你的图11和图12好像是一个图啊？

TNC

可能是传错了吧

chbo76

我这里怎么看不到图啊？

knights

看不到图呀