diff --git a/src/Exports.jl b/src/Exports.jl
index 68d15f3..78170f3 100644
--- a/src/Exports.jl
+++ b/src/Exports.jl
@@ -68,6 +68,7 @@ end
 @publish PhysicalModels MagnetoMechModel
 @publish PhysicalModels ARAP2D
 @publish PhysicalModels ARAP2D_regularized
+@publish PhysicalModels NonlinearARAP2D
 @publish PhysicalModels HessianRegularization
 @publish PhysicalModels Hessian∇JRegularization
 @publish PhysicalModels ViscousIncompressible
diff --git a/src/PhysicalModels/MechanicalModels.jl b/src/PhysicalModels/MechanicalModels.jl
index dc6aeac..6cac6fc 100644
--- a/src/PhysicalModels/MechanicalModels.jl
+++ b/src/PhysicalModels/MechanicalModels.jl
@@ -1004,6 +1004,48 @@ struct ARAP2D <: IsoElastic
 end
 
 
+
+struct NonlinearARAP2D <: IsoElastic
+  μ::Float64
+  p::Float64
+  ρ::Float64
+  function NonlinearARAP2D(; μ::Float64, p::Float64, ρ::Float64=0.0)
+    new(μ, p, ρ)
+  end
+
+  function (obj::NonlinearARAP2D)(Λ::Float64=1.0)
+    μ = obj.μ
+    p = obj.p
+  
+    J(F) = det(F)
+    H(F) = det(F) * inv(F)'
+    f(F) = 0.5 * J(F)^(-1) * (tr((F)' * F)) - 1.0
+    g(x) = x^p
+    Ψ(F) = (μ/p) * g(f(F))
+
+    ∂f_∂F(F)   = F * J(F)^(-1)
+    ∂f_∂J(F)   = -1.0 / 2.0 * (tr((F)' * F)) * J(F)^(-2)
+    ∂2f_∂J2(F) = J(F)^(-3) * (tr((F)' * F))
+    ∂2f_∂FJ(F) = -J(F)^(-2) * F
+    ∂2f_∂FF(F) = J(F)^(-1) * I4
+
+    ∂g_∂x(x)   = p * x^(p-1)
+    ∂2g_∂x2(x) = p * (p-1) * x^(p-2)
+
+    ∂fu(F)  = ∂f_∂F(F) + ∂f_∂J(F) * H(F)
+    ∂fuu(F) = ∂2f_∂FF(F) + ∂2f_∂J2(F) * (H(F) ⊗ H(F)) + ∂2f_∂FJ(F) ⊗ H(F) + H(F) ⊗ ∂2f_∂FJ(F) + ∂f_∂J(F) * _∂H∂F_2D()
+
+    ∂Ψu(F)   =  (μ/p) * (∂g_∂x(f(F))* ∂fu(F))
+    ∂Ψuu(F)  =  (μ/p) * (∂g_∂x(f(F))* ∂fuu(F) + ∂2g_∂x2(f(F)) * (∂fu(F) ⊗ ∂fu(F)))
+
+    return (Ψ, ∂Ψu, ∂Ψuu)
+  end
+end
+
+
+
+
+
 struct IsochoricNeoHookean3D <: IsoElastic
   μ::Float64
 end
diff --git a/src/PhysicalModels/PhysicalModels.jl b/src/PhysicalModels/PhysicalModels.jl
index d23847e..31eb4f7 100644
--- a/src/PhysicalModels/PhysicalModels.jl
+++ b/src/PhysicalModels/PhysicalModels.jl
@@ -24,6 +24,7 @@ export IncompressibleNeoHookean2D_CV
 export IncompressibleNeoHookean3D_2dP
 export ARAP2D
 export ARAP2D_regularized
+export NonlinearARAP2D
 export VolumetricEnergy
 export MooneyRivlin3D
 export MooneyRivlin2D
diff --git a/test/TestConstitutiveModels/PhysicalModelTests.jl b/test/TestConstitutiveModels/PhysicalModelTests.jl
index 10e2b11..0f428b8 100644
--- a/test/TestConstitutiveModels/PhysicalModelTests.jl
+++ b/test/TestConstitutiveModels/PhysicalModelTests.jl
@@ -956,6 +956,13 @@ end
   test_equilibrium_at_rest_2D(model)
 end
 
+@testset "NonlinearARAP2D" begin
+  #  Memory estimate: 0 bytes, allocs estimate: 0.
+  model = NonlinearARAP2D(μ=μParams[1], p=3.0)
+  test_derivatives_2D_(model, Kinematics(Mechano, Solid), rtol=1e-13)
+  test_equilibrium_at_rest_2D(model)
+end
+
 
 @testset "HessianRegularization" begin
   # 3.56 μs      Histogram: log(frequency) by time      9.21 μs <