How can I use a fourth-order strain energy?

variational formulation
1

Hello,
I am writing a bistable material model and am trying to use a 4th order polynomial as the strain energy. I am also open to using piecewise function with conditional() or heaviside functions. The following code (a single time step) runs with a quadratic strain energy but not the higher order polynomial. The error I get is:
UFLException: Found Argument in Power id=4830777408, this is an invalid expression.
This comes from the solve() line, and I get a similar error trying to use a piecewise strain-energy with conditional()

I have applied an initial strain with interpolate, but the strain should remain in the first ā€˜wellā€™ of the fourth-order polynomial

from dolfin import *
from ufl import replace
dt=1/50
beta = 1e3
L, H = 2, 1
mesh = RectangleMesh(Point(0., 0.), Point(L, H), 40, 20)
W = VectorFunctionSpace(mesh, 'Lagrange', degree=2)

Here I define the initial strain to test if it will return to the minimum at strain = 0

u_0 = Expression(("0.08*x[0]", "0.08*x[1]") ,degree=1)
u = interpolate(u_0,W)
u_old = interpolate(u_0,W)
u_ = TestFunction(W)
du = TrialFunction(W)

def eps(disp):
    return sym(grad(disp))

I have tried to define the strain energy as 3rd order or 4th order and both return the same error. If I try to define a piecewise energy using conditional(le(eps_vol, root), 0.5*b0*eps_vol**2, 0.5*b1*(eps_vol-eps_s)**2+phi_s) I get a similar error, but with ā€œConditional idā€ instead of Power.

def strain_energy(eps_vol):
    return 1296*pow(eps_vol,4)-633*pow(eps_vol,3)+80*pow(eps_vol,2)
    #return pow(eps_vol,2)

def dissipation_potential(u, u_old):
    return 0.5*beta*inner(eps(u)-eps(u_old),eps(u)-eps(u_old))

incremental_potential =  strain_energy(tr(eps(u))/2)*dx + dt*dissipation_potential(u,u_old)*dx

F = derivative(incremental_potential, u, u_)
form = replace(F, {u: du})

solve(lhs(form) == rhs(form), u)

I think the issue is with integrating lhs(form), but Iā€™m not totally sure.

2

Hello, your problem is non-linear when using a higher order strain energy, therefore you should not use the
solve(a == L, u) syntax, use instead:

F = derivative(incremental_potential, u, u_)
J = derivative(F, u, du)
solve(F==0, u, J=J)
3

Thank you, it worked!
Can I ask why you eliminated the line form = replace(F, {u: du}) ā€“ I had based my code on one of the examples but am still not totally sure why this line was necessary (and is no longer, in my case).

4

This is for linear problems when your residual F is linear i.e. of the form F(v)=a(u, v)-L(v) where v is a TestFunction and a a bilinear form. If you want to get the associated bilinear form you can do replace(F, {u, du}) which will give you a(du, v) with du a TrialFunction. But for non-linear problem you donā€™t have a bilinear form a in the first place.