dolfin team mailing list archive

[Alexander Jarosch <alexanj@xxxxx>]

Johan Hoffman wrote:

> > On Tue, May 30, 2006 at 10:07:47AM +0000, Alexander Jarosch wrote:
> >    
> >
> > > Hi everybody.
> > >
> > > There was a short discussion about Stokes flow on complex geometries. I
> > > would like to report that a ffc form like the one below works well on
> > > complex geometries, I still only have a factor of 1.2 offset in
> > > magnitude if I compare it with another fem model.  I would like to get
> > > some feedback about this form, whether it is okay or not from the
> > > developers point of view.
> > >
> > > much appreciated,
> > >
> > > Alex
> > >
> > > ffc form -------------------
> > >
> > > scalar = FiniteElement("Lagrange", "triangle", 1)
> > > vector = FiniteElement("Vector Lagrange", "triangle", 1)
> > > system = vector + scalar
> > >
> > > (v, q) = TestFunctions(system)
> > > (u, p) = TrialFunctions(system)
> > >
> > > f = Function(vector)
> > > h = Function(scalar)
> > > nu = Function(scalar)
> > >
> > > beta  = 0.2
> > > delta = beta*h*h
> > >
> > > a = (nu*dot(grad(v), grad(u)) - div(v)*p + nu*q*div(u) +
> > > delta*dot(grad(q), grad(p)))*dx
> > > L = dot(v + mult(delta, grad(q)), f)*dx
> > >      
> > I think it looks a little strange. The difference from what you had
> > before seems to be that you have included the viscosity in the
> > incompressibility term:
> >
> >    nu*q*div(u)
> >
> > Is this what makes the difference? If so, I don't know why this
> > helps. It's certainly ok to put it there since you're solving for
> > div(u) = 0 anyway, so nu*div(u) = 0 is also correct.
> >
> > In some way, it corresponds to a smaller penalty for compressibility
> > when the viscosity is small which might be reasonable.
> >
> > /Anders
> >    
>
> I think this term is somewhat strange, and I haven't seen it before. I do
> not see why the viscosity should multiply the incompressibility.
>
> If the problem is connected to the stability of the incompressibility
> constraint, I propose to instead add to the left hand side of the form a
> stabilization term for the incompressibility, of the form:
>
> delta2*div(v)*div(u)
>
> with delta2 = C2*h*h, and C2 of unit order.
>
> Also, what is the value of your parameter beta in the definition of your
> delta? It might be worthhile to try to experiment with higher values of
> beta.
>
> /Johan
>
>
>
>  
>
> > _______________________________________________
> > DOLFIN-dev mailing list
> > DOLFIN-dev@xxxxxxxxxx
> > http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> >
> >    
>
>
>
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
>
>  
Thanks for all the input. The idea of multiplying the viscosity with the 
incompressibility was also strange to me, but it seems to make the 
difference.
I tried a form proposed by Johan like this (is that what you had in mind?):

beta  = 0.4
beta2 = 0.4
delta = beta*h*h
delta2 = beta2*h*h

a = (nu*dot(grad(v), grad(u)) - div(v)*p + q*div(u) + delta*dot(grad(q), 
grad(p)) + delta2*div(v)*div(u))*dx
L = dot(v + mult(delta, grad(q)), f)*dx

and also played around with beta and beta2 values. It did not fix the 
problem, which is: I work at very high viscosities, 1e13 Pa s and my 
surface is curved. Now if I use the form above or a normal stokes form, 
the pressure on the surface, where it should be 0 is not 0, the pressure 
is 0 on the endpoints, but not along the curved surface. The pressure 
seems not to realize that the surface is curved. Now if I multiply the 
viscosity into the incompressibility term like

nu*q*div(u)

suddenly the pressure is correct on the surface boundary, and the 
velocities as well, as mentioned.

Any idea why this could be?

Thanks in advance

/Alex