Strain energy function with additional integral term

variational formulation
1

Hi
I am solving a Hyperelasticity problem with a strain energy function in this form:

A2
A21124×199 7.24 KB

The \mu and \lambda are Lame constants. The last integral is an additional term. The term p(\theta) is defined as:

A1
where \alpha is a constant. The parameter \beta is defined as \beta=a_0 \cdot C \cdot a_0

C is the right Cauchy deformation tensor and a_0 is function of \theta:

A6

The parameter \beta could be rewritten:

A8
A81146×62 2.91 KB

The \beta in integrand (last term in the first equation) includes the components of C tensor . If I ignore the \beta and keep only p(\theta) in the integrand and use quad from scipy.integrate, it works fine but I need to have \beta included! In this case scipy.integrate does NOT work anymore.

Is there any idea how this integration should be implemented based on the above explanation?

Here is my implementation:

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

n = 1
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), n, n, n)

#define boundaries
class left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

class right(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return (on_boundary and (abs(x[0] - 1.0) < tol))

Left = left()
Right = right()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)

#Marking boundaries
Left.mark(boundaries, 1)
Right.mark(boundaries, 2)

domains = MeshFunction('size_t', mesh, mesh.topology().dim())
ds = Measure("ds", subdomain_data=boundaries)
dx = Measure('dx',subdomain_data=domains)

n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, 'CG', degree=1)

bc = DirichletBC(V, Constant((0., 0.,0.)),boundaries, 1)

du  = TrialFunction(V)           # Trial function
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement field

d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor
b = F * F.T                 # Left Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

E, nu = 1E6, 0.46
mu = 0.4e6
lmbda= Constant(E*nu/((1 + nu)*(1 - 2*nu)))

force  = Constant((1.0,  0.0, 0.0))
alpha = 1.26

# DefinE BETA
a_0 = as_vector([[1, 0,0]])
BETA = inner(dot(a_0,C),a_0)

def integrand(Theta):

    BETA = pow(np.cos(Theta),2)*C[0,0] + pow(np.sin(Theta),2)*C[1,1] + np.cos(Theta)*np.sin(Theta)*C[0,1] + np.cos(Theta)*np.sin(Theta)*C[1,0]

    return ((1. / (np.pi * alpha)) * np.exp(np.cos(2. * (Theta))))*BETA  # It works if BETA is removed

P_int = quad(integrand, - np.pi/2., np.pi/2.)
VAL = list(P_int)[0]
#print (VAL)

psi = ((lmbda / 2.)*pow(ln(J),2)+0.5*mu*(Ic-3) + VAL)*dx - dot(force, u)*ds(2)
F1 = derivative(psi, u, v)

# Compute Jacobian
Jac = derivative(F1, u, du)
problem = NonlinearVariationalProblem(F1, u, bc, Jac)

solver = NonlinearVariationalSolver(problem)
solver.solve()
2

You might try SciPy’s fixed_quad function instead of quad, and find a suitable quadrature degree empirically. The following quick test looks promising:

from dolfin import *
from scipy.integrate import fixed_quad

# Non-numerical function returning UFL, for testing:
def func(theta):
    return theta*Constant(1.0)

# For unclear reasons, fixed_quad returns a tuple, where the second component
# is always(?) None, but the first argument is what one expects, and does
# not try to cast the output of func or its own return value to type float.  
quad_deg = 2
int_func = fixed_quad(func,0,1,n=quad_deg)[0]

# The result remains a UFL type, not a numerical value:
print(int_func)
1 Like
3

Thank you. It was really promising and solved this issue.
I was wondering if there is a way to define the parameter \beta in terms of an exponential integral function. To be more specific I am trying to include Ei (\beta) in the strain energy function. The exponential integral function is available in scipy by calling: from scipy.special import expi
Here is an example:

from dolfin import *
from scipy.integrate import fixed_quad
from scipy.integrate import quad
import numpy as np
from scipy.special import expi

n = 1
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), n, n, n)

#define boundaries
class left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

class right(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return (on_boundary and (abs(x[0] - 1.0) < tol))

Left = left()
Right = right()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundaries.set_all(0)

#Marking boundaries
Left.mark(boundaries, 1)
Right.mark(boundaries, 2)

domains = MeshFunction('size_t', mesh, mesh.topology().dim())
ds = Measure("ds", subdomain_data=boundaries)
dx = Measure('dx',subdomain_data=domains)

n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, 'CG', degree=1)

bc = DirichletBC(V, Constant((0., 0.,0.)),boundaries, 1)

du  = TrialFunction(V)           # Trial function
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement field

d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor
b = F * F.T                 # Left Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

E, nu = 1E6, 0.46
mu = 0.4e6
lmbda= Constant(E*nu/((1 + nu)*(1 - 2*nu)))

force  = Constant((1.0,  0.0, 0.0))
alpha = 1.26

def func(Theta):
    BETA = expi(np.cos(Theta)*C[0,0])  # I have simplified BETA just as an example
    return ((1. / (np.pi * alpha)) * np.exp(np.cos(2. * (Theta))))*BETA

quad_deg = 2

int_func = fixed_quad(func,- np.pi/2., np.pi/2.,n=quad_deg)[0]

It returns: ufunc 'expi' not supported for the input types
It is because of the term C[0,0] in the definition of \beta which is an exponential integral function.
Do you know how this issue might be resolved?
Once thanks again for your help.

4

Some of the approximations of the exponential integral listed on the Wikipedia page could be coded directly in UFL, although I don’t know the details of their accuracy or range of validity.