二维粘性不可压流域中模拟源、汇(页 1) - A16：FEPG - SimWe仿真论坛 |Simwe.com|仿真|设计|有限元|虚拟仪器

saint_suan 发表于 2006-12-22 23:24

二维粘性不可压流域中模拟源、汇

我想用飞箭模拟二维平面区域中存在源和汇时的流场分布情况。几何模型如图，x0y平面区域（10×5）中挖两个孔，用给定的垂直于孔边界的速度来模拟源和汇的强度。
看了一个飞箭作流体的例子，就是方腔流，顶盖匀速拖动的例子（单一流场，只求解N-S方程，不考虑热什么的），觉得里边的文件可以不加修改地直接用，我用了五个文件，描述速度u,v和压强p的cav.pde文件，求解流线的文件dnsulxy.vde,nellf.nfe,还有系统文件cavity.gcn,cavity.gio，如下面都给出。就是生成.pre文件后修改a场的初值为1,然后在GID中给整个平面赋初值u=5,v=2,p=100。其他边界条件也如下给出。J(F8ML6X.G?0d
但DOS批命令求解时DOS窗口都是unknow，没算出数来。工作目录下disp0，elem0，coor0文件都有，就是空的。_*_6R*V]2T a
我不知道是怎么回事，烦请哪位高人指点一下。%d0X2SQ'K C7[
几何和网格如图，结构化网格，quadratic9单元。O1R#R!yY

[b]边界条件如下：[/b] q&V7{7ZRg
上下两边:
u-I=1, u-D=0,Q)wvX T;Qfd-\
v-I=－1,v-D=0,
p-I=1,p-D=0;
左右两边：
u-I=－1, u-D=0,
v-I=1,v-D=0,
p-I=1,p-D=0;
左边的圆孔模拟源（圆心坐标（0，0，0））：
u-I=－1, u-D=x,g#c-Q \~ Q!uWS
v-I=－1,v-D=y,dejr9KvH
p-I=1,p-D=0;)bV3w:zzDx0a
右边的圆孔模拟汇（圆心坐标（5，0，0））：Y#m T5PnK(a
u-I=－1, u-D=x*(－1)+5,
v-I=－1,v-D=y*(－1),~tO'jB4p,Pd
p-I=1,p-D=0;t l Esg1q\0~;f
8WM!v-Jf/F|
[b]结构化网格如图[/b](B4?$GstX8S:x

juu ~T5hb~K
[b]cav.pde如下：[/b]
\this method can only applied to structured rectangle mesh.Ta$]o;?k!T^;j
\pde file for solving u,v p in momentum equation
\non-linear iteration excuted by newton-raphson method
\gls stabilized terms are added to account for dominative convective term
\we solve the navier-stokes equation by Galerkin\least square method4B/nK9ep8X,L
\that isy!U1hk?
\we take the virtual displacement as u+(u/x*un+u/y*vn+p/x/rou)*taus e+j1f5Ln:Mvd
\ and v+(v/x*un+v/y*vn+p/y/rou)*tausN(Z2x_K'VI1@a
\ so the eqation is (note if we take linear-function as base-function on
\ a structure mesh, we can deduce that -u/x,x=0,-v/x,x=0) Y'D0mur
\ (u/x;u/x)*amu+(u/y;u/y)*amu+(v/x;v/x)*amu+(v/y;v/y)*amu
\ +(u/x*un+u/y*vn;u)*rho+(v/x*un+v/y*vn;v)*rho-(p;div)+(div;p)N,A-}}rE
\ +(u/x*un+u/y*vn;u/x*un+u/y*vn)*rho*taus
\ +(v/x*un+v/y*vn;v/x*un+v/y*vn)*rho*taus?s,r4L+|Wq;T
\ +(u/x*un+u/y*vn;p/x)*taus+(v/x*un+v/y*vn;p/y)*tausw"@"? nO.p+}
\ +(p/x;u/x*un+u/y*vn)*taus+(p/y;v/x*un+v/y*vn)*taus0\!}p\r
\ +(p/x;p/x)*taus/rho+(p/y;p/y)*taus/rho
\= (fx;u)+(fx;(u/x*un+u/y*vn+p/x/rou)*taus)
\ +(fy;v)+(fx;(u/x*un+u/y*vn+p/x/rou)*taus).A1Uom&|9E]
disp u v p+eq!EOa3?(BO
coor x y
shap q 4{9u,M#Y t GlW&h+y
gaus q0u)DaFn];Y Z
coef un vn pn
mate amu rho fx fy 1.0d-3 1.0d0 0.0 0.0aCSqt0i)kz
func udsumx udsumy udux uduy gradpx gradpy eux euy evx evy div5CU``N$P#P^*h

