michael 2007-9-4 11:11
Moving boundaries
[size=2]I want to model, in 2-d, a rigid moving circle through a tube of 2-d Poiseuille flow. I have attached the code that I have written to try and accomplish this task. FlexPDE5 complains that a number I use inside of the boundaries section is not a constant. I realize FlexPDE5 may not be written to support nonconstant functions, but ideally I'd like FlexPDE5 to evaluate x0 at whatever time step it is on (it is only a function of time and as such always remains constant during any given time step) and use that constant to construct the boundary of my circle at that time step. How can I do this without calling FlexPDE5 from some other program? I want to use FlexPDE5's time integrator. In addition, I would love it if there was some way to make a condition which, when met, the time integration would stop (rather than the time integration having to be defined beforehand between two constants). Hopefully either FlexPDE5 can do this or can be made to do this. If FlexPDE5 can't do this yet, I imagine it would not be hard for them to be incorporated. [/size]
5e)J r.w0e
E mZ.MSs
[code] { Fill in the following sections, or delete those that are unused.}
a$W9H7B9fv6xivj-B m
TITLE 'Traveling Cylinder '
:M
[#] A+U V
A)gK`Rd%M9k
SELECTr'L?|$sO eBQ
ngrid = 10
*`(_h!fQv V9]0k
VARIABLES
;a"|T ? ENW `9lF
u(0.1) {standard stokes flow variables}
[ SG F4R~
v(0.01)
EHw0[ jzV
p(1)X#G)B
PQ]wj
@ O
*{ A)?m!@
SCALAR VARIABLES
M g"Cr2S
x0 O&]_ lS\y{
Y0v
u0 {horizontal velocity of cylinder which is immersed in the tube}9gLS+K?6y
v0 {vertical velocity of cylinder which is immersed in the tube}V x6i x P.^ ^Q1pd
w0 {angular velocity of the cylinder}ou8b4K#l EMN ~+?
;h8v j ] f OmG!m(g
DEFINITIONS
1Sb;R;DT
mu = 1 {viscosity}"c!T%fNv
n
L = 20 {Length of tube}
|
O;@)B X/Q*w
r = 4 {"radius" of tube, across tube = 2r}7Fp1h2C9dZ@7V.Q
Q = 8 {total "volume" flux in tube (8 square length units per time unit) }
s a3~6a9N
n@
rc = 2 {radius of cylinder}#H)|r-k.]
\
p0 = 0 {A reference pressure for one end of the tube}4a6g {g
ba"S c!S
K = 100 {Penalty constant for "penalty method" for enforcing incompressibility}XJ:dx db|5V x
{ x0 = 5 {current x position of cell}};[%Nv_%z-TY~P
y0 = 0 {current y position of cell}7S;`7j(|:K bMsa
gradp = -3*mu*Q/(2*r^3) {gradient of pressure for background pois flow}
3s5f7_ Az)^.| [U
z X c5z |)p%T
sigma11 = 2*mu*dx(u)-p {stress vector}/GY
t A(~YL%Q3_
h
sigma21 = mu*(dy(u)+dx(v))X-{#OQ~E
sigma12 = sigma21
6T? C
f,?+j0Gm }
sigma22 = 2*mu*dy(v)-p Wb3S"E:z
.nE
{W8t5lx+vB
INITIAL VALUES&b+w }]\
u = 3*Q*(r^2-y^2)/(4*r^3) {unidirectional Poiseuille flow right to left}
g;F0y7`d~7` J n
v = 0 ,U&d5_+y(Wx
p = p0-x*gradp+UZq H9s,@
O
x0 = 5
!MV+X\y1jG'S
u0 = 3*Q*(r^2-y0^2)/(4*r^3) {unidirectional Poiseuille flow right to left}~5fq4k:WT:` nkx
v0 = 0BNf x A4U#a6yIE
w0 = 0
6Y1t+Hh [t
jw]W
o&H
EQUATIONS
)c_4a/C,T~O
u: dx(sigma11)+dy(sigma12) = 0 {Stokes flow}"n:Fm\n(\ ZP*R
P~
v: dx(sigma21)+dy(sigma22) = 0
^ ?%oO |+SO ?+s?,O
p: div(grad(p)) = K*(dx(u)+dy(v)) {Penalty method}
`n4^$|j
Mn/M
kzG:M
x0: dt(x0) = u0
`Wyt5Q
{ u0: u0 =0}3c2r,lha
u0: line_integral(normal(sigma11,sigma12),'cyl_1') = 0 {Force on cylinder must be zero in x and y direction}i$Jc T:G-M
{ v0: v0 = 0}"j yE.G Tvb0I
v0: line_integral(normal(sigma21,sigma22),'cyl_1') = 0 {Torque (next line) also zero}
1]7Z6| j&F7f$B
{ w0: w0 = 0}
{ avoGQ I%o/e
w0: line_integral(cross(vector(normal(sigma11,sigma12),normal(sigma21,sigma22)),vector(x-x0,y-y0)),'cyl_1') = 0
[6z"VAp)y O
X*m.s?irN_
BOUNDARIES2t2F6d3QxA*h
REGION 1
|t X'f-@
START(0,-r) { Bottom, no slip and assumed no gradient of p normal to wall }n"tRd%F2v
value(u) = 0 value(v) = 0 load(p) = 07['T1X@*Cf/b%m1c d
line to (L,-r )
)Fw V^#l t
eRR
value(u) = 3*(r^2-y^2)*Q/(4*r^3) value(v) = 0 load(p) = -3*mu*Q/(2*r^3)
yJ~.R3tcj,z N1s
line to (L,r)
]A/uNP/G2F2R0z{ C!v
value(u) = 0 value(v) = 0 load(p) =0
l$d b["X/h(mA
line to (0,r)
W.T m.GT
load(u) = p0 value(v) = 0 value(p) = p0
U"p+FLKd&S6q
line to close!a a
aKBz(B
5yWP h O7Y ?
EXCLUDE*bw/J!a(`y8r.M1v2OH
start 'cyl_1' (x0+rc,y0)
5t(zp\P"_*?
value(u) = u0-(y-y0)*w0(?Ewn7dV
value(v) = v0-(x-x0)*w0
t I |
@ ~2X/F$\
load(p) = 02M+O+K!S7Gm?8u1Nu&M
arc (center = x0,y0) angle = 360p!^!Y8Sjw
sN:K|D l
MONITORS
l'h"jR4nD g$ur1X1n
contour(u) as 'u'8g)^` ?]:Qz5_
contour(v) as 'v'
r}Y,DO7mW
contour(p) as 'p'c/h0Y1j Y*??"H2~.h3Q
contour(dx(u)+dy(v)) as 'div(velocity)'2u B_-M0d
t r
report u0%K P%H9vD_q])x
report v0
#@RA v.O*x.T(`2r
report w0.R-]EVp~ E f`(PJ
n"M
{({*^;iH!J!D-n(z
PLOTS
q4g2z3V9J~y)~R
@ ~
contour(u) as 'u'
N&Q,\D&@ P
o h
contour(v) as 'v'
0B/tY"N VF*[2u
contour(p) as 'p' f uu p @H;}f
C1T
contour(dx(u)+dy(v)) as 'div(velocity)'
a]N(_
c wx6o|9s1I
report u0
3j1Z.z%z$utudM:a|xx7C
report v02mj,Xo[9C
report w0 g2_3c8onhg9l
(D1A2E-U4nM6]
END|)Xe]%},n
1{X OS M@ ^
j wzDO4q6]#[B
[/code]
michael 2007-9-4 11:11
The BOUNDARIES layout section describes what the INITIAL positions of the boundaries are, and these have to be reducible to concrete numbers at layout time. E7^Kt#B['f@)W^
You can't say that the position depends on what the position is.