added engineering constants feature for anisotropic material #216

src/Core/Solver/Verlet.jl

@@ -220,13 +220,34 @@ function compute_crititical_time_step(
         if mechanical
             bulk_modulus =
                 datamanager.get_property(iblock, "Material Model", "Bulk Modulus")
+            g_xy = datamanager.get_property(iblock, "Material Model", "Shear Modulus XY")
+            g_yz = datamanager.get_property(iblock, "Material Model", "Shear Modulus YZ")
+            g_zx = datamanager.get_property(iblock, "Material Model", "Shear Modulus ZX")
+            c_44 = datamanager.get_property(iblock, "Material Model", "C44")
+            c_55 = datamanager.get_property(iblock, "Material Model", "C55")
+            c_66 = datamanager.get_property(iblock, "Material Model", "C66")
             if !isnothing(bulk_modulus)
                 t = compute_mechanical_critical_time_step(
                     block_nodes[iblock],
                     datamanager,
                     bulk_modulus,
                 )
                 critical_time_step = test_timestep(t, critical_time_step)
+            elseif !isnothing(g_xy) && !isnothing(g_yz) && !isnothing(g_zx)
+                t = compute_mechanical_critical_time_step(
+                    block_nodes[iblock],
+                    datamanager,
+                    maximum([g_xy, g_yz, g_zx]),
+                )
+                critical_time_step = test_timestep(t, critical_time_step)
+                #TODO: temporary solution!!!
+            elseif !isnothing(c_44) && !isnothing(c_55) && !isnothing(c_66)
+                t = compute_mechanical_critical_time_step(
+                    block_nodes[iblock],
+                    datamanager,
+                    maximum([c_44 / 2, c_55 / 2, c_66 / 2]),
+                )
+                critical_time_step = test_timestep(t, critical_time_step)
             else
                 @error "No time step for material is determined because of missing properties."
                 return nothing

src/Models/Material/Material_Basis.jl

@@ -131,6 +131,36 @@ function get_all_elastic_moduli(
         nu = nu_fixed
         poissons = true
     end
+    if haskey(parameter, "Symmetry") &&
+       occursin("anisotropic", lowercase(parameter["Symmetry"]))
+        E_x = haskey(parameter, "Young's Modulus X")
+        E_y = haskey(parameter, "Young's Modulus Y")
+        E_z = haskey(parameter, "Young's Modulus Z")
+        nu_xy = haskey(parameter, "Poisson's Ratio XY")
+        nu_yz = haskey(parameter, "Poisson's Ratio YZ")
+        nu_zx = haskey(parameter, "Poisson's Ratio ZX")
+        g_xy = haskey(parameter, "Shear Modulus XY")
+        g_yz = haskey(parameter, "Shear Modulus YZ")
+        g_zx = haskey(parameter, "Shear Modulus ZX")
+        if E_x + E_y + E_z + nu_xy + nu_yz + nu_zx + g_xy + g_yz + g_zx > 0
+            if !E_x || !E_y || !E_z || !nu_xy || !nu_yz || !nu_zx || !g_xy || !g_yz || !g_zx
+                @error "Non isotropic material requires Young's Modulus X, Y, Z, Poisson's Ratio XY, YZ, ZX, Shear Modulus XY, YZ, ZX"
+                return nothing
+            end
+            parameter["Compute_Hook"] = true
+        else
+            for iID = 1:6
+                for jID = iID:6
+                    if !("C" * string(iID) * string(jID) in keys(parameter))
+                        @error "C" * string(iID) * string(jID) * " not defined"
+                        return nothing
+                    end
+                end
+            end
+            parameter["Compute_Hook"] = false
+        end
+        return
+    end
 
     # tbd non isotropic material check
     if bulk + youngs + shear + poissons < 2
@@ -193,21 +223,59 @@ Returns the Hooke matrix of the material.
 """
 function get_Hooke_matrix(parameter::Dict, symmetry::String, dof::Int64, ID::Int64 = 1)
     """https://www.efunda.com/formulae/solid_mechanics/mat_mechanics/hooke_plane_stress.cfm"""
-
-    if occursin("anisotropic", symmetry)
+    if occursin("anisotropic", lowercase(symmetry))
         aniso_matrix = get_MMatrix(36)
-        for iID = 1:6
-            for jID = iID:6
-                if "C" * string(iID) * string(jID) in keys(parameter)
+
+        if parameter["Compute_Hook"]
+
+            E_x = parameter["Young's Modulus X"]
+            E_y = parameter["Young's Modulus Y"]
+            E_z = parameter["Young's Modulus Z"]
+            nu_xy = parameter["Poisson's Ratio XY"]
+            nu_yz = parameter["Poisson's Ratio YZ"]
+            nu_zx = parameter["Poisson's Ratio ZX"]
+            g_xy = parameter["Shear Modulus XY"]
+            g_yz = parameter["Shear Modulus YZ"]
+            g_zx = parameter["Shear Modulus ZX"]
+
+            nu_yx = nu_xy * E_y / E_x
+            nu_xz = nu_zx * E_x / E_z
+            nu_zy = nu_yz * E_z / E_y
+
+            delta =
+                (
+                    1 - nu_xy * nu_yx - nu_yz * nu_zy - nu_zx * nu_xz -
+                    2 * nu_xy * nu_yz * nu_zx
+                ) / E_x *
+                E_y *
+                E_z
+
+            aniso_matrix[1, 1] = (1 - nu_yz * nu_zy) / E_y * E_z * delta
+            aniso_matrix[2, 2] = (1 - nu_zx * nu_xz) / E_z * E_x * delta
+            aniso_matrix[3, 3] = (1 - nu_xy * nu_yx) / E_x * E_y * delta
+
+            aniso_matrix[1, 2] = (nu_yx + nu_zx * nu_yz) / E_y * E_z * delta
+            aniso_matrix[2, 1] = (nu_xy + nu_xz * nu_zy) / E_z * E_x * delta
+
+            aniso_matrix[1, 3] = (nu_zx + nu_yx * nu_zy) / E_y * E_z * delta
+            aniso_matrix[3, 1] = (nu_xz + nu_xy * nu_yz) / E_x * E_y * delta
+
+            aniso_matrix[2, 3] = (nu_zy + nu_zx * nu_xy) / E_z * E_x * delta
+            aniso_matrix[3, 2] = (nu_yz + nu_xz * nu_yx) / E_x * E_y * delta
+
+            aniso_matrix[4, 4] = 2 * g_yz
+            aniso_matrix[5, 5] = 2 * g_zx
+            aniso_matrix[6, 6] = 2 * g_xy
+        else
+            for iID = 1:6
+                for jID = iID:6
                     value = parameter["C"*string(iID)*string(jID)]
-                else
-                    @error "C" * string(iID) * string(jID) * " not defined"
-                    return nothing
+                    aniso_matrix[iID, jID] = value
+                    aniso_matrix[jID, iID] = value
                 end
-                aniso_matrix[iID, jID] = value
-                aniso_matrix[jID, iID] = value
             end
         end
+
         if dof == 3
             return aniso_matrix
         elseif occursin("plane strain", symmetry)
@@ -218,12 +286,16 @@ function get_Hooke_matrix(parameter::Dict, symmetry::String, dof::Int64, ID::Int
             matrix[3, 3] = aniso_matrix[6, 6]
             return matrix
         elseif occursin("plane stress", symmetry)
-            matrix = get_MMatrix(36)
+            @info aniso_matrix
             inv_aniso = invert(aniso_matrix, "Hooke matrix not invertable")
+            @info inv_aniso
+            matrix = get_MMatrix(36)
+            @info aniso_matrix
             matrix[1:2, 1:2] = inv_aniso[1:2, 1:2]
             matrix[3, 1:2] = inv_aniso[6, 1:2]
             matrix[1:2, 3] = inv_aniso[1:2, 6]
             matrix[3, 3] = inv_aniso[6, 6]
+            @info matrix
             return invert(matrix, "Hooke matrix not invertable")
         else
             @error "2D model defintion is missing; plane stress or plane strain "

src/Support/Helpers.jl

@@ -112,15 +112,18 @@ function add_in_place!(
     return C
 end
 function get_MMatrix(len::Int64)
+    global A2x2
+    global A3x3
+    global A6x6
 
     if len == 4
-        global A2x2
+        A2x2 .= 0
         return A2x2
     elseif len == 9
-        global A3x3
+        A3x3 .= 0
         return A3x3
     elseif len == 36
-        global A6x6
+        A6x6 .= 0
         return A6x6
     else
         @error "MMatrix length $len not pre-allocated. Please add your size to helper.jl in get_MMatrix."

src/Support/Parameters/parameter_handling.jl

@@ -202,9 +202,24 @@ global expected_structure = Dict(
                                     "Bond Associated" => [Bool, false],
                                     "Bond Horizon" => [Union{Float64,Int64}, false],
                                     "Poisson's Ratio" => [Union{Float64,Int64}, false],
+                                    "Poisson's Ratio XY" =>
+                                        [Union{Float64,Int64}, false],
+                                    "Poisson's Ratio YZ" =>
+                                        [Union{Float64,Int64}, false],
+                                    "Poisson's Ratio ZX" =>
+                                        [Union{Float64,Int64}, false],
                                     "Young's Modulus" => [Union{Float64,Int64}, false],
+                                    "Young's Modulus X" =>
+                                        [Union{Float64,Int64}, false],
+                                    "Young's Modulus Y" =>
+                                        [Union{Float64,Int64}, false],
+                                    "Young's Modulus Z" =>
+                                        [Union{Float64,Int64}, false],
                                     "Bulk Modulus" => [Union{Float64,Int64}, false],
                                     "Shear Modulus" => [Union{Float64,Int64}, false],
+                                    "Shear Modulus XY" => [Union{Float64,Int64}, false],
+                                    "Shear Modulus YZ" => [Union{Float64,Int64}, false],
+                                    "Shear Modulus ZX" => [Union{Float64,Int64}, false],
                                     "Yield Stress" => [Union{Float64,Int64}, false],
                                     "Zero Energy Control" => [String, false],
                                     "C11" => [Union{Float64,Int64}, false],

test/fullscale_tests/test_aniso/Reference/aniso_eng.e.license

@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: 2023 Christian Willberg <[email protected]>, Jan-Timo Hesse <[email protected]>
+
+SPDX-License-Identifier: BSD-3-Clause

test/fullscale_tests/test_aniso/Reference/aniso_hooke.e.license

@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: 2023 Christian Willberg <[email protected]>, Jan-Timo Hesse <[email protected]>
+
+SPDX-License-Identifier: BSD-3-Clause

test/fullscale_tests/test_aniso/aniso_eng.yaml

@@ -0,0 +1,80 @@
+# SPDX-FileCopyrightText: 2023 Christian Willberg <[email protected]>, Jan-Timo Hesse <[email protected]>
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+PeriLab:
+  Discretization:
+    Node Sets:
+      Node Set 1: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
+    Type: "Text File"
+    Input Mesh File: "test_mesh.txt"
+  Models:
+    Material Models:
+      Mat_1:
+        Material Model: "Correspondence Elastic"
+        Symmetry: Anisotropic
+        Young's Modulus X: 2.5e3
+        Young's Modulus Y: 2.5e3
+        Young's Modulus Z: 2.5e3
+        Poisson's Ratio XY: 0.3
+        Poisson's Ratio YZ: 0.3
+        Poisson's Ratio ZX: 0.3
+        Shear Modulus XY: 2.15e3
+        Shear Modulus YZ: 2.15e3
+        Shear Modulus ZX: 2.15e3
+        Zero Energy Control: "Global"
+      Mat_2:
+        Material Model: "Correspondence Elastic"
+        Symmetry: Anisotropic
+        Young's Modulus X: 2.5e3
+        Young's Modulus Y: 2.5e3
+        Young's Modulus Z: 2.5e3
+        Poisson's Ratio XY: 0.3
+        Poisson's Ratio YZ: 0.3
+        Poisson's Ratio ZX: 0.3
+        Shear Modulus XY: 2.15e3
+        Shear Modulus YZ: 2.15e3
+        Shear Modulus ZX: 2.15e3
+        Zero Energy Control: "Global"
+  Blocks:
+    block_1:
+      Material Model: "Mat_1"
+      Density: 2.7e-9
+      Horizon: 2.1
+    block_2:
+      Material Model: "Mat_2"
+      Density: 2.7e-9
+      Horizon: 2.1
+  Boundary Conditions:
+    BC_1:
+      Variable: "Displacements"
+      Node Set: "Node Set 1"
+      Coordinate: "x"
+      Value: "0.1*x"
+      Type: Dirichlet
+    BC_2:
+      Variable: "Displacements"
+      Node Set: "Node Set 1"
+      Coordinate: "y"
+      Value: "0.0"
+      Type: Dirichlet
+    BC_3:
+      Variable: "Displacements"
+      Node Set: "Node Set 1"
+      Coordinate: "z"
+      Value: "0.0"
+      Type: Dirichlet
+  Solver:
+    Initial Time: 0.0
+    Final Time: 5.0e-6
+    Verlet:
+      Safety Factor: 1.0
+  Outputs:
+    Output1:
+      Output Filename: "aniso_eng"
+      Output File Type: Exodus
+      Output Frequency: 1
+      Output Variables:
+        Displacements: True
+        Forces: True
+        Volume: True