206
views
0
3 months ago by
hi there

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)
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.
Thanks.

Community: FEniCS Project

2
3 months ago by
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 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.

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.

written 3 months ago by Ovais
2
3 months ago by
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
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.
def quasistatic_solver(N, active_tension, t):

flags = ["-O3", "-ffast-math", "-march=native"]
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 = 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)
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

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
written 3 months ago by Emek