func@}A.c3hS#T7_lL6a
$c6 h2=abs(coorr(2,1)-coorr(2,4))
$c6 h1=abs(coorr(1,1)-coorr(1,2));n7y'vP"?-RC
$c6 h2a=abs(coorr(2,1)-coorr(2,2))2? j){;c a df
$c6 h1a=abs(coorr(1,1)-coorr(1,4))IK"E6]/n @*CVP
$c6 if(h2.eq.0.0) then+w8sy} b+KI/n\
$c6    h2=h2a;
$c6 endif;
$c6 if(h1.eq.0.0) thenJ:c\H6e1t
$c6    h1=h1a;-a f_].[*{II1wu
$c6 endif;
$c6 unorm=sqrt(un*un+vn*vn)
$c6 if(unorm.eq.0.0) then
$c6 taps=rho*h2**2/amu/3.0!@ Qd#uiz.M
$c6 elseUY{*Ct
$c6 udni1=(un/h1+vn/h2)bQ8@J{Nj
$c6 udni2=(un/h1-vn/h2)
$c6 udni=2.0*(abs(un/h1+vn/h2)+abs(un/h1-vn/h2))D6cxFve
$c6 hugn=2.0*unorm/udni6{^Wh\4\
$c6 pe=0.5*rho*unorm*hugn/amuG^6p5G2Ay M
$c6 if(pe.lt.0.1) then ^l:`QR
$c6    fnpe=pe*4.0/3.0
$c6 elsem] B(_C'y
$c6 if(pe.lt.3.0.and.pe.gt.0.1)thenF(Q.h1S D
$c6    fnpe=pe/3.0\fAI;zV {-S&b
$c6 else/[!@7Y$nTOBo;B
$c6    fnpe=1.0IO#aW/[m
$c6 endif
$c6 endif
$c6 tasu=0.5*hugn/unorm*fnpe#N,W[j,yT G&}A#x}N
$c6 taps=tasu#NT!` ^n.d CK
$c6 tals=0.5* hugn*unorm*fnpe2O`]Jq?x4P \*i9`}
$c6 endifz#`g/Y"WXQ
udsumx=+[u/x]*un+[u/y]*vnh}]6`(c"K4c1L6Q
)n:Tj t/RN8r
udsumy=+[v/x]*un+[v/y]*vn
6n bIYp3Ov
eux=+[u/x]9K#} OA1L1n(N

euy=+[u/y]K/jV+E$`4evB
tzlVN4eBp3V
evx=+[v/x]

evy=+[v/y]9qB9i nC-V'C
3g-_5o5gg{6] Y4Q.cq
div=+[u/x]+[v/y]+x)s}3U y
"r t7Q Q9k2L*Gm
gradpx=+[p/x]

gradpy=+[p/y]

udux=+[u/x]*un+[u/y]*vnBG'dq(j`*}-a
*Zg1Lk4cI6b0r3s
uduy=+[v/x]*un+[v/y]*vn

stif w0V(YP-e#L!|'p
dist=+[eux;eux]*amu+[euy;euy]*amu+[evx;evx]*amu+[evy;evy]*amu
+[udsumx;u]*rho+[udsumy;v]*rho+[div;p]-[p;div]
+[udux;udux]*rho*tasu+[uduy;uduy]*rho*tasu
+[gradpx;udux]*tasu+[gradpy;uduy]*tasu,mZoS-wJ
+[udux;gradpx]*taps+[uduy;gradpy]*taps)N*f"Q"L#i!b/k
+[gradpx;gradpx]/rho*taps+[gradpy;gradpy]/rho*taps
+[div;div]*rho*tals
:IVp&mM?7s)~"zhR
load=+*fx*rho+[v]*fy*rho.J:}L3Y1xxd
+[udux]*rho*fx*tasu+[uduy]*rho*fy*tasu9D;_j-Zh2]
+[gradpx]*fx*taps+[gradpy]*fy*taps

