Eigenvectors of 2D matrix field for Hyperelasticity problem

variational formulation
1

Hi
I solved a non linear variational problem as follows

U=VectorFunctionSpace(mesh,"P",1,2)
Unidfs=FunctionSpace(mesh,"P",1)
Matr=TensorFunctionSpace(mesh, "CG", 1)

def F(u): return grad(u) + Identity(2)
def C(u):return (F(u).T)*F(u)
def B(u): return F(u)*(F(u).T)
def lambda_1(u): return (tr(C(u)) + sqrt((tr(C(u))**2)-(4*det(C(u)))+ DOLFIN_EPS))/2
def lambda_2(u): return (tr(C(u)) - sqrt((tr(C(u))**2)-(4*det(C(u)))+ DOLFIN_EPS))/2

def Ogd_Energy(u):
    t1=(mu_1/alfa_1)*((lambda_1(u)**(alfa_1st))+(lambda_2(u)**(alfa_1st))-2)
    t2=(mu_2/alfa_2)*((lambda_1(u)**(alfa_2st))+(lambda_2(u)**(alfa_2st))-2)
    td=D_1*((det(C(u))-1)**2)
    return t1+t2+td

Energy=Ogd_Energy(u)*dx
var1=derivative(Energy,u,Test_u)
var2=derivative(var1,u,Trial_u)
bc_left=DirichletBC(U, Constant((0.0,0.0)), boundary_left)
bc_right=DirichletBC(U,Constant((uL,0.0)),boundary_right)
bc=[bc_left,bc_right]
problem = NonlinearVariationalProblem(var1,u,bc,var2)
solver = NonlinearVariationalSolver(problem)

solver.solve()
C_=C(u)
B_=B(u)
l1=lambda_1(u)
l2=lambda_2(u)
lF1=sqrt(l1)
lF2=sqrt(l2)

#eigenvectors?? second piola stress?

I obtained the displacement field “u” as I expected.
Now I have to find the eigenvectors of C_ and B_ and the stress tensors (Cauchy, first Piola, Second Piola). I tried to use the following codes:

1:

e1=interpolate(Constant((1,0)),U)
e2=interpolate(Constant((0,1)),U)
e11=project(outer(e1,e1),Matr)
e12=project(outer(e1,e2),Matr)
e21=project(outer(e2,e1),Matr)
e22=project(outer(e2,e2),Matr)

def eig_vect(Cm,l1,l2):
    l1v=project(l1,Unidfs)
    i=0
    imax=mesh.num_vertices()
    C=project(Cm,Matr)
    bf=project(inner(C,e11),Unidfs)
    cf=project(inner(C,e12),Unidfs)
    df=project(inner(C,e22),Unidfs)
    vers1=interpolate(Constant((1,0)),U)
    vers2=interpolate(Constant((0,1)),U)
    u1_f=Function(Unidfs)
    u2_f=Function(Unidfs)
    while i<imax +tollb:
        b=bf.vector()[i]
        c=cf.vector()[i]
        d=df.vector()[i]
        l1p=l1v.vector()[i]
        u1=conditional(le(abs(c),tollb),Constant(1.),abs(c)/sqrt(((l1p-b)**2)+(c**2)))
        u2=conditional(le(abs(c),tollb),Constant(0.0),abs(l1p-b)/sqrt(((l1p-b)**2)+(c**2)))
        u1_f.vector()[i]=u1
        u2_f.vector()[i]=u2
        i=i+1
    vect1=(u1_f*vers1)+(u2_f*vers2)
    vect2=(-u2_f*vers1)+(u1_f*vers2)
    return vect1,vect2
#first piola
def ti_ogd(a): return (mu_1*(a**(alfa_1-1)))+(mu_2*(a**(alfa_2-1)))+(2*D_1*lambda_1(u)*lambda_2(u)/a)
ein_vC1 , ein_vC2=eig_vect(C_,lambda1_C,lambda2_C)
ein_vB1 , ein_vB2=eig_vect(B_,lambda1_C,lambda2_C)
t1=ti_ogd(lF1)
t2=ti_ogd(lF2)
m1=outer(ein_vC1,ein_vB1)
m2=outer(ein_vC2,ein_vB2)
P=(t1*m1)+(t2*m2)#first piola
S=P*(inv(F(u)).T)
S=project(S,Matr)

but the function for eigenvectors is very slow and it isn’t very accurate. The second Piola stress is not as I expected it.
2 the second try is:

def eig_vect(Cm,l1,l2):
    C=project(Cm,Matr)
    b=C[0,0]
    c=C[0,1]
    d=C[1,1]
    u1=ufl.conditional(ufl.le(ufl.abs(c),tollb),Constant(1.),c/sqrt(((l1-b)**2)+(c**2)))
    u2=conditional(le(abs(c),tollb),Constant(0.0),(l1-b)/sqrt(((l1-b)**2)+(c**2)))
    vect1=as_vector([u1,u2])
    vect2=as_vector([-u2,u1])
    return vect1,vect2
#...

but i didn’t obtain a correct form for eigenvectors
3 in the third try, I find only the second Piola stress tensor:

def second_piola(T):
    psi=Ogd_Energy(u)
    C_=C(u)
    C_=variable(C_)
    S=diff(psi, C_)
    return S

but it get a constant zero matrix.
Excuse me for the length of the post and thank you for the answers.

2

In the third try, you would need to order the steps differently, e.g.,

C_ = variable(C(u))
psi = Ogd_Energy_C(C_)
S = 2*diff(psi,C_)

where Ogd_Energy_C is a function defining the energy density in terms of the deformation tensor instead of the displacement. (I also added a factor of two, assuming psi is the standard energy density with \mathbf{S} = \partial\psi/\partial\mathbf{E} = 2\partial\psi/\partial\mathbf{C}.)

To differentiate with respect to a variable using diff, the function being differentiated needs to be explicitly defined in terms of that variable, after it is registered with variable(). In the MWE, since C_ is defined after psi, psi is not considered to depend on C_, so diff(psi,C_) is zero. (UFL is not smart enough to figure out the implicit dependence of psi on C_, based on the separate definitions of psi and C_ in terms of u.)

1 Like