testcase to implement surface tension in biFluids

testcases/liqDrop/liqDrop.py

+# MigFlow - Copyright (C) <2010-2018>
+# <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
+from migflow import fluid
+from migflow import scontact
+
+import numpy as np
+import os
+import time
+import shutil
+import random
+
+outputdir = "output"
+if not os.path.isdir(outputdir) :
+    os.makedirs(outputdir)
+
+t = 0
+ii = 0
+
+
+#physical parameters
+g =  -9.81                
+rhof = 1030                                      # fluid density
+rhop = 2450                                     # grains density
+r=154e-6                                        # grains radii
+R = 3.3e-3                                      # drop radius
+compacity = .20                                 # solid volume fraction in the drop
+phi = 1 - compacity
+nuf = 1.17/rhof                                   # kinematic viscosity
+muf = nuf*rhof                                     # dynamic viscosity
+
+##mixture properties
+rhom = (1-phi)*rhop+phi*rhof                   #mixture density
+num = nuf*(1+5/2*phi)                          #Einstein viscosity
+
+print("rhom = %g, num = %g" % (rhom,num))
+
+
+#numerical parameters
+tEnd = 200                                     # final time
+dt = 1e-2                                       # time step
+outf = 10                                       # number of iterations between output files
+
+def outerBndV(x) :
+    print(.4*max(np.sin(t*np.pi*2./1),0))
+    return 0.4*max(np.sin(t*np.pi*2./1),0)
+
+fluid = fluid.FluidProblem(2,g,[num*rhom,nuf*rhof],[rhom,rhof])
+fluid.load_msh("mesh.msh")
+fluid.set_strong_boundary("Top",2,0)
+fluid.set_strong_boundary("Top",1,0)
+fluid.set_strong_boundary("Bottom",1,0)
+fluid.set_strong_boundary("Lateral",0,0)
+fluid.set_weak_boundary("Bottom","Wall")
+fluid.set_weak_boundary("Lateral","Wall")
+fluid.set_weak_boundary("Top","Wall")
+ii = 0
+t = 0
+
+s = fluid.solution()
+x = fluid.coordinates()
+
+for i in range(len(x[:,0])):
+    if (x[i,0])**2+(x[i,1]-.52)**2<R**2:
+        s[i,3] = 1
+    else:
+        s[i,3] = 0
+
+
+#set initial_condition
+
+fluid.export_vtk(outputdir,0,0)
+
+ii = 0
+tic = time.time()
+while t < tEnd : 
+    #Fluid solver
+    fluid.implicit_euler(dt)
+    if (ii%11==0 and ii != 0):
+       fluid.adapt_mesh(1e-3,1e-4,10000)
+    t += dt
+    #Output files writting
+    if ii %outf == 0 :
+        ioutput = int(ii/outf) + 1
+        fluid.export_vtk(outputdir, t, ioutput)
+    ii += 1
+    print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))

testcases/liqDrop/mesh.geo

+L = 0.05;
+H = 0.6;
+y = 0;
+lc = 0.001;
+
+Point(1) = {-L, H, 0};
+Point(2) = {-L, -H, 0};
+Point(3) = {L, -H, 0};
+Point(4) = {L, H, 0};
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Point(5) = {0,0.52,0};
+
+Line Loop(1) = {1:4};
+
+Plane Surface(1) = {1};
+Physical Line("Bottom") = {2};
+Physical Line("Lateral") = {1,3};
+Physical Line("Top") = {4};
+Physical Surface("Domain") = {1};
+Physical Point("PtFix") = {1};
+
+Field[1] = Attractor;
+Field[1].NodesList = {5};
+
+Field[2] = Threshold;
+Field[2].DistMax = 0.08;
+Field[2].DistMin = 0.02;
+Field[2].LcMax = 0.008;
+Field[2].LcMin = 0.001;
+Field[2].IField = 1;
+
+Background Field = 2;
+
+Mesh.CharacteristicLengthFromPoints = 0;
+Mesh.CharacteristicLengthExtendFromBoundary = 0;