| 1 |
from esys.escript import * |
| 2 |
import esys.finley |
| 3 |
from esys.escript.models import StokesProblemCartesian |
| 4 |
from esys.finley import finley |
| 5 |
from esys.finley import Rectangle |
| 6 |
from LevelSet import * |
| 7 |
|
| 8 |
#physical properties |
| 9 |
rho1 = 1000 #fluid density on bottom |
| 10 |
rho2 = 1010 #fluid density on top |
| 11 |
eta1 = 100.0 #fluid viscosity on bottom |
| 12 |
eta2 = 100.0 #fluid viscosity on top |
| 13 |
g=10.0 |
| 14 |
|
| 15 |
#solver settings |
| 16 |
dt = 0.001 |
| 17 |
t_step = 0 |
| 18 |
t_step_end = 2000 |
| 19 |
TOL = 1.0e-5 |
| 20 |
max_iter=400 |
| 21 |
verbose=True |
| 22 |
useUzawa=True |
| 23 |
|
| 24 |
#define mesh |
| 25 |
l0=0.9142 |
| 26 |
l1=1.0 |
| 27 |
n0=200 |
| 28 |
n1=200 |
| 29 |
|
| 30 |
mesh=Rectangle(l0=l0, l1=l1, order=2, n0=n0, n1=n1) |
| 31 |
#get mesh dimensions |
| 32 |
numDim = mesh.getDim() |
| 33 |
#get element size |
| 34 |
h = Lsup(mesh.getSize()) |
| 35 |
|
| 36 |
#level set parameters |
| 37 |
tolerance = 1.0e-6 |
| 38 |
reinit_max = 30 |
| 39 |
reinit_each = 3 |
| 40 |
alpha = 1 |
| 41 |
smooth = alpha*h |
| 42 |
|
| 43 |
#boundary conditions |
| 44 |
x = mesh.getX() |
| 45 |
#left + bottom + right + top |
| 46 |
b_c = whereZero(x[0])*[1.0,0.0] + whereZero(x[1])*[1.0,1.0] + whereZero(x[0]-l0)*[1.0,0.0] + whereZero(x[1]-l1)*[1.0,1.0] |
| 47 |
|
| 48 |
velocity = Vector(0.0, ContinuousFunction(mesh)) |
| 49 |
pressure = Scalar(0.0, ContinuousFunction(mesh)) |
| 50 |
Y = Vector(0.0,Function(mesh)) |
| 51 |
|
| 52 |
#define initial interface between fluids |
| 53 |
xx = mesh.getX()[0] |
| 54 |
yy = mesh.getX()[1] |
| 55 |
func = Scalar(0.0, ContinuousFunction(mesh)) |
| 56 |
h_interface = Scalar(0.0, ContinuousFunction(mesh)) |
| 57 |
h_interface = h_interface + (0.02*cos(math.pi*xx/l0) + 0.2) |
| 58 |
func = yy - h_interface |
| 59 |
func_new = func.interpolate(ReducedSolution(mesh)) |
| 60 |
|
| 61 |
#Stokes Cartesian |
| 62 |
solution=StokesProblemCartesian(mesh,debug=True) |
| 63 |
solution.setTolerance(TOL) |
| 64 |
solution.setSubProblemTolerance(TOL**2) |
| 65 |
|
| 66 |
#level set |
| 67 |
levelset = LevelSet(mesh, func_new, reinit_max, reinit_each, tolerance, smooth) |
| 68 |
|
| 69 |
while t_step <= t_step_end: |
| 70 |
#update density and viscosity |
| 71 |
rho = levelset.update_parameter(rho1, rho2) |
| 72 |
eta = levelset.update_parameter(eta1, eta2) |
| 73 |
|
| 74 |
#get velocity and pressure of fluid |
| 75 |
Y[1] = -rho*g |
| 76 |
solution.initialize(fixed_u_mask=b_c,eta=eta,f=Y) |
| 77 |
velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,useUzawa=useUzawa) |
| 78 |
|
| 79 |
#update the interface |
| 80 |
func = levelset.update_phi(velocity, dt, t_step) |
| 81 |
|
| 82 |
print "##########################################################" |
| 83 |
print "time step:", t_step, " completed with dt:", dt |
| 84 |
print "Velocity: min =", inf(velocity), "max =", Lsup(velocity) |
| 85 |
print "##########################################################" |
| 86 |
|
| 87 |
#save interface, velocity and pressure |
| 88 |
saveVTK("phi2D.%2.4i.vtu"%t_step,interface=func,velocity=velocity,pressure=pressure) |
| 89 |
#courant condition |
| 90 |
dt = 0.4*Lsup(mesh.getSize())/Lsup(velocity) |
| 91 |
t_step += 1 |
Parent Directory
| 
Revision Log