Hello,
Currently, i am facing an issue while working on a hyperelasticity problem.In a growth problem the deformation gradient is splited into elastic and growth term. I defined a deformation gradient driven problem. So, the growth term in X direction is changing of the deformation gradient. However, I noticed that when I increase the parameter of XX term of growth deformation gradient it fails if the geometry rotate. It does not follow the growth tensor on deformed shape but it still tries to deform in in global X direction even the geometry is rotated. Is there any idea to solve that problem ? There is an example I applied growth in X direction but it does not follow the deformed X… Thanks
I also added a minimal code:
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
Type ="T2"
load_step=400#Total load step will be used in N-R
load_start=0 #Initial load
load_stop=4#*Inner_Area #End load in traction
inc_load = np.linspace(load_start,load_stop,load_step+1)
E, nu = 200, 0.49
Element = "CG"
order_u=2
mesh = RectangleMesh(Point(0., 0.), Point(1, 0.1), 40,4)
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
V2 = VectorFunctionSpace(mesh, Element, order_u)
du = TrialFunction(V2)
v = TestFunction(V2)
u = Function(V2)
quad_degree=3
metadata = {"quadrature_degree": quad_degree, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
# Kinematics
d = 3
grad_u = grad(u)
grad_u3 = as_tensor([[grad_u[0, 0], grad_u[0, 1], 0.0 ], [grad_u[1, 0], grad_u[1, 1],0.0], [0.0, 0.0,0.0]])
l_fac = Expression(('pr'),degree=1,pr=0)
Y = Expression(('x[1]'),degree=1)
g1=l_fac*Y
I = Identity(d)
F_1 = I + grad_u3
F_g = as_tensor([[1+g1, 0.0,0.0], [0.0,1.0,0.0], [0.0, 0.0,1.0]])
F_= inv(F_g)*F_1
J_g=det(F_g)
C_ = F_.T*F_
I1 = tr(C_)
I3 = det(C_)
tol =1e-14
class sliding_boundary_left(SubDomain): #Region definition to be fixed
def inside(self,x,on_boundary):
return on_boundary and near(x[0],0,tol)
class CornerPoint(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0., tol) and near(x[1], 0., tol)
left_bc = sliding_boundary_left()
corner_bc = CornerPoint()
boundary_part = MeshFunction("size_t", mesh, 1)
boundary_part.set_all(0)
bc1 = DirichletBC(V2.sub(0), Constant((0)), left_bc)
bc3 = DirichletBC(V2, Constant((0, 0)), corner_bc,method='pointwise')
bc = [bc1,bc3]
psi_iso = (lmbda/4)*(I3-1)-((lmbda/2)+mu)*ln((I3)**(0.5))+(mu/2)*(I1-3)
psi = psi_iso*J_g
# Total potential energy
Pi = psi*dxm
F = derivative(Pi, u, v)
# Compute Jacobian of F
J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, bc, J=J)
solver = NonlinearVariationalSolver(problem)
file_results = XDMFFile("Result.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
for t in inc_load:
l_fac.pr=t
print("t is",t)
solver.solve()
file_results.write(u, t)I think there may be a mix-up here in the convention for ordering indices of deformation gradients. F_1 and F_g are defined assuming the (standard) ordering
F_{iJ} = \frac{\partial x_i}{\partial X_J}\text{ ,}
where \mathbf{x} and \mathbf{X} are coordinates in the spatial and reference configurations (i.e., the convention of UFL’s grad operator). However, the product inv(F_g)*F_1 appears to assume a transposed definition (i.e., the nabla_grad convention). As a sanity check, consider the case where the intermediate configuration is compatible after growth, with coordinates \hat{\mathbf{X}}, and you want to use it as the new reference configuration for hyperelasticity. Then what you would want to feed into the hyperelastic model (corresponding to F_) is
\hat{F}_{i\alpha} = \frac{\partial x_i}{\partial\hat{X}_\alpha} = \frac{\partial x_i}{\partial X_K}\frac{\partial X_K}{\partial \hat{X}_\alpha} = F_{iK}(\mathbf{F}^{-1}_g)_{K\alpha}\text{ ,}
where I used lowercase indices on the spatial configuration, uppercase indices on the reference configuration, and Greek-letter indices on the intermediate configuration, for clarity. The product inv(F_g)*F_1 would be \hat{F}_{IJ} = (\mathbf{F}_g^{-1})_{I\alpha}F_{\alpha J}, which you can immediately spot as incorrect, since it forces you to break the conventions for what type of letter to use for indices on different configurations.
Hello kamensky,
Thanks a lot, you are totally right ! Last days I really struggled to find the mistake but I did not think that is about the convention. Now, it works.
Regards
hi, any chance you can upload the correct code with right adjustment ?