SimWe仿真论坛's Archiver

zaza8202 发表于 2005-8-2 10:26

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

femab介绍里有这么一条
3M8y)Y[tX '4、FEMLAB 在科研方面:
H!V(X&d'hf7[R9y 定义和耦合任意数量偏微分方程的能力使得FEMLAB成为一个强大的分析工具。其灵活性和基于方程的建模方式可以帮助用户深入在MEMS、纳米技术、燃烧室、光子学、生物工程和许多其它领域内的研究。' cC ZT)]3~q4{ y
请问各位大侠怎么把自己的方程输入呢?
2|NCD m ~){Z E YV/{a
[[i] 本帖最后由 zaza8202 于 2006-11-16 16:46 编辑 [/i]]

临江仙 发表于 2005-8-3 00:18

Re:怎么在femlab里写入自己得方城

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

zaza8202 发表于 2005-8-6 21:40

Re:怎么在femlab里写入自己得方城

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

yycgy 发表于 2005-8-8 10:55

Re:怎么在femlab里写入自己得方城

临江仙,
)O8t }A}\
eOW3P+z!HLA      建议你搞个初步讲座,教会大家怎么在femlab中输入自己的方程,谢谢!

临江仙 发表于 2005-8-8 16:41

Re:怎么在femlab里写入自己得方城

既然版主要求,有时间的话我一定会写上的。
`6{:[A+S0o ]._ajM6w g7fg
至于matlab中能不能将方程写入,我不是太清楚。因为matlab的偏微分工具箱比较弱,我一般都不使用。不过好像matlab的方程形式跟femlab的coeffient form的方程形式类似,所以如果可以写成femlab的coeffient form形式,稍加变换应该可以写成matlab的形式。Qyb!x x(PeX{*i H
不过如果方程的形式比较复杂,例如有非线性项的话,可能matlab无法解决。

yycgy 发表于 2005-8-9 09:01

Re:怎么在femlab里写入自己得方城

临江仙, xo\y'm@
     你说的不错,对于非线性问题,matlab只能通过编程处理了(一般采用有限差分)。希望早日看到你的ppt!

zaza8202 发表于 2005-8-10 10:03

Re:怎么在femlab里写入自己得方城

强烈期待 临江仙 的讲座

临江仙 发表于 2005-8-12 22:45

Re:怎么在femlab里写入自己得方城

[b]今天先讲一部分,若有错误请各位老师同学指正,欢迎补充。[/b]
x)N6|kr)I5EA!h!X 1i$^jq&d
femlab的自定义方程主要有三种形式:$l/O7uO&x8F
1、参数形式(coefficient form)(H"FA Z7L&Qt&m(n^{
2、普通形式(general form)
zF%~e{ 3、弱解形式(weak form)
&{A&l#EF3f&z     其中参数形式和普通形式十分类似,参数形式主要解决线性问题,而普通形式可以解决线性和弱的线性问题。而弱解形式主要用于解决非线性问题以及一些无法用前两者表达的线性问题。弱解形式是功能最为强大的一种求解方法,前两者可以解决的问题都可以用弱解形式解决,只是要费一些脑筋。
p7JF;T aF{     既然是初级讲座就只讲前两者吧,因为弱解形式较为复杂,需要有有限元解法的一些基础知识,而且使用上需要用到分步积分法。
q"ve!d }1? ~"bu_ _3e D5qYq8H
一、参数形式
c,r(n|+| 1、稳态问题5MTj~w
    稳态问题是指求解量不随时间变化。对于静力学,求解量确实是与时间无关的静态量,例如位移、应力、应变等。而对于动力学、波动力学、电磁学等问题,求解量不随时间变化只是指其幅值与时间无关。
Yd|7sOnhA_+o6b     femlab指定的参数形式的方程如下(图1):
3x/}L%{[d)v 1k"{?*K!|vC}G p
图1:

临江仙 发表于 2005-8-12 22:47

Re:怎么在femlab里写入自己得方城

其中第一个式子是求解域上的方程。第二、三个式子是边界条件,分别为第三类边界条件和第一类边界条件,即罗宾边界和狄利克边界条件,两者只能选取一个。)rNp9[:zs3B
注:罗宾边界条件在femlab里常常被成为(广义的)牛曼边界条件。
JhA |0Y"lna } 方程里各参量的意义,在不同的问题里分别有不同的物理意义,但由于内在数学形式的一致性,因此其意义也十分类似,现以连续介质力学即扩散力学为例进行说明。
w1V ]L'HZ@ c-扩散系数"YCb(|O2thz Y B
a-吸收系数9@M#zH O}
f-源项
H McA @@ α-保守通量对流系数l`I+W9?VV
β-对流系数8t*j#mwl
γ-保守通量源项/fd yL`:?M&Y?
q-边界上的吸收系数
4Y|;rk-Pz g-边界上的源项1]#R [/R}
这里需要注意的是源项f,对于弹性力学,对应于加在弹性体上的荷载。对于许多的物理场(声场、电磁场、温度场)问题,对应于场源(声源、点电荷、热源)。 @'a#k_eN[n^4I5L
而其他系数大都用于定义材料和媒质的特性。其中α、β可以是矢量,c可以是矩阵,用于表示各项异性的材料特性。2\}8P M-b
现在给出一个简单地判断问题是线性还是非线性的方法,若这些参量是求解变量u的函数,则问题是非线性的,反之为线性。当然还有个别的问题不符合这条准则,但要准确判断一个问题是否线性是比较困难的,详细说明可能要花上一章的功夫。呵呵,不说了,大家就用这条简单的准则判断吧。*C1m#k'e,c X.@g
femlab指定的PDE的参数形式已经包含了对空间坐标的二阶、一阶、0阶导数,因此大部分常见的线性PDE都可以转换成femlab指定的参数形式。其方法是将求解变量对空间坐标的导数按照阶数的高低依次排列,然后对照写出相应的参量。现举一例如下:
'Q|iLq:z8k@{6M 例:稳态下的声学波动方程如下(图2):
3LhV%e qU 图2:

临江仙 发表于 2005-8-12 23:12

Re:怎么在femlab里写入自己得方城

其中ρ为媒质密度,cs为声速,ω为声波角频率,p为声压。因此对比femlab的参数形式方程可以得出相应的参量的值为:
+M,IIZK{/G"^ e c=-1/ρ,a=ω[sup]2[/sup]/ρ/cs[sup]2[/sup],其余的参量为0。Yd.L `7G'd?4C6p$ELi$r
具体操作步骤如下:
/T9@,Oy!p4S m!H (1)在Model Navigator 窗口里选择PDE modes-PDE coefficient form-stationary analysis
:o!\AD"WGo 设置求解变量的名称为p,如下图3:#H]@qMZ9}bF%y
图3:

