### how to transform grad(u)

206

views

0

hi there

i want to compute deformation gradient so that it always comes out to be of the following form :-

i want to compute deformation gradient so that it always comes out to be of the following form :-

I have been obtaining deformation gradient through F = I+grad(u) but can anyone kindly suggest a way to obtain this form of deformation gradient

```
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(5.0, 2.0, 5.0), 5, 2,5)
U = VectorFunctionSpace(mesh,"Lagrange", 2)
up = Function(W)
(u, p) = split(up)
f_0 = Constant((1, 0, 0))
s_0 = Constant((0, 1, 0))
n_0 = Constant((0, 0, 1))
d = u.geometric_dimension()
I = Identity(d)
X = grad(u)
F = I+(outer(f_0,f_0))*X[0,0] + (outer(s_0,s_0))*X[1,1] + (outer(n_0,n_0))*X[1,1]
```

I have tried through the above pasted code. Is the correct way to obtain what i want ? .. i am not sure because when i try to solve the nonlinear elasticity problem , the problem blows up. So i have asked this question just to be sure that i am computing deformation gradient in the correct manner.

Please do share your comments.

Thanks.

Community: FEniCS Project

### 2 Answers

2

1.) The deformation gradient tensor is diagonal only in basis composed of its eigenvectors.

2.) Even if you diagonalise it, there is no guarantee that two eigenvalues are the same (i.e. one eigenvalue `\lambda_t` with algebraic multiplicity 2) for whatever mechanical deformation the material undergoes.

3.) When you write tensor `\mathbb{F}` in a spectral decomposition form you assume, that

Do the classical `\mathbb{F} = \mathbb{I} + \text{grad}(\mathbf u)` computation and later diagonalise the tensor.

And the last thing, this question is too general - continuum mechanics related. Next time please post questions that cannot be answered on other solid mechanics forum, otherwise i'll close it.

2.) Even if you diagonalise it, there is no guarantee that two eigenvalues are the same (i.e. one eigenvalue `\lambda_t` with algebraic multiplicity 2) for whatever mechanical deformation the material undergoes.

3.) When you write tensor `\mathbb{F}` in a spectral decomposition form you assume, that

`f_0`

, `s_0`

and `n_0`

are its eigenvectors. But these eigenvectors must change from point to point and also as body deforms. The problem blows up because your assumptions are wrong.Do the classical `\mathbb{F} = \mathbb{I} + \text{grad}(\mathbf u)` computation and later diagonalise the tensor.

And the last thing, this question is too general - continuum mechanics related. Next time please post questions that cannot be answered on other solid mechanics forum, otherwise i'll close it.

2

Hi Ovais,

can you ask the question including your goal? What kind of a physical problem force you to that trial? I simply cannot understand the reason why you want to use anything else than F = Identity + grad(u)

Best, Emek

can you ask the question including your goal? What kind of a physical problem force you to that trial? I simply cannot understand the reason why you want to use anything else than F = Identity + grad(u)

Best, Emek

1

Thanks for your interest in my question emak.

I am trying to validate the FEM solution for a transversely isotropic cardiac tissue model with an analytical solution written in matlab or scipy. The geometry is cube geometry. I have applied stress my self and calculating stretches and pressure as output using a mixed function space formulation.

I realize , now , that the way deformation gradient was expressed in an earlier thesis (from which i took the picture pasted above) led me to believe that i should use some alternate form of deformation gradient.

But now i have been able to reproduce same result with F=I+grad(u) and using a vectorfunctionspace for fibers. I still have few queries.

If you can possibly share your email (because the moderator in his comments mentions this questions as 'not' a typical fenics question),i ll be able to send the source which i initially consulted. I am sharing the code that i am using to calculate stretches and pressure now.

I am trying to validate the FEM solution for a transversely isotropic cardiac tissue model with an analytical solution written in matlab or scipy. The geometry is cube geometry. I have applied stress my self and calculating stretches and pressure as output using a mixed function space formulation.

I realize , now , that the way deformation gradient was expressed in an earlier thesis (from which i took the picture pasted above) led me to believe that i should use some alternate form of deformation gradient.

But now i have been able to reproduce same result with F=I+grad(u) and using a vectorfunctionspace for fibers. I still have few queries.

If you can possibly share your email (because the moderator in his comments mentions this questions as 'not' a typical fenics question),i ll be able to send the source which i initially consulted. I am sharing the code that i am using to calculate stretches and pressure now.

