Element Orientation (Hyperelasticity)

variational formulation
1

Hi
I am working on a hyperelasticity problem.The material model is same as the hyperelasticity demo.
Now I want to add a new term to the strain energy function which is : cos(\theta) .tr(C)
Where C is the Right Cauchy-Green tensor and \theta is defined in this figure:

333

In my implementation the strain energy function is defined as:

psi = ((lmbda / 2.)*pow(ln(J),2)-mu*ln(J)+0.5*mu*(Ic-3))*x[0]*dx - inner(-p*n, u)*x[0]*ds(3) + cos(theta) * Ic * dx

The last term is what I want to add to the \psi. As you can see, \theta is calculated with respect to the center.
I want to add the term cos(\theta) .tr(C) to the strain energy function in a way that it is calculated with the orientation of the elements. For example I want \theta to be zero at the bottom edge, 45 degree at the middle of the geometry and 90 degree at the left edge, etc. I mean it is supposed to change with orientation of the elements (With respect to the center).
What is the best way to do it in FEniCS?

Here is my minimal work:
Note: The geometry and problem definition has been taken from Numerical tours of continuum mechanics

from dolfin import *
from mshr import *

Re = 11.
Ri = 9.

rect = Rectangle(Point(0., 0.), Point(Re, Re))
domain = Circle(Point(0., 0.), Re, 100) - Circle(Point(0., 0.), Ri, 100)
domain = domain - Rectangle(Point(0., -Re), Point(-Re, Re)) \
                - Rectangle(Point(0., 0.), Point(Re, -Re))

mesh = generate_mesh(domain, 10)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0) and on_boundary
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0) and on_boundary
class Outer(SubDomain):
    def inside(self, x, on_boundary):
        return near(sqrt(x[0]**2+x[1]**2), Re, 1e-1) and on_boundary

facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Bottom().mark(facets, 1)
Left().mark(facets, 2)
Outer().mark(facets, 3)
ds = Measure("ds", subdomain_data=facets)

x = SpatialCoordinate(mesh)

def GRAD(v):

    return  as_tensor([[v[0].dx(0), 0, v[0].dx(1)],
                        [0, v[0]/x[0], 0],
                        [v[1].dx(0), 0, v[1].dx(1)]])
E = Constant(1e6)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

n = FacetNormal(mesh)
p = Constant(10.)

V = VectorFunctionSpace(mesh, 'CG', degree=1)
du  = TrialFunction(V)           # Trial function
v = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement field

F = Identity(3) + GRAD(u)
C = F.T*F                   # Right Cauchy-Green tensor
Ic = tr(C)
J  = det(F)

#How should the last term (theta) be defined?
psi = ((lmbda / 2.)*pow(ln(J),2)-mu*ln(J)+0.5*mu*(Ic-3))*x[0]*dx - inner(-p*n, u)*x[0]*ds(3) + cos(theta) * Ic * dx

bcs = [DirichletBC(V.sub(1), Constant(0), facets, 1),
       DirichletBC(V.sub(0), Constant(0), facets, 2)]

F1 = derivative(psi, u, v)
Jac = derivative(F1, u, du)
problem = NonlinearVariationalProblem(F1, u, bcs, Jac)
solver = NonlinearVariationalSolver(problem)
solver.solve()
2

Cant you use the fact that cos(theta)=x[0]/sqrt(x[0]**2+x[1]**2) and directly insert this in your variational form?

1 Like
3

Thanks for the idea. It seems to work in this case that was a simplified form of my problem. In my real problem I want to calculate an integral that the \theta appears on the integral bounds:

c

The \theta in the integral bounds of the above integral should change with the element orientation. By using a trapezoidal rule, here is how this may be implemented into the variational form:

def trapezoidal(f, a, b, n):

    h = float(b - a) / n
    s = 0.0
    s += f(a)/2.0
    for i in range(1, n):
        s += f(a + i*h)
    s += f(b)/2.0
    return s * h

def G(omega):

    return 1

# How to define theta in integral boundas? (that changes with element orientation with respect to the center)
new_term = trapezoidal(G, theta - np.pi / 2., theta + np.pi / 2., 10)

psi = ... + new_term * Ic * dx

However I was not able to figure out how theta should be defined based on the first figure I provided above. Do you know how \theta in the integral bounds may be defined in the variational form?

4

I want to bring it up again in another way. I want to calculate this integral:
1
Where \theta is defined in the figure above. The \theta could be defined as:
4
So the integral at the bottom boundary where \theta=0(For example coordinate (9,0)) is:
2
On the left boundary where \theta=90 (For example coordinate (0,9)) :
3
In other words the integral has a different value at each coordinate as \theta changes with respect to the center (Please see the figure). I am using trapezoidal rule for integral calculation. In addition I tried to define the \theta as: acos(x[0]/sqrt(x[0]**2+x[1]**2))
Here is my minimal work:

from dolfin import *
from mshr import *
import numpy as np

Re = 11.
Ri = 9.

rect = Rectangle(Point(0., 0.), Point(Re, Re))
domain = Circle(Point(0., 0.), Re, 100) - Circle(Point(0., 0.), Ri, 100)
domain = domain - Rectangle(Point(0., -Re), Point(-Re, Re)) \
                - Rectangle(Point(0., 0.), Point(Re, -Re))

mesh = generate_mesh(domain, 40)

# Integral calculation using Trapezoidal rule
def trapezoidal(f, lower_limit, upper_limit, intervals):

    h = (upper_limit - lower_limit) / intervals
    s = 0.0
    s += f(lower_limit)/2.0
    for i in range(1, intervals):
        s += f(lower_limit + i*h)
    s += f(upper_limit)/2.0
    return s * h

V = FunctionSpace(mesh, "CG", 1)

x = SpatialCoordinate(mesh)

def cos(alpha):

    return np.cos(alpha)

# If theta = 0 everywhere
#A = trapezoidal(cos, -np.pi/2. ,np.pi/2., 100)

# variable theta based on the coordinate
A = trapezoidal(cos, acos(x[0]/sqrt(x[0]**2+x[1]**2)) - np.pi/2. , acos(x[0]/sqrt(x[0]**2+x[1]**2)) + np.pi/2., 100)

B = project(A, V)

# Theta = 0
print (B(9,0))

# Theta = 90
print (B(0,9))

The above code, returns this error:

AttributeError: 'Sum' object has no attribute 'cos'

Is there any way to resolve this issue?