Incompressible Viscoelastic Flows and the Upper Convected Derivative (d3)

[Mike Norton]

Hi Folks,

I'm modeling an incompressible viscoelastic flow at low Re.  The Left Cauchy strain tensor B=F transpose(F), (where F is the deformation gradient) is evolved by the upper convected derivative, some "diffusion", and strain relaxation (with rate kB):

 problem.add_equation("dt(B) -Diff*lap(B) +kB*B = kB*I -dot(vp,grad(B)) + dot(transpose(grad(vp)),B) + dot(B,grad(vp))")

Ideally, because of incompressibililty, it should be the case that J=sqrt(det(B))=det(F)=1. However, over time (and sometimes rapidly), I find that J > 1 even though div(v)=0 continues to be satisfied at all time steps.  My interpretation is that I'm accumulating errors from div(v) not being exactly zero over time, though there is some evidence to the contrary.**

I'm curious if anyone else has used the upper convicted derivate in their models and found the same.

One solution might be to break the dynamics of B up component wise, and only evolve Bxx and Bxy while using the definitions Byx=Bxy and Byy=(J**2-Bxy*Byx)/Bxx to enforce the constraint. I haven't done this yet because I don't know how to create a substitution rule for a vector/tensor. That is, the following returns an error:

    B = dist.TensorField((coords, coords), name='B', bases=(xbasis,ybasis))
    B_xx = dist.Field(name='B_xx', bases=(xbasis,ybasis))
    B_xy = dist.Field(name='B_xy', bases=(xbasis,ybasis))

    B['g'][0,0]=B_xx

    problem = de3.IVP([B_xx, B_xy, vp, tau_p, p], namespace=locals())

This prevents me from using B to create a force density ~div(B) in my momentum equation. I could go back to component form for everything but that just felt like moving backwards!

To summarize, I have a three questions:

1. Did I implement the upper-convected derivative correctly?

2. What are some options for constraining B?

3. Is it possible to define substitution rules for tensors?

Any thoughts would be appreciated. Thanks!

Some additional info from troubleshooting that might help narrow things down:

1a. If I decompose grad(vp) and consider only the antisymmetric part (rotations), then J=1 is satisfied.

**1b.If I consider symmetric and antisymmetric parts of the flow gradient and make the symmetric part tracless (E- (1/2)Tr(E), where Tr(E)=div(vp)), the problem returns.

2. Bxy=Byx is always satisfied even though I haven't explicitly enforced this symmetry.

 

[Jeffrey S. Oishi]

Hi Mike,

I haven't moved my viscoelastic scripts to d3 yet, but I will be doing it in the next few weeks (teaching makes for hard deadlines...).

It's a little hard to tell from just this snippet, but I don't see anything obvious. Your upper convected derivative looks fine. 

Are you doing Navier Stokes with very low Re or a pure stokes type equation? Have you tried forcing with some VP fixed in time? 

As an aside, I wouldn't suggest using components in your equations as the ultimate solution. 

If you're willing to send a script I would be happy to see if I can offer some more help. In particular, by seeing the whole system you are using would help figure out the best way to maintain the constraints on B. 

Jeff

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/da74bc82-1bd8-4ae1-8e67-5508eb33800an%40googlegroups.com.