临江仙 发表于 2005-8-12 23:21

Re:怎么在femlab里写入自己得方城

(2)画出求解域图形,这里我用的是长为2,宽为1.2的矩形,中心在原点。选择physics-subdomain settings,在c一栏出填入-1/rho;在a一栏处填入omega^2/rho/cs^2,如下图4:$] SVw*Zs'~+X ^k{V
图4:

临江仙 发表于 2005-8-12 23:24

Re:怎么在femlab里写入自己得方城

(3)在options-constants处根据实际情况分别填入rho、omega、cs的值,如图5:
'jtkE1zsa 图5:

临江仙 发表于 2005-8-12 23:25

Re:怎么在femlab里写入自己得方城

(4)最后可以设置边界条件求解了。
uHbU2~CJI)W 假设边界条件是声学硬边界,即∂p/∂n=0。此为牛曼边界,且g=q=0。在physics-boundary setting中选择要设置的边界,并选择Neumann boundary condition,设置q、g都为0。如图6:&x3Ur'U\/[
图6:

临江仙 发表于 2005-8-13 00:47

Re:怎么在femlab里写入自己得方城

下表列出了各边界的边界条件,大家可以根据步骤(4)的方法依次设置各边界。
_e(srR!w4l[Y 边界  边界条件  类型  设置项jh2w1nPqy!WF
1     ∂p/∂n=0.05    neumann   q=0,g=0.05
(Kb6Oii_k 2,3  ∂p/∂n=0    neumann  q=0,g=0W%R,m)@e U;c
4     ∂p/∂n=i*ω*p/100  neumann  q=-i*omega/100,g=0
{d1`pyY&`D
H[ kI!Ar}r-_)K 设置好后就可以网格化,求解了。p的幅值求解结果如图7:i]/O BaB1p8z
图7

zaza8202 发表于 2005-8-13 15:49

Re:怎么在femlab里写入自己得方城

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

临江仙 发表于 2005-8-14 00:53

Re:怎么在femlab里写入自己得方城

(5)现在给一个复杂一点的情况。
N8a`_&_V 其实方程和边界条件里的各参量都可以是空间坐标的函数,对于本例,假设媒质的密度是x的函数,即ρ=ρ(x),原方程变为(图8):
2We%}b9mP!A1Y 图8:

临江仙 发表于 2005-8-14 01:04

Re:怎么在femlab里写入自己得方城

这里假设ρ(x)=1.2*abs(x)
f2r e XQi 我们增加一个名为rhox的变量用以代替原来的rho常量。
#t j4^z1S+Tj1a 在option-expressions-scalar expressions中加入一个名为rhox、表达式为rho*abs(x)的变量。如图9:
/e T:bmB*rD"LW 图9:

临江仙 发表于 2005-8-14 01:09

Re:怎么在femlab里写入自己得方城

然后在subdomain settings的页面中,出现rho的地方全部用rhox代替。'Dif.j?f^p-k
边界条件中各参量也可以是空间坐标的函数,假设边界1的条件变为:x _r9e%D
∂p/∂n=0.05*y,则只要在boundary setting中将 g一项设置为0.05*y即可。更改后的求解结果如图10,注意图中显示的是幅值,即abs(p)。
&j&Q Z9i$v2H 图10:

renhy115 发表于 2005-8-14 09:38

Re:怎么在femlab里写入自己得方城

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

yycgy 发表于 2005-8-22 10:42

Re:怎么在femlab里写入自己得方城

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

optoelec 发表于 2005-8-23 21:20

Re:怎么在femlab里写入自己得方城

强烈支持临江仙英雄!

临江仙 发表于 2005-8-24 21:38

Re:怎么在femlab里写入自己得方城

多谢版主,不过现在我放假回家了,没多少机会用机。开学后将此讲座写完。

yycgy 发表于 2005-8-24 21:49

Re:怎么在femlab里写入自己得方城

谢谢临江仙,希望你多结合实例进行讲解。

贾济榕 发表于 2005-8-30 10:21

Re:怎么在femlab里写入自己得方城

看看你所研究的问题属于那种情况!线性或非线性,稳态或非稳态!然后根据情况,将自己的方程写成coefficient form或general form指定的形式,然后在subdomain settings 、boundary settings设置相应的参数和数值。 Nur^ N!a%V

]1^/A.y h/u6d(l] 你最好将方程传上来!大家相互学习,共同进步!!!!!!!!!!!!!!!!!!!

贾济榕 发表于 2005-9-6 09:45

Re:怎么在femlab里写入自己得方城

希望临江仙讲解一下"general form和weak form"的求解过程!

临江仙 发表于 2005-9-8 21:18

Re:怎么在femlab里写入自己得方城

今天抽了点时间把剩余的部分写完。5k/OFCH:c"m
2、参数形式的动态问题。
1Oq\ b`SB 所谓动态问题是指求解量与时间相关,反映在方程中就是含有时间t的项,通常是关于时间的一阶或二阶的导数。如图11:5KvoHjd j,s
图11

临江仙 发表于 2005-9-8 21:29

Re:怎么在femlab里写入自己得方城

比较稳态问题和动态问题的方程,可以发现后者只是比前者多了一个d[sub]a[/sub]*∂u/∂t 一项,其中的d[sub]a[/sub]被称为质量系数,这是因为这个系数通常都要被研究物体的质量或密度有关。
&c o1ivur*km 当然,很多的动态问题都含有对时间t的二阶导数的项(例如波动方程),这类动态问题的方程可表示为图12所示:
s/qGWE.y 图12:

临江仙 发表于 2005-9-8 21:38

Re:怎么在femlab里写入自己得方城

其实含有时间二阶导数项的方程,可以适当转换为方程组,使每个方程都只含有时间一阶导数项。例如图12所示的方程就可以通过引入一个中间变量v,从而建立一个只含时间一阶导数项的方程组,如图13:
VK;n"N)H'|4R[ 图13:

临江仙 发表于 2005-9-8 21:56

Re:怎么在femlab里写入自己得方城

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

临江仙 发表于 2005-9-8 22:14

Re:怎么在femlab里写入自己得方城

注意填写方程组的参量时,各参量都应该是矩阵或数列的形式,会比较复杂。当然femlab已经内置了含时间二阶导数项方程的模式,称为"time-dependent analysis,wave extension",大家可以直接使用和以之为参考。"E.Z1im+A2y/zg#WzE
当然,我们也可以用弱解形式来求解这类问题,在这里我就不谈了。8y9@Qa9[(SapE5R
%\'x+By0cE6Py&s
3、特征值问题
j1W Y9AAq+Q 所谓特征值问题,其方程的形式如图15。
~ Dy,g]$sU#t,^ S 图15:

临江仙 发表于 2005-9-8 22:32

Re:怎么在femlab里写入自己得方城

与稳态问题的方程相比较,多了 一项λ*d[sub]a[/sub]*u。其中λ就是特征值,由特征值我们还可以算出特征频率。
5D?*X6|M&@|-UVE 求解特征值问题可以得到各特征值的值,以及不同特征值所对应的变量u的值。.]6qf`V2c0yHI
u8{\#m[#K
二、普通形式 |$?nM5z/t
普通形式和参数形式十分类似,也同样分为稳态、动态、特征值问题。为简练起见,以下仅以稳态问题为例。 n p/V Z9G2hFt
普通形式的稳态问题的方程如下:%f RZ m M-Y1T
图16:

临江仙 发表于 2005-9-9 01:31

Re:怎么在femlab里写入自己得方城

比较参数形式和普通形式的方程,可以发现存在以下关系:
Zor6i_$}o b Γ=-c▽u-α*u+γd+Y:j0D/g-} g9I2]L
F=f-β▽u-a*u Bh1Ug#m&F4vq
G=g-q*u
#KW/u!W)agxx7B G R=r-h*uAM-LDS9Y*f9xs"d
因此,参数形式的方程都可以通过以上公式转换为普通形式的方程。
wG"ZT#Q#W.X6_ 例如对于前面的第一个例子,因为c=-1/ρ,a=ω[sup]2[/sup]/ρ/cs[sup]2[/sup],所以Γ=▽u/ρ,F=-ω[sup]2[/sup]*u/ρ/cs[sup]2[/sup]t-]4} Q*Ao3l
在femlab中的设置如图17。d rW7}j@@/{(]
图17:

临江仙 发表于 2005-9-9 01:37

Re:怎么在femlab里写入自己得方城

边界条件同样可以根据以上公式进行转换,边界1~4的设置见下表:
8e._ Pi)LUn 边界  类型  g  U/cIr$H#?S(jrV
1 Neumann 0.05
+G4V~._1^w 2 Neumann 0r |G.R2kYT|
3 Neumann 0
3s'B pd&vX!p"f6\`p7d 4 Neumann i*omega/100*u(J7q$M~T3[IR g*e
图18为边界4的设置,其他边界类似,略。d wi:o s4T!SG/p
图18:

临江仙 发表于 2005-9-9 01:55

Re:怎么在femlab里写入自己得方城

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

临江仙 发表于 2005-9-9 01:58

Re:怎么在femlab里写入自己得方城

这个初级讲座就完成了。大家有什么心得欢迎补充。Yu8r W:o!]
h!gxzvQ
[color=red]讲座很精彩,谢谢老兄,已经加分鼓励 了。希望能看到老兄的中级和高级讲座,偶会大力支持你的,希望结合工程实例进行说明,方便大家交流。[/color]

zaza8202 发表于 2005-9-9 09:45

Re:怎么在femlab里写入自己得方城

强烈建议再给临江仙 加两分
1p0f3H+zY0L_~ 斑竹我提出这么好的问题是不是奖励一分呀^_^
(\;bdffV"_;O.V
"NnGw$N#s)A [color=red]奖励该奖励的,希望共同努力,多出一些好贴。[/color]

dqakl 发表于 2005-9-13 14:37

Re:怎么在femlab里写入自己得方城

刚入门femlab,不懂的实在太多!慢慢学习!!

dqakl 发表于 2005-9-13 14:38

Re:怎么在femlab里写入自己得方城

斑竹能不能支持一下讲解一些温度场随时间变化的实际例子呀?谢谢啦!

horfman 发表于 2005-9-19 09:58

Re:怎么在femlab里写入自己得方城

大家好,。欢迎进入FEMLAB论坛电磁专场,QQ群号码:13546386

jdlu 发表于 2005-9-26 20:58

Re:怎么在femlab里写入自己得方城

临江仙,下面的讲座什么时候开始阿,盼望啊…………

页: [1] 2 3
 

Powered by Discuz! Archiver 6.1.0  © 2001-2007 Comsenz Inc.