Merge pull request #70 from MultiSimOLab/anisomagneto

Implements anisotropic magneto-mechanical coupling

Author: jmartfrut

src/PhysicalModels/MagnetoMechanicalModels.jl

@@ -15,7 +15,7 @@ struct MagnetoMechModel{MG<:Magneto,MC<:Mechano} <: MagnetoMechano
     new{MG,MC}(magneto, mechano)
   end
 
-  function (obj::MagnetoMechModel)(Λ::Float64=1.0)
+  function (obj::MagnetoMechModel{<:Magneto,<:IsoElastic})(Λ::Float64=1.0)
     Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
     Ψmm, ∂Ψmm_u, ∂Ψmm_φ, ∂Ψmm_uu, ∂Ψmm_φu, ∂Ψmm_φφ = _getCoupling(obj.magneto, obj.mechano, Λ)
 
@@ -28,9 +28,26 @@ struct MagnetoMechModel{MG<:Magneto,MC<:Mechano} <: MagnetoMechano
 
     return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
   end
+
+    function (obj::MagnetoMechModel{<:Magneto,<:AnisoElastic})(Λ::Float64=1.0)
+    Ψm, ∂Ψm_u, ∂Ψm_uu = obj.mechano(Λ)
+    Ψmm, ∂Ψmm_u, ∂Ψmm_φ, ∂Ψmm_uu, ∂Ψmm_φu, ∂Ψmm_φφ = _getCoupling(obj.magneto, obj.mechano, Λ)
+
+    Ψ(F, ℋ₀, N) = Ψm(F,N) + Ψmm(F, ℋ₀, N)
+    ∂Ψu(F, ℋ₀, N) = ∂Ψm_u(F,N) + ∂Ψmm_u(F, ℋ₀, N)
+    ∂Ψφ(F, ℋ₀, N) = ∂Ψmm_φ(F, ℋ₀, N)
+    ∂Ψuu(F, ℋ₀, N) = ∂Ψm_uu(F,N) + ∂Ψmm_uu(F, ℋ₀, N)
+    ∂Ψφu(F, ℋ₀, N) = ∂Ψmm_φu(F, ℋ₀, N)
+    ∂Ψφφ(F, ℋ₀, N) = ∂Ψmm_φφ(F, ℋ₀, N)
+
+    return (Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ)
+  end
+
+
 end
 
 
+
 function _getCoupling(mag::Union{IdealMagnetic,IdealMagnetic2D}, ::Mechano, Λ::Float64)
   Ψmm, ∂Ψmm_∂u, ∂Ψmm_∂φ, ∂Ψmm_∂uu, ∂Ψmm_∂φu, ∂Ψmm_∂φφ = mag(Λ)
 

test/TestConstitutiveModels/PhysicalModelTests.jl

@@ -644,6 +644,47 @@ end
 end
 
 
+@testset "HardMagnetic_SoftMaterial3D_aniso" begin
+
+  ∇u = TensorValue(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0) * 1e-3
+  ∇φ = VectorValue(1.0, 2.0, 3.0)
+  N = VectorValue(0.0, 0.0, 1.0)
+
+  modelMRiso = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
+  modelMRaniso = MooneyRivlin3D(λ=3.0, μ1=1.0, μ2=2.0)
+  modelMR=modelMRiso+modelMRaniso
+  modelID = HardMagnetic(μ0=1.2566e-6, αr=40e-3, χe=0.0, χr=8.0;  βmok=1.0, βcoup=1.0)
+  modelmagneto=modelMR+modelID
+  Ψ, ∂Ψu, ∂Ψφ, ∂Ψuu, ∂Ψφu, ∂Ψφφ = modelmagneto()
+  K=Kinematics(Mechano,Solid)
+  F, _, _ = get_Kinematics(K)
+  Km=Kinematics(Magneto,Solid)
+  H0 = get_Kinematics(Km)
+
+
+  # ∂Ψu_(F) =TensorValue(ForwardDiff.gradient(x -> Ψ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F)))
+  # ∂Ψφ_(H) =VectorValue(ForwardDiff.gradient(x -> Ψ(get_array( F(∇u)), x,get_array(N) ), get_array(H)))
+  # ∂Ψuu_(F) =TensorValue(ForwardDiff.hessian(x -> Ψ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F)))
+  # ∂Ψφu_(H) =TensorValue(ForwardDiff.jacobian(x -> ∂Ψφ(x, get_array( H0(∇φ)),get_array(N) ), get_array(F(∇u))))
+  # ∂Ψφφ_(H) =TensorValue(ForwardDiff.jacobian(x -> ∂Ψφ(get_array( F(∇u)), x,get_array(N) ), get_array(H)))
+
+  # norm(∂Ψu_(F(∇u))) 
+  # norm(∂Ψφ_(H0(∇φ))) 
+  # norm(∂Ψuu_(F(∇u))) 
+  # norm(∂Ψφu_(H0(∇φ))) 
+  # norm(∂Ψφφ_(H0(∇φ))) 
+
+
+  @test Ψ(F(∇u), H0(∇φ), N) == 0.003187725760804994 
+  @test norm(∂Ψu(F(∇u), H0(∇φ), N)) == 0.4966662732306754
+  @test norm(∂Ψφ(F(∇u), H0(∇φ), N)) == 4.660348298920368e-6 
+  @test norm(∂Ψuu(F(∇u), H0(∇φ), N)) ==60.735729041745294 
+  @test norm(∂Ψφu(F(∇u), H0(∇φ), N)) == 1.2369035467980284e-5 
+  @test norm(∂Ψφφ(F(∇u), H0(∇φ), N)) ==  2.1878750641250348e-6
+end
+
+
+
 
 @testset "HardMagnetic_SoftMaterial3D" begin