square.py 3.32 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# MigFlow - Copyright (C) <2010-2020>
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
# 	List of the contributors to the development of MigFlow: see AUTHORS file.
# Description and complete License: see LICENSE file.
# 	
# This program (MigFlow) is free software: 
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
# Public License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
# 
# You should have received a copy of the GNU Lesser General Public License
# along with this program (see COPYING and COPYING.LESSER files).  If not, 
# see <http://www.gnu.org/licenses/>.

#!/usr/bin/env python
Matthieu Constant's avatar
Matthieu Constant committed
22
23
24
25

# TESTCASE DESCRIPTION
# Transformation of a square of ethanol into a circle due to the tension surface force exerted by the air surrounding the square

26
27
from migflow import fluid
from migflow import scontact
Matthieu Constant's avatar
Matthieu Constant committed
28
from migflow import time_integration
29
30
31
32
33
34
35
36
37
38
39
40

import numpy as np
import os
import time
import shutil
import random

outputdir = "output"
if not os.path.isdir(outputdir) :
    os.makedirs(outputdir)

#physical parameters
Matthieu Constant's avatar
Matthieu Constant committed
41
g = np.array([0,0])               
Matthieu Constant's avatar
Matthieu Constant committed
42
43
44
45
46
47
rho0 = 1.117                                            # air density
rho1 = 785.92                                           # ethanol density
nu0 = 1.57e-5                                           # air kinematic viscosity
nu1 = 1.2e-3/rho1                                       # ethanol kinematic viscosity
R = 0.01*np.pi**0.5/2*10                                # square lenght
sigma = 22.39e-3                                        # surface tension
48
49

#numerical parameters
Matthieu Constant's avatar
Matthieu Constant committed
50
51
52
53
54
55
56
57
58
59
tEnd = 200                                              # final time
dt = 1e-3                                               # time step
outf = 10                                               # number of iterations between output files
ii = 0                                                  # initial iteration
t = 0                                                   # initial time

#
# FLUID PROBLEM
#
fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1],sigma=sigma)
60
61
62
63
fluid.load_msh("mesh.msh")
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
fluid.set_wall_boundary("Top",pressure=0)
Matthieu Constant's avatar
Matthieu Constant committed
64
65
66
67
68
69
70
71
fluid.set_strong_boundary("Top",0,0)
fluid.set_strong_boundary("Top",1,0)
fluid.set_strong_boundary("Bottom",0,0)
fluid.set_strong_boundary("Bottom",1,0)
fluid.set_strong_boundary("Lateral",0,0)
fluid.set_strong_boundary("Lateral",1,0)

# Solution fields and coordinates
72
73
74
75
s = fluid.solution()
c = np.ndarray((fluid.n_nodes()))
x = fluid.coordinates()

Matthieu Constant's avatar
Matthieu Constant committed
76
77
78
# Initial concentration field
c[:] = 0
c[np.logical_and(np.abs(x[:,0])<R,np.abs(x[:,1])<R)] = 1
79
fluid.set_concentration_cg(c)
Nathan Coppin's avatar
Nathan Coppin committed
80
fluid.write_vtk(outputdir,0,0)
81

Matthieu Constant's avatar
Matthieu Constant committed
82

83
tic = time.time()
Matthieu Constant's avatar
Matthieu Constant committed
84
85
86
#
# COMPUTATION LOOP
#
87
while t < tEnd : 
Matthieu Constant's avatar
Matthieu Constant committed
88
    time_integration.iterate(fluid,None,dt)
89
90
91
92
93
    t += dt
    #Output files writting
    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
    if ii %outf == 0 :
        ioutput = int(ii/outf) + 1
Nathan Coppin's avatar
Nathan Coppin committed
94
        fluid.write_vtk(outputdir, ioutput, t)
95
    ii += 1