怎么在femlab里写入自己的方程

zaza8202

femab介绍里有这么一条
'4、FEMLAB 在科研方面：
定义和耦合任意数量偏微分方程的能力使得FEMLAB成为一个强大的分析工具。其灵活性和基于方程的建模方式可以帮助用户深入在MEMS、纳米技术、燃烧室、光子学、生物工程和许多其它领域内的研究。'
请问各位大侠怎么把自己的方程输入呢?

[ 本帖最后由 zaza8202 于 2006-11-16 16:46 编辑 ]

临江仙

将自己的方程写成coefficient form或general form指定的形式，然后在subdomain设置中填入相应的项。

zaza8202

请问有没有办法在matlab里把自己的方程写入呢.

yycgy

临江仙，

     建议你搞个初步讲座，教会大家怎么在femlab中输入自己的方程，谢谢！

临江仙

既然版主要求，有时间的话我一定会写上的。

至于matlab中能不能将方程写入，我不是太清楚。因为matlab的偏微分工具箱比较弱，我一般都不使用。不过好像matlab的方程形式跟femlab的coeffient form的方程形式类似，所以如果可以写成femlab的coeffient form形式，稍加变换应该可以写成matlab的形式。
不过如果方程的形式比较复杂，例如有非线性项的话，可能matlab无法解决。

yycgy

临江仙，
     你说的不错，对于非线性问题，matlab只能通过编程处理了（一般采用有限差分）。希望早日看到你的ppt！

zaza8202

强烈期待 临江仙 的讲座

临江仙

今天先讲一部分，若有错误请各位老师同学指正，欢迎补充。

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

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

图1：

临江仙

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

临江仙

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

临江仙

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

临江仙

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

临江仙

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

zaza8202

呵呵，强烈建议给临江仙同志加积分

临江仙

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

临江仙

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

临江仙

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

renhy115

建议斑竹给临江仙加分，我觉得如此热心公益的同志应该鼓励一下啊！呵呵

yycgy

已经给临江仙加分，希望老兄多多发好贴！