Julia Implementation of Neo-Hookean Model

				
					using LinearAlgebra, Printf, Optim
function NH_uniax_incomp(mu, lam)
    return mu * (lam^2 - 1/lam)
end
function NH_uniax_comp(mu, kappa, lam1, lam2, wholeStressTensor)
    F = [lam1 0 0; 0 lam2 0; 0 0 lam2]
    J = det(F)
    B = F * transpose(F)
    Bstar = J^(-2/3) * B
    devBstar = Bstar - tr(Bstar)/3 * I
    S = mu/J * devBstar + kappa*(J-1)*I
    #@printf("mu=%g, kappa=%g, lam1=%g, lam2=%g -> S11=%g, S22=%g\n", mu, kappa, lam1, lam2, S[1,1], S[2,2])
    if (wholeStressTensor)
        return S
    else
        return abs(S[2,2])
    end
end
# parameters
mu = 0.99
kappa = 200
lam = 1.4
# incompressible
S1_incomp = NH_uniax_incomp(mu, lam)
@printf("S1_incomp=%.6g\n", S1_incomp)
# compressible
opt = optimize( L2 -> NH_uniax_comp(mu, kappa, lam, L2, false), 0.5/sqrt(lam), 1.5/sqrt(lam), GoldenSection())
Optim.minimum(opt)
lam2 = Optim.minimizer(opt)
S1_comp = NH_uniax_comp(mu, kappa, lam, lam2, true)
println("S1_comp=", S1_comp)