16. Fluid-Structure-interaction#
Michael Neunteufel + JS ‘21
coupling of non-linear elasticity and Navier-Stokes on moving domain
strong coupling of velocities via Lagrange parameter
mesh deformation keeps div-free functions div-free
Piola transform causes another term like mesh-velocity
from netgen.geom2d import CSG2d, Circle, Rectangle
from ngsolve import *
from ngsolve.webgui import Draw
tau = 0.005
tend = 5.5
order = 3
alpha = Parameter(5)
inverse = "umfpack"
truecompile = True
(rhos, nus, mus, rhof, nuf, U) = (1e4, 0.4, 0.5*1e6, 1e3, 1e-3, 1)
ls = 2*mus*nus/(1-2*nus)
uinflow = CF( (4*U*1.5*y*(0.41-y)/(0.41*0.41), 0) )
def GenerateMesh(order, maxh=0.12, L=2.5, H=0.41):
r = 0.05
D = 0.1
geo = CSG2d()
circle = Circle( center=(0.2,0.2), radius=0.05, bc="circ" )
fluid = Rectangle( pmin=(0,0), pmax=(L,H), mat="fluid", left="inlet", top="wall", bottom="wall", right="outlet" )
solid = Rectangle( pmin=(0.248989794855664,0.19), pmax=(0.6,0.21), mat="solid", left="circ_inner", top="interface", bottom="interface", right="interface" )
domain_fluid = (fluid - circle) - solid
domain_solid = (solid - circle)
geo.Add(domain_fluid)
geo.Add(domain_solid)
# generate mesh
mesh = Mesh(geo.GenerateMesh(maxh=maxh))
mesh.Curve(order)
return mesh
mesh = GenerateMesh(order=order)
Draw(mesh)
V1 = VectorH1( mesh, order=order, dirichlet = "wall|outlet|inlet|circ|circ_inner" )
V2 = HDiv( mesh, order=order, dirichlet = "wall|inlet|circ|circ_inner", definedon = "fluid", hodivfree=True )
V3 = TangentialFacetFESpace(mesh, order=order, dirichlet="wall|inlet|circ|circ_inner|outlet", definedon = "fluid" )
V4 = VectorH1( mesh, order=order, dirichlet="circ|circ_inner", definedon = "solid" )
Q = L2( mesh, order=0, definedon = "fluid", lowest_order_wb=True )
V5 = SurfaceL2( mesh, order=order, dirichlet="wall|outlet|inlet|circ|circ_inner" )
V = V2*V3*Q*V1*V4*V5*V5
Y = V2*V3*Q
H = VectorH1( mesh, order=order, dirichlet = "wall|inlet|outlet|circ|circ_inner" )
u, uhat, p, d, v, fx, fy = V.TrialFunction()
ut,uhatt,pt,dt,vt,fxt,fyt = V.TestFunction()
f, ft = CF( (fx,fy) ), CF( (fxt,fyt) )
gfsol = GridFunction(V) #solution
gfsolold = GridFunction(V) #oldsolution
gfset = GridFunction(H) #for SetDeformation
gfu, gfvelhat, pressure, gfd, gfv, gffx, gffy = gfsol.components
uold, uhatold, pold, dold, vold, fxold, fyold = gfsolold.components
I = Id(mesh.dim)
F = grad(d) + I
C = F.trans*F
Fold = grad(dold) + I
Cold = Fold.trans*Fold
n,h = specialcf.normal(2),specialcf.mesh_size
def tang(vec):
return vec - (vec*n)*n
def norma(vec):
return (vec*n)*n
def Stress(mat):
return mus*mat+ls/2*Trace(mat)*I
#For Stokes problem
stokes = BilinearForm( Y, symmetric = True, check_unused=False )
stokes += (nuf*rhof*InnerProduct( 2*Sym(grad(u)), Sym(grad(ut))) \
- div(u)*pt - div(ut)*p - 1e-8*p*pt )*dx("fluid")
stokes += (nuf*rhof*(InnerProduct( 2*Sym(grad(u))*n, tang(uhatt-ut) ) \
+InnerProduct( 2*Sym(grad(ut))*n, tang(uhat-u) ) \
+alpha*order*order/h * InnerProduct( tang(uhatt-ut), tang(uhat-u) ))) \
*dx("fluid", element_boundary=True)
stokes.Assemble()
######################################### SOLID #############################################
a = BilinearForm(V, symmetric=False, check_unused=False, condense=True)
a += ( (1/tau*(d-dold)-0.5*(v+vold))*dt +rhos/tau*(v-vold)*vt +2*InnerProduct(F*Stress(0.5*(0.5*(C+Cold)-I)), grad(vt)) ).Compile(truecompile,True)*dx("solid")
factor = 1e-22*mus
gfdist = GridFunction(H1(mesh, order=2))
def minCF(a,b) : return IfPos(a-b, b, a)
gfdist.Set( minCF( (x-0.6)**2+(y-0.19)**2, (x-0.6)**2+(y-0.21)**2 ) )
def NeoHookExt (C, mu=1,lam=1):
return 0.5 * mu * (Trace(C-I) + 2*mu/lam * Det(C)**(-lam/2/mu) - 1)
a += Variation( factor/sqrt(gfdist*gfdist+1e-12)*NeoHookExt(C)*dx("fluid")).Compile(truecompile,True)
################################# FLUID ##############################################
graddp = 1/(tau)*(grad(d)-grad(dold))
relv = u-1/tau*(d-dold)
relvn = relv*n
upwind = (u*n)*n + (IfPos(relvn, tang(u), tang(uhat)))
cfinner = nuf*rhof*InnerProduct(Sym(grad(u))+Sym(grad(uold)),Sym(grad(ut))) - div(u)*pt-div(ut)*p -rhof*InnerProduct(0.5*(u+uold), grad(ut).trans*relv) + rhof*InnerProduct(0.5*graddp*(u+uold),ut) + rhof/tau*(u-uold)*ut
cfouter = nuf*rhof*(InnerProduct(Sym(grad(u)+grad(uold))*n,tang(uhatt-ut)) + InnerProduct(2*Sym(grad(ut))*n,0.5*tang(uhat-u + uhatold-uold)) + alpha*order*order/h*InnerProduct(tang(uhatt-ut),0.5*tang(uhat-u+uhatold-uold)))+rhof*(relvn*upwind*ut+IfPos(relvn,1,0)*0.5*relvn*tang(uhat-u+uhatold-uold)*uhatt)
a += cfinner.Compile(truecompile,True)*dx("fluid",deformation=gfset)
a += cfouter.Compile(truecompile,True)*dx("fluid",element_boundary=True, deformation=gfset)
#interface
a += ((norma(v-u.Trace())+tang(v-uhat.Trace()))*ft + (norma(vt-ut.Trace())+tang(vt-uhatt.Trace()))*f).Compile(truecompile,True)*ds("interface", deformation=gfset)
bts = Y.FreeDofs() & ~Y.GetDofs(mesh.Materials("solid"))
bts &= ~Y.GetDofs(mesh.Boundaries("wall|inlet|circ|interface|circ_inner"))
invstoke = stokes.mat.Inverse(bts, inverse = "umfpack")
tmp = GridFunction(Y)
rstokes = tmp.vec.CreateVector()
tmp.components[0].Set( uinflow, definedon=mesh.Boundaries("inlet") )
rstokes.data = stokes.mat*tmp.vec
tmp.vec.data -= invstoke*rstokes
gfsol.components[0].vec.data = tmp.components[0].vec
gfsol.components[1].vec.data = tmp.components[1].vec
gfsol.components[2].vec.data = tmp.components[2].vec
w = gfsol.vec.CreateVector()
res = gfsol.vec.CreateVector()
gfsol.Load("fsisol4.sol")
scene = Draw(mesh.MaterialCF( {"fluid" :gfu}, default=None), mesh.Materials("fluid"), \
deformation=gfd, order=3)
t = 3
tend=4
SetNumThreads(4)
with TaskManager():
while t < tend-tau/2.0:
t += tau
print("\rt = ", t, end="")
gfsolold.vec.data = gfsol.vec
for step in range(3):
gfset.Set(gfsol.components[3])
a.AssembleLinearization(gfsol.vec)
inv = a.mat.Inverse(a.space.FreeDofs(a.condense), inverse=inverse)
a.Apply (gfsol.vec, res)
if a.condense:
res.data += a.harmonic_extension_trans * res
w.data = inv*res
if a.condense:
w.data += a.harmonic_extension * w
w.data += a.inner_solve * res
err = InnerProduct(w, res)
# print ("err^2 = ", err)
gfsol.vec.data -= w
if abs(err) < 1e-20:
break;
scene.Redraw()
# gfsol.Save("fsisol4.sol")
t = 3.005
t = 3.01
t = 3.0149999999999997
t = 3.0199999999999996
t = 3.0249999999999995
t = 3.0299999999999994
t = 3.0349999999999993
t = 3.039999999999999
t = 3.044999999999999
t = 3.049999999999999
t = 3.054999999999999
t = 3.0599999999999987
t = 3.0649999999999986
t = 3.0699999999999985
t = 3.0749999999999984
t = 3.0799999999999983
t = 3.084999999999998
t = 3.089999999999998
t = 3.094999999999998
t = 3.099999999999998
t = 3.1049999999999978
t = 3.1099999999999977
t = 3.1149999999999975
t = 3.1199999999999974
t = 3.1249999999999973
t = 3.1299999999999972
t = 3.134999999999997
t = 3.139999999999997
t = 3.144999999999997
t = 3.149999999999997
t = 3.1549999999999967
t = 3.1599999999999966
t = 3.1649999999999965
t = 3.1699999999999964
t = 3.1749999999999963
t = 3.179999999999996
t = 3.184999999999996
t = 3.189999999999996
t = 3.194999999999996
t = 3.1999999999999957
t = 3.2049999999999956
t = 3.2099999999999955
t = 3.2149999999999954
t = 3.2199999999999953
t = 3.224999999999995
t = 3.229999999999995
t = 3.234999999999995
t = 3.239999999999995
t = 3.2449999999999948
t = 3.2499999999999947
t = 3.2549999999999946
t = 3.2599999999999945
t = 3.2649999999999944
t = 3.2699999999999942
t = 3.274999999999994
t = 3.279999999999994
t = 3.284999999999994
t = 3.289999999999994
t = 3.2949999999999937
t = 3.2999999999999936
t = 3.3049999999999935
t = 3.3099999999999934
t = 3.3149999999999933
t = 3.319999999999993
t = 3.324999999999993
t = 3.329999999999993
t = 3.334999999999993
t = 3.3399999999999928
t = 3.3449999999999926
t = 3.3499999999999925
t = 3.3549999999999924
t = 3.3599999999999923
t = 3.364999999999992
t = 3.369999999999992
t = 3.374999999999992
t = 3.379999999999992
t = 3.384999999999992
t = 3.3899999999999917
t = 3.3949999999999916
t = 3.3999999999999915
t = 3.4049999999999914
t = 3.4099999999999913
t = 3.414999999999991
t = 3.419999999999991
t = 3.424999999999991
t = 3.429999999999991
t = 3.4349999999999907
t = 3.4399999999999906
t = 3.4449999999999905
t = 3.4499999999999904
t = 3.4549999999999903
t = 3.45999999999999
t = 3.46499999999999
t = 3.46999999999999
t = 3.47499999999999
t = 3.4799999999999898
t = 3.4849999999999897
t = 3.4899999999999896
t = 3.4949999999999894
t = 3.4999999999999893
t = 3.5049999999999892
t = 3.509999999999989
t = 3.514999999999989
t = 3.519999999999989
t = 3.524999999999989
t = 3.5299999999999887
t = 3.5349999999999886
t = 3.5399999999999885
t = 3.5449999999999884
t = 3.5499999999999883
t = 3.554999999999988
t = 3.559999999999988
t = 3.564999999999988
t = 3.569999999999988
t = 3.5749999999999877
t = 3.5799999999999876
t = 3.5849999999999875
t = 3.5899999999999874
t = 3.5949999999999873
t = 3.599999999999987
t = 3.604999999999987
t = 3.609999999999987
t = 3.614999999999987
t = 3.619999999999987
t = 3.6249999999999867
t = 3.6299999999999866
t = 3.6349999999999865
t = 3.6399999999999864
t = 3.6449999999999863
t = 3.649999999999986
t = 3.654999999999986
t = 3.659999999999986
t = 3.664999999999986
t = 3.6699999999999857
t = 3.6749999999999856
t = 3.6799999999999855
t = 3.6849999999999854
t = 3.6899999999999853
t = 3.694999999999985
t = 3.699999999999985
t = 3.704999999999985
t = 3.709999999999985
t = 3.7149999999999848
t = 3.7199999999999847
t = 3.7249999999999845
t = 3.7299999999999844
t = 3.7349999999999843
t = 3.7399999999999842
t = 3.744999999999984
t = 3.749999999999984
t = 3.754999999999984
t = 3.759999999999984
t = 3.7649999999999837
t = 3.7699999999999836
t = 3.7749999999999835
t = 3.7799999999999834
t = 3.7849999999999833
t = 3.789999999999983
t = 3.794999999999983
t = 3.799999999999983
t = 3.804999999999983
t = 3.8099999999999827
t = 3.8149999999999826
t = 3.8199999999999825
t = 3.8249999999999824
t = 3.8299999999999823
t = 3.834999999999982
t = 3.839999999999982
t = 3.844999999999982
t = 3.849999999999982
t = 3.8549999999999818
t = 3.8599999999999817
t = 3.8649999999999816
t = 3.8699999999999815
t = 3.8749999999999813
t = 3.8799999999999812
t = 3.884999999999981
t = 3.889999999999981
/Applications/Netgen.app/Contents/Resources/lib/python3.13/site-packages/netgen/webgui.py:20: RuntimeWarning: overflow encountered in cast
values = np.array(data.flatten(), dtype=dtype)
t = 3.894999999999981
t = 3.899999999999981
---------------------------------------------------------------------------
NgException Traceback (most recent call last)
Cell In[3], line 18
16 for step in range(3):
17 gfset.Set(gfsol.components[3])
---> 18 a.AssembleLinearization(gfsol.vec)
19 inv = a.mat.Inverse(a.space.FreeDofs(a.condense), inverse=inverse)
20 a.Apply (gfsol.vec, res)
NgException: Inverse matrix: Matrix singularin AssembleLinearization