Author: Jukka Aho
Email: jukka.aho@kapsi.fi
In [1]:
tic()
Out[1]:
We consider here linear quadrangle element and evaluate solution for it. It has been verified using Elmer that the correct solution for minimum displacement in tip is $u_\mathrm{min,nonlinear}$ = -2.22224475 for geometrically nonlinear case and $u_\mathrm{min,linear} = -2.18056991$ for linear case.
Solution domain is $\Omega = [0, 10] \times [0, 1] \in \mathbb{R}^2$ and we are having $F=-2 \mathrm{N}$ nodal force in $(10, 1)$. $E=90$, $\nu=0.25$. Left side is fixed, $u|_{x=0}=0$.
Geometry, material definitions and shape functions:
In [2]:
X = [[0.0 10.0 10.0 0.0],
[0.0 0.0 1.0 1.0]]
# Partial derivatives of bilinear Lagrange polynomials
dNdξ(ξ) = [[-(1-ξ[2])/4.0 -(1-ξ[1])/4.0],
[ (1-ξ[2])/4.0 -(1+ξ[1])/4.0],
[ (1+ξ[2])/4.0 (1+ξ[1])/4.0],
[-(1+ξ[2])/4.0 (1-ξ[1])/4.0]]
E = 90
ν = 0.25
μ = E/(2*(1+ν))
λ = E*ν/((1+ν)*(1-2*ν))
λ = 2*λ*μ/(λ + 2*μ)
μ, λ
Out[2]:
In [3]:
function f_int(X, u, dNdξ, λ, μ, dim=2)
T = zeros(size(X))
I = eye(dim)
function J(ξ)
Jᵀ = X*dNdξ(ξ)
∇N = inv(Jᵀ)*dNdξ(ξ)'
∇u = u*∇N'
F = I + ∇u # Deformation gradient
E = 1/2*(∇u' + ∇u + ∇u'*∇u) # Green-Lagrange strain tensor
P = λ*trace(E)*I + 2*μ*E # PK1 stress tensor
S = F*P # PK2 stress tensor
return S*∇N*det(Jᵀ)
end
a = 1/sqrt(3)
ipoints = [[-a -a], [a -a], [a a], [-a a]]
iweights = [1 1 1 1]
for m = 1:length(iweights)
w = iweights[m]
ξ = ipoints[m, :]
T += w*J(ξ)
end
return T
end
Out[3]:
In [4]:
using NLsolve
function solve()
free_dofs = [3, 4, 5, 6]
u = zeros(2, 4)
F = zeros(2, 4)
F[2, 3] = -2.0
function f(u_)
u[free_dofs] = u_
T = f_int(X, u, dNdξ, λ, μ)
R = T-F
return R[free_dofs]
end
# Go!
sol = nlsolve(not_in_place(f), zeros(4), store_trace=true, show_trace=true)
println(sol)
u[free_dofs] = sol.zero
return u
end
u = solve()
Out[4]:
Difference to Elmer solution is
In [5]:
abs(u[2, 3] - -2.22224475)
Out[5]:
Linear solution is achieved by assuming small deformation, i.e. :
F = I + ∇u ≈ I
E = 1/2*(∇u' + ∇u + ∇u'*∇u) ≈ 1/2*(∇u' + ∇u)OFFTOPIC: Here we should put something like
using JuliaFEM
test_allclose(u[2, 3], -2.22224457, tol = ...)
and get a nice green "PASSED" if everything is ok at this point, otherwise assertionerror
For more robust solution, we shoud get analytical Jacobian for nonlinear system of equations. If not given it's probably some sort of finite difference approximation, i.e.
In [6]:
function Lin(f, h=1.0e-6)
function D(x)
J_approx = zeros((8, 8))
for i=1:8
Δx = zeros(2, 4)
Δx[i] += h
Δf = f(x+Δx) - f(x)
J_approx[i, :] = Δf[:]/h
end
return J_approx
end
return D
end
J_approx = Lin(u -> f_int(X, u, dNdξ, λ, μ))
J_approx(u)
Out[6]:
A more elegant solution is to use automatic differentiation to get analytical tangent stiffness matrix
In [7]:
using ForwardDiff
function T!(u, T)
T[:] = f_int(X, reshape(u, 2, 4), dNdξ, λ, μ)
end
J_analytical = forwarddiff_jacobian(T!, Float64, fadtype=:dual, n=8, m=8)
J_analytical(reshape(u, 8))
Out[7]:
Difference is in smooth situations quite small, but analytical is of course faster.
In [8]:
norm(J_approx(u) - J_analytical(reshape(u, 8)))
Out[8]:
If we use this analytical Jacobian we should get quadratic convergence when using Newton-Rhapson method:
In [9]:
function solve2()
free_dofs = [3, 4, 5, 6]
u = zeros(2, 4)
F = zeros(2, 4)
F[2, 3] = -2.0
function f(u_)
u[free_dofs] = u_
T = f_int(X, u, dNdξ, λ, μ)
R = T-F
return R[free_dofs]
end
function g(u_)
u[free_dofs] = u_
J = J_analytical(reshape(u, 8))
return J[free_dofs, free_dofs]
end
# Go!
sol = nlsolve(not_in_place(f, g), zeros(4), method = :newton, store_trace=true, show_trace=true)
println(sol)
u[free_dofs] = sol.zero
return u
end
u = solve2()
Out[9]:
Clearly the last iterations are converging quadratic speeds. Own naive implementation of solver would be:
In [10]:
function solve3()
free_dofs = [3, 4, 5, 6]
u = zeros(2, 4)
F = zeros(2, 4)
F[2, 3] = -2.0
function f(u_)
u[free_dofs] = u_
T = f_int(X, u, dNdξ, λ, μ)
R = T-F
return R[free_dofs]
end
function g(u_)
u[free_dofs] = u_
J = J_analytical(reshape(u, 8))
return J[free_dofs, free_dofs]
end
u2 = zeros(4)
for i=1:7
du = g(u2) \ -f(u2)
u2 += du
println("norm = ",norm(du))
end
u[free_dofs] = u2
return u
end
u = solve3()
Out[10]:
I don't know what NLsolve does because it takes a bit longer to converge?
At end we do compare results between linear and nonlinear systems and put all together.
In [11]:
function f_int(X, u, dNdξ, λ, μ, dim=2; nlgeom=true)
T = zeros(size(X))
I = eye(dim)
function J(ξ)
Jᵀ = X*dNdξ(ξ)
∇N = inv(Jᵀ)*dNdξ(ξ)'
∇u = u*∇N'
F = I # Deformation gradient for small deformations
E = 1/2*(∇u' + ∇u) # Small strain tensor
if nlgeom # .. if nonlinear geometry
F += ∇u # .. add nonlinear part to deformation gradient, and
E += 1/2*(∇u'*∇u) # add nonlinear part to strain tensor => Green-Lagrange strain tensor
end
P = λ*trace(E)*I + 2*μ*E # PK1 stress tensor
S = F*P # PK2 stress tensor
return S*∇N*det(Jᵀ)
end
a = 1/sqrt(3)
ipoints = [[-a -a], [a -a], [a a], [-a a]]
iweights = [1 1 1 1]
for m = 1:length(iweights)
w = iweights[m]
ξ = ipoints[m, :]
T += w*J(ξ)
end
return T
end
function Lin(f, h=1.0e-6)
function D(x)
J_approx = zeros((8, 8))
for i=1:8
Δx = zeros(2, 4)
Δx[i] += h
Δf = f(x+Δx) - f(x)
J_approx[i, :] = Δf[:]/h
end
return J_approx
end
return D
end
function solve4!(F, u; nlgeom=true, analytical_jacobian=false)
free_dofs = [3, 4, 5, 6]
# Tangent stiffness matrix
function T!(u, T)
T[:] = f_int(X, reshape(u, 2, 4), dNdξ, λ, μ, nlgeom=nlgeom)
end
if analytical_jacobian
Jacobian = ForwardDiff.forwarddiff_jacobian(T!, Float64, fadtype=:dual, n=8, m=8)
else
Jacobian = Lin(u -> f_int(X, u, dNdξ, λ, μ))
end
# Assembly of residual vector R(u) = T(u) - F(u)
function f(u_)
u[free_dofs] = u_
T = f_int(X, u, dNdξ, λ, μ, nlgeom=nlgeom)
return (T-F)[free_dofs]
end
# Derivative of residual vector w.r.t displacement field u, "tangent stiffness"
function g(u_)
u[free_dofs] = u_
if analytical_jacobian
J = Jacobian(reshape(u, 8))
else
J = Jacobian(u)
end
return J[free_dofs, free_dofs]
end
# Go!
#sol = nlsolve(not_in_place(f, g), zeros(4), method=:newton) # not working as expected
#sol = nlsolve(not_in_place(f), zeros(4)) # working, but not quadratic convergence
# working, very well :)
u2 = zeros(4)
for i=1:100
du = g(u2) \ -f(u2)
u2 += du
if norm(du) < 1.0e-6
#println("norm = ",norm(du), " converged in ", i, " iterations")
break
end
end
#u0[free_dofs] = sol.zero
u0[free_dofs] = u2
return u0
end
u0 = zeros(2, 4)
F = zeros(2, 4)
F[2, 3] = -20.0
solve4!(F, u0; nlgeom=true, analytical_jacobian=true)
Out[11]:
In [12]:
F = zeros(2, 4)
N = 41
F[2, 3] = -20.0
uall = zeros(2, 4, 2, N)
u0 = zeros(2, 4)
t = linspace(0, 1, N)
for i = 1:N
@printf "Solving %2i F = %8.3f " i t[i]*F[2, 3]
u0 = solve4!(t[i]*F, u0; nlgeom=false, analytical_jacobian=false)
uall[:, :, 1, i] = u0
@printf "%8.3f" u0[2, 3]
u0 = solve4!(t[i]*F, u0; nlgeom=true, analytical_jacobian=true)
uall[:, :, 2, i] = u0
@printf "%8.3f" u0[2, 3]
println()
end
In [13]:
import PyPlot
In [14]:
PyPlot.plot(-reshape(uall[2, 3, 1, :], N), -t*F[2, 3], color="red", linewidth=1.0, label="linear")
PyPlot.plot(-reshape(uall[2, 3, 2, :], N), -t*F[2, 3], color="blue", linewidth=1.0, label="nonlinear")
PyPlot.legend(loc="best")
PyPlot.grid()
PyPlot.xlabel("Displacement [m]")
PyPlot.ylabel("Force [N]")
Out[14]:
Last remarks. For some unknown reason this notebooks takes occasionally a very long time to complete for 40 different solutions. Sometimes NLsolve takes a very long time when using Newton algorithm. Without Newton all passes well. Actually my implementation of Newton method does not have any problems with speed.
In [15]:
toc()
Out[15]: