| 240 |
\begin{figure} |
\begin{figure} |
| 241 |
\center |
\center |
| 242 |
\scalebox{0.7}{\includegraphics{figures/RT2Dsetup.eps}} |
\scalebox{0.7}{\includegraphics{figures/RT2Dsetup.eps}} |
| 243 |
\caption{Parameters, initial interface and boundary conditions for the Rayleigh-Taylor instability problem. The interface is defined as $\phi=0.02cos(\frac{\pi x}{\lambda}) + 0.2$. The fluids have been assigned different densities and equal viscosity (isoviscous).} |
\caption{Parameters, initial interface and boundary conditions for the Rayleigh-Taylor instability problem. The interface is defined as $\phi=0.02cos(\frac{\pi x}{\lambda}) + 0.2$. The fluids have been assigned different densities and equal viscosity (isoviscous) \cite{BOURGOUIN2006}.} |
| 244 |
\label{RT2DSETUP} |
\label{RT2DSETUP} |
| 245 |
\end{figure} |
\end{figure} |
| 246 |
|
|
| 248 |
|
|
| 249 |
\begin{python} |
\begin{python} |
| 250 |
|
|
| 251 |
|
from esys.escript import * |
| 252 |
|
import esys.finley |
| 253 |
|
from esys.escript.models import StokesProblemCartesian |
| 254 |
|
from esys.finley import finley |
| 255 |
|
from LevelSet import * |
| 256 |
|
|
| 257 |
|
#physical properties |
| 258 |
|
rho1 = 1000 #fluid density on bottom |
| 259 |
|
rho2 = 1010 #fluid density on top |
| 260 |
|
eta1 = 100.0 #fluid viscosity on bottom |
| 261 |
|
eta2 = 100.0 #fluid viscosity on top |
| 262 |
|
penalty = 100.0 |
| 263 |
|
g=10.0 |
| 264 |
|
|
| 265 |
|
#solver settings |
| 266 |
|
dt = 0.001 |
| 267 |
|
t_step = 0 |
| 268 |
|
t_step_end = 2000 |
| 269 |
|
TOL = 1.0e-5 |
| 270 |
|
max_iter=400 |
| 271 |
|
verbose=True |
| 272 |
|
useUzawa=True |
| 273 |
|
|
| 274 |
|
#define mesh |
| 275 |
|
l0=0.9142 |
| 276 |
|
l1=1.0 |
| 277 |
|
n0=100 |
| 278 |
|
n1=100 |
| 279 |
|
|
| 280 |
|
mesh=esys.finley.Rectangle(l0=l0, l1=l1, order=2, n0=n0, n1=n1) |
| 281 |
|
#get mesh dimensions |
| 282 |
|
numDim = mesh.getDim() |
| 283 |
|
#get element size |
| 284 |
|
h = Lsup(mesh.getSize()) |
| 285 |
|
print "element size",h |
| 286 |
|
|
| 287 |
|
#level set parameters |
| 288 |
|
tolerance = 1.0e-6 |
| 289 |
|
reinit_max = 30 |
| 290 |
|
reinit_each = 3 |
| 291 |
|
alpha = 1 |
| 292 |
|
smooth = alpha*h |
| 293 |
|
|
| 294 |
|
#boundary conditions |
| 295 |
|
x = mesh.getX() |
| 296 |
|
#left + bottom + right + top |
| 297 |
|
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] |
| 298 |
|
|
| 299 |
|
velocity = Vector(0.0, ContinuousFunction(mesh)) |
| 300 |
|
pressure = Scalar(0.0, ContinuousFunction(mesh)) |
| 301 |
|
Y = Vector(0.0,Function(mesh)) |
| 302 |
|
|
| 303 |
|
#define initial interface between fluids |
| 304 |
|
xx = mesh.getX()[0] |
| 305 |
|
yy = mesh.getX()[1] |
| 306 |
|
func = Scalar(0.0, ContinuousFunction(mesh)) |
| 307 |
|
h_interface = Scalar(0.0, ContinuousFunction(mesh)) |
| 308 |
|
h_interface = h_interface + (0.02*cos(math.pi*xx/l0) + 0.2) |
| 309 |
|
func = yy - h_interface |
| 310 |
|
func_new = func.interpolate(ReducedSolution(mesh)) |
| 311 |
|
|
| 312 |
|
#Stokes cartesian |
| 313 |
|
solution=StokesProblemCartesian(mesh,debug=True) |
| 314 |
|
solution.setTolerance(TOL) |
| 315 |
|
solution.setSubProblemTolerance(TOL**2) |
| 316 |
|
|
| 317 |
|
#level set |
| 318 |
|
levelset = LevelSet(mesh, func_new, reinit_max, reinit_each, tolerance, smooth) |
| 319 |
|
|
| 320 |
|
while t_step <= t_step_end: |
| 321 |
|
#update density and viscosity |
| 322 |
|
rho = levelset.update_parameter(rho1, rho2) |
| 323 |
|
eta = levelset.update_parameter(eta1, eta2) |
| 324 |
|
|
| 325 |
|
#get velocity and pressue of fluid |
| 326 |
|
Y[1] = -rho*g |
| 327 |
|
solution.initialize(fixed_u_mask=b_c,eta=eta,f=Y) |
| 328 |
|
velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,useUzawa=useUzawa) |
| 329 |
|
|
| 330 |
|
#update the interface |
| 331 |
|
func = levelset.update_phi(velocity, dt, t_step) |
| 332 |
|
|
| 333 |
|
print "##########################################################" |
| 334 |
|
print "time step:", t_step, " completed with dt:", dt |
| 335 |
|
print "Velocity: min =", inf(velocity), "max =", Lsup(velocity) |
| 336 |
|
print "##########################################################" |
| 337 |
|
|
| 338 |
|
#save interface, velocity and pressure |
| 339 |
|
saveVTK("phi2D.%2.4i.vtu"%t_step,interface=func,velocity=velocity,pressure=pressure) |
| 340 |
|
#courant condition |
| 341 |
|
dt = 0.4*Lsup(mesh.getSize())/Lsup(velocity) |
| 342 |
|
t_step += 1 |
| 343 |
|
|
| 344 |
\end{python} |
\end{python} |
Parent Directory
| 
Revision Log
| 
Patch