```
def quasistatic_solver(N, active_tension, t):
flags = ["-O3", "-ffast-math", "-march=native"]
dolfin.parameters["form_compiler"]["quadrature_degree"] = 2
dolfin.parameters["form_compiler"]["representation"] = "uflacs"
dolfin.parameters["form_compiler"]["cpp_optimize"] = True
dolfin.parameters["form_compiler"]["cpp_optimize_flags"] = " ".join(flags)
mesh = UnitCubeMesh(N,N,N)
class X0_Plane(SubDomain):
def inside(self, x, on_boundry):
return (x[0] < 1e-4) and on_boundry
class X1_Plane(SubDomain):
def inside(self, x, on_boundry):
return (x[1] < 1e-4) and on_boundry
class X2_Plane(SubDomain):
def inside(self, x, on_boundry):
return (x[2] < 1e-4) and on_boundry
mf = FacetFunction("size_t", mesh)
mf.set_all(4)
x0_plane = X0_Plane(); x0_plane.mark(mf, 0)
x1_plane = X1_Plane(); x1_plane.mark(mf, 1)
x2_plane = X2_Plane(); x2_plane.mark(mf, 2)
#plot(mf,interactive=True)
# Space and functions
V = VectorFunctionSpace(mesh, "Lagrange", 2)
P = FunctionSpace(mesh,"Lagrange",1)
W = V*P
# Define functions
w = Function(W)
(u, p) = split(w)
(du, dp) = TestFunctions(W)
# Kinematics
I = Identity(3)
F = I + grad(u)
F = variable(F)
C = F.T*F
B = F*F.T
# Elasticity/material parameters
eta = 0
alpha = 0.2 # scaling parameter
a = 2280.0 *alpha
a_f = 1168.5 *alpha
b = 9.726
b_f = 15.779
# Unit vectors (needed to set up active stress)
FiberSpace = VectorFunctionSpace(mesh,'Quadrature',2)
f_0 = interpolate(Expression(('1','0','0')), FiberSpace) # Ref. config.
s_0 = interpolate(Expression(('0','1','0')), FiberSpace)
n_0 = interpolate(Expression(('0','0','1')), FiberSpace)
# Current config.
f = F*f_0
s = F*s_0
n = F*n_0
# Invariants
I_4f = inner(C*f_0,f_0)
I_1 = tr(C)
J = det(F)
# Initial conditions and boundary conditions
bc1 = DirichletBC(W.sub(0).sub(0), Constant(0), x0_plane)
bc2 = DirichletBC(W.sub(0).sub(1), Constant(0), x1_plane)
bc3 = DirichletBC(W.sub(0).sub(2), Constant(0), x2_plane)
bcs = [bc1,bc2,bc3]
T_a = Constant(0) # initial traction
## Constitutive Law ##
# Passive stress
cauchy_stress = a*exp(b*(I_1 - 3))*B + 2*a_f*(I_4f-1)*exp(b_f*pow(I_4f - 1, 2))*outer(f,f) - p*I
P_p = J*cauchy_stress*inv(F).T
# Active stress
active_cauchy_stress = T_a*(outer(f,f)+eta*outer(s,s)+eta*outer(n,n))
P_a = J*active_cauchy_stress*inv(F.T)
# Total stress
P = P_p + P_a
eq = inner(P,grad(du))*dx + inner(J-1,dp)*dx
Jac = derivative(eq, w)
lam_f = np.zeros(len(t))
lam_t = np.zeros(len(t))
pressure = np.zeros(len(t))
ufile = File('quasistatic_displacement.pvd')
for i in range(0,len(t)):
print 'Time-step:',i,'/',len(t)
T_a.assign(Constant(active_tension[i]))
solve(eq==0,w,bcs=bcs,J=Jac)
(u,p)=w.split(deepcopy=True)
ufile << u
lam_f[i] = u[0](np.array([1,0,0])) + 1
lam_t[i] = u[1](np.array([0,1,0])) + 1
pressure[i] = assemble(p*dx)/assemble(1.*dx(mesh))
return lam_f, lam_t, pressure
N = 5 # number of elements
t = np.linspace(0,1,100)
active_tension = 5000*np.exp(-100*(t-0.5)**2)
lam_f, lam_t, pressure = quasistatic_solver(N,active_tension,t)
```

written
3 months ago by
Ovais

Hi, I am happy that you brought your code to the track of continuum mechanics, it looks way more possible at the moment. You can find my e-mail address in my papers, have a look at in my web site: http://bilenemek.abali.org

By the way, Michal Habera has a valid point stating that you should ask about FEniCS and not about continuum mechanics. Your question was not clear at all. Believe me that no one is paid for helping you, FEniCS is not a commercial software.

Emek

By the way, Michal Habera has a valid point stating that you should ask about FEniCS and not about continuum mechanics. Your question was not clear at all. Believe me that no one is paid for helping you, FEniCS is not a commercial software.

Emek

written
3 months ago by
Emek

Please login to add an answer/comment or follow this question.

1. I am using f = F*f_0 so the fiber directions do change from point to point. I did not post the complete code because then this question would have become too much of continuum mechanics

2. We , the user are not developers , hence we expect 'fenics' to solve our solid mechanics issues. Which means some discussion related to solid/ fluid mechanics would essentially form part of the questions that are posted here.