end4g-{;F!ddE2Z
\,D9J|EJwR3M4?
[b]nellf.nfe文件如下defi[/b] Km'X{T3N+X;?
stif s
mass m
load f
type e
mdty lj+~ pA2i
step 0

coef u!e\#KO2z5`:m
7pugJhkD
equation
vect u T3H'^RH+OFf
read(s,unod) u;T+J Z-prR:z+M B|
matrix = [s]
forc=[f]K{BIga|
J9})zfr:NQ@
solution u
vect u1,u,ue,du
$c6 open(99,file='elem.it',form='formatted',status='old')$Cn6i;F1DmY
$c6 read(99,*) itn,cc{2yv"C&U8~o#~1p
read(s,unod) u1,duJqG(u;B,eHx
write(s,unodf) u1
write(s,unodr) u
[ue]=-[u1]
$c6 aa = 0.0$`7{h6zinf
$c6 ab = 0.0*G'a'Y(E0k#k*mY
$c6 bb = 0.0.PR.q'ak5[8l;O5Q
%nod
%dof
aa = aa+[ue]*[ue]
   ab = ab+[ue]*[du]5o!L3ZOF1Qq
   bb = bb+[du]*[du]*Ey+P@$jB:g|
%dof
%nodW-K[UH!`
$c6 if (bb.gt.1.0e-20) then3ha1S$f9I
$c6 rab = sqrt(aa)*sqrt(bb)5@(djK8\,xY6K
$c6 if (ab.gt.0.8*rab) cc = cc*1.2\0h.LJQ v;qbJx{nC
$c6 if (ab.lt.0.0) cc = cc/1.2
$c6 endif
$c6 if (cc.gt.1.0) cc = 1.0
$c6 err = 0.0en^0sGg#y
$c6 ul = 0.0
%nod
%dof{)S6BAuE9l
   err = err+[ue]**2 WR8i @ M6L(U}-x
   [ue] = [ue]*cc+T9f]w^7v6f
= [u1]+[ue]
ul = ul + **2
%dof&K{)Yr4WC$bD{}!_
%nodWl1lN9Z#@3}`
[du] = 0.0
$c6 write(*,*) 'cc,err,iteration times =',cc,err,itnS2AF(O'a6M s:D@
$c6 write(*,*) 'ab,rab =',ab,rab_@sm![2Ag\
$c6 if (err.lt.1.0e-8 .or. err.lt.1.0e-8*ul .or. itn.gt.500) thenn8lg ^+WmEd
$c6 open(11,file='end',status='new')rS'o8d,W'`?(j
$c6 close(11)
write(s,unod) u,du
$c6 itn=1
$c6 cc=1.0R#u`iZ*B
$c6 else'YA @Q&wj8c
write(s,unod) u,ue
$c6 itn=itn+1
$c6 endif nq1UQ0f5_x1re-|
$c6 rewind(99)
$c6 write(99,*) itn,ccA(DL.i,ub
$c6 close(99)
uoBj'K)X,x(O R a8X-N
endg7Rn"g~
Ml(O/c){H0U;z
[b]cav.gio文件如下：[/b]%}9] jN;yme {
cav[;{4{)Q'A
dnsulxy(fi2Qu^0I

#elemtype q40J SY8~0]q*G
2dxybMtD.I m1G
\?5m{gKO
[b]cav.gcn文件如下：[/b]rw/pcZa
defi
a nellf &
b ell a

startc a
startc b+nJh.T'hVY G
IF EXIST END DEL ENDAE6g-qe/hA%cA
:2
solvc aSR8z+dU2V
IF NOT EXIST END GOTO 2
solvc b

页: [1]

SimWe仿真论坛's Archiver

二维粘性不可压流域中模拟源、汇