title '3D flow  in  muzzle'
                                 {uint g-mm-ms}
coordinates
cartesian3
 
variables
vx(1e-8)
vy(1e-8)
vz(1e-8)
p
 
select
ngrid=20
stages=2
regrid=off
 
 
definitions    {define the solution domain and field}
 
r1=5    {define the radius of the arc}
r2=217.7
r3=180
r4=186.78
r5=5
r6=5
r7=15.5
 
ox=0            oy=0
o1x=0          o1y=108
 
ax=-r7*cos(16*PI/180)     ay=o1y+r7*sin(16*PI/180)
                                       {precisely define the point on the circle perimeter}
bxx=-17.2     byy=117.9
cx=0             cy=135.2
d0x=98.7      d0y=150.5
ex=134.8      ey=119.3
fx=20.8         fy=119.3
gx=17.2        gy=117.9
 
 
 hx=r7*cos(16*PI/180)         hy=o1y+r7*sin(16*PI/180)
            {precisely define the point on the circle perimeter}
 
h1=20   {define the extrusion length of 3 parts}
h2=90
h3=20
 
 
dens=staged(0, 150e-6 )     { fluid density }
 
visc= 1.8e-8                       { fluid viscosity}
 
v=vector(vx,vy,vz)
 
vm=magnitude(v)
 
div_v = dx(vx) + dy(vy) + dz(vz)
 
PENALTY = 150e-6*1000
 
Pin = 40                                  {input pressure 40mpa}
Pout = 0.1                               {input pressure 0.1mpa}
 
INITIAL VALUES
vx=0
vy=0
vz=0
p=Pin+(Pout-Pin)*z/(h1+h2+h3)
                                                  { start with a rough approximation to the pressure }
                                                  {give the initial value for pressure}
 
EQUATIONS
vx:   dens*(vx*dx(vx) + vy*dy(vx) + vz*dz(vx)) + dx(p) -visc*div(grad(vx)) = 0
vy:   dens*(vx*dx(vy) + vy*dy(vy) + vz*dz(vy)) + dy(p) -visc*div(grad(vy)) = 0
vz:   dens*(vx*dx(vz) + vy*dy(vz) + vz*dz(vz)) + dz(p) -visc*div(grad(vz)) = 0
p:     div(grad(p)) = PENALTY*div_v
 
extrusion         {define the extrusion of the geometry of the body}
     !  surface 1   z=0
     !   layer    1
     ! surface  2  z=h1
     !  layer   2
     ! surface  3  z=(h1+h2)
     ! layer   3
     !surface  4   z=(h1+h2+h3)
      z =0, h1, h1+h2, h1+h2+h3
 
boundaries
Region 1     { whole cylinder }
 
     surface 2 value(vx)=0 value(vy)=0 value(vz)=0    natural(p)=0
 
     surface 3 value(vx)=0 value(vy)=0 value(vz)=0    natural(p)=0
     layer  1  void
     layer  3  void
 
   start(hx,hy)
 
                 value(vx)=0 value(vy)=0 value(vz)=0 natural(p)=0
 
   arc(center=o1x,o1y)  angle=-212  {from H to the point of A}
 
                 value(vx)=0 value(vy)=0 value(vz)=0  natural(p)=0
 
   arc(radius=r1)  to(bxx,byy)   {counter-clockwise  description the arc}
 
                 value(vx)=0 value(vy)=0 value(vz)=0   natural(p)=0
 
   line to(cx,cy)
 
                value(vx)=0 value(vy)=0 value(vz)=0   natural(p)=0
 
   arc(radius=r2)  to(d0x,d0y)
 
              layer 2  natural(vx)=0 natural(vy)=0 natural(vz)=0   value(p)=Pout
               { duct sidewall drag }
   arc(radius=r3)  to(ex,ey)
 
              value(vx)=0 value(vy)=0 value(vz)=0   natural(p)=0
 
   arc(radius=r4)  to(fx,fy)
 
             value(vx)=0 value(vy)=0 value(vz)=0   natural(p)=0
 
   arc(radius=r5)  to(gx,gy)
 
              value(vx)=0 value(vy)=0 value(vz)=0   natural(p)=0
 
   arc(radius=r6)  to(hx,hy)
 
   to close
 
 
limited Region 2               { main hole }
  layer  1
  surface 1
   !value(vz)=800
  natural(vx)=0  natural(vy)=0  natural(vz)=0    value(p)=Pin       { input duct opening }
  start(r7,o1y)
   layer 1  value(vx)=0 value(vy)=0 value(vz)=0    natural(p)=0         { duct sidewall drag }
   arc(center=o1x,o1y)  angle=360
   to close
 
 limited Region 3     {sidewall  nozzle}
 layer 3
   !surface 4 natural(vx)=0 natural(vy)=0 natural(vz)=0 value(p)=Pout { output duct opening }
 surface 4    natural(vx)=0 natural(vy)=0 natural(vz)=0  value(p)=Pout    { input duct opening }
 start(r7,o1y)
   layer 3    value(vx)=0 value(vy)=0 value(vz)=0    natural(p)=0             { duct sidewall drag }
   arc(center=o1x,o1y)  angle=360
 to close
 
 
monitors
 
contour(vz)    on z=25   report dens report pin
contour(vz)    on z=65   report dens report pin
contour(vx)    on z=25   report dens report pin
contour(vx)    on z=65   report dens report pin
 
contour(vx)    on z=130  zoom(-20,90,40,40)  report dens report pin
contour(vz)    on z=130  zoom(-20,90,40,40)  report dens report pin
contour(vz)    on x= 0       report dens report pin
 
 
vector(vx,vz) on y=o1y     report dens report pin
contour(vz)   on y=o1y    report dens report pin
contour(vy)   on y=o1y     report dens report pin
contour(p)     on z=65
 
 
plots
 
contour(vz)    on z=25   report dens report pin
contour(vz)    on z=65   report dens report pin
 
contour(vx)    on z=25   report dens report pin
contour(vx)    on z=65   report dens report pin
 
contour(vx)    on z=130  zoom(-20,90,40,40)  report dens report pin
contour(vz)    on z=130 zoom(-20,90,40,40)  report dens report pin
contour(vz)    on x= 0     report dens report pin
 
vector(vx,vz) on  y=o1y     report dens report pin
contour(vz)    on  y=o1y     report dens report pin
contour(vy)    on  y=o1y     report dens report pin
contour(p)      on  z=65
 
end