Merge pull request #221 from adtzlr/add-solidbody-kirchhoff-stress

Add SolidBody kirchhoff stress

README.md

@@ -48,8 +48,11 @@ All notable changes to this project will be documented in this file. The format
 
 ## [Unreleased]
 
+### Added
+- Add `SolidBody.evaluate.kirchhoff_stress()` method. Contrary to the Cauchy stress method, this gives correct results in incompressible plane stress.
+
 ### Fixed
-- Fix `tovoigt()` helper for data with more or less than two trailing axes.
+- Fix `tovoigt()` helper for data with more or less than two trailing axes and 2D tensors.
 
 ## [3.1.0] - 2022-05-02
 

felupe/math/_tensor.py

@@ -279,8 +279,12 @@ def ddot(A, B, n=2, parallel=False):
 
 def tovoigt(A):
     "Convert (3, 3) tensor to (6, ) voigt notation."
-    B = np.zeros((6, *A.shape[2:]))
-    ij = [(0, 0), (1, 1), (2, 2), (0, 1), (1, 2), (0, 2)]
+    if A.shape[0] * A.shape[1] == 4:
+        B = np.zeros((3, *A.shape[2:]))
+        ij = [(0, 0), (1, 1), (0, 1)]
+    else:
+        B = np.zeros((6, *A.shape[2:]))
+        ij = [(0, 0), (1, 1), (2, 2), (0, 1), (1, 2), (0, 2)]
     for i6, (i, j) in enumerate(ij):
         B[i6] = A[i, j]
     return B

felupe/mechanics/_helpers.py

@@ -37,12 +37,13 @@ def __init__(self, vector, matrix):
 class Evaluate:
     "A class with evaluate methods of a SolidBody."
 
-    def __init__(self, gradient, hessian, cauchy_stress=None):
+    def __init__(self, gradient, hessian, cauchy_stress=None, kirchhoff_stress=None):
         self.gradient = gradient
         self.hessian = hessian
 
         if cauchy_stress is not None:
             self.cauchy_stress = cauchy_stress
+            self.kirchhoff_stress = kirchhoff_stress
 
 
 class Results:
@@ -51,6 +52,7 @@ class Results:
     def __init__(self, stress=False, elasticity=False):
 
         self.force = None
+        self._force = None
         self.stiffness = None
         self.kinematics = None
 

felupe/mechanics/_solidbody.py

@@ -28,7 +28,7 @@
 from .._field import Field, FieldMixed, FieldsMixed, FieldAxisymmetric
 from .._assembly import IntegralForm, IntegralFormMixed, IntegralFormAxisymmetric
 from ..constitution import AreaChange
-from ..math import inv, dot
+from ..math import inv, dot, transpose, det
 from ._helpers import Assemble, Evaluate, Results
 
 
@@ -54,6 +54,7 @@ def __init__(self, umat, field):
             gradient=self._gradient,
             hessian=self._hessian,
             cauchy_stress=self._cauchy_stress,
+            kirchhoff_stress=self._kirchhoff_stress,
         )
 
         self._area_change = AreaChange()
@@ -141,6 +142,19 @@ def _hessian(self, field=None, args=(), kwargs={}):
 
         return self.results.elasticity
 
+    def _kirchhoff_stress(self, field=None):
+
+        self._gradient(field)
+
+        if len(self.results.kinematics) > 1:
+            P = self.results.stress[0]
+        else:
+            P = self.results.stress
+
+        F = self.results.kinematics[0]
+
+        return dot(P, transpose(F))
+
     def _cauchy_stress(self, field=None):
 
         self._gradient(field)
@@ -150,6 +164,7 @@ def _cauchy_stress(self, field=None):
         else:
             P = self.results.stress
 
-        JiFT = self._area_change.function(self.results.kinematics[0])
+        F = self.results.kinematics[0]
+        J = det(F)
 
-        return dot(P, inv(JiFT))
+        return dot(P, transpose(F)) / J

tests/test_math.py

@@ -127,6 +127,7 @@ def test_math():
     fe.math.cdya(F, F)
 
     fe.math.tovoigt(C)
+    fe.math.tovoigt(C[:2, :2])
     fe.math.eigvals(C)
     fe.math.eigvals(C[:2, :2])
     fe.math.eigvals(C, shear=True)

tests/test_mechanics.py

@@ -46,13 +46,15 @@ def test_simple():
     F = b.results.kinematics[0]
     P = b.results.stress
     s = b.evaluate.cauchy_stress()
+    t = b.evaluate.kirchhoff_stress()
     C = b.results.elasticity
 
     assert K.shape == (81, 81)
     assert r.shape == (81, 1)
     assert F.shape == (3, 3, 8, 8)
     assert P.shape == (3, 3, 8, 8)
     assert s.shape == (3, 3, 8, 8)
+    assert t.shape == (3, 3, 8, 8)
     assert C.shape == (3, 3, 3, 3, 8, 8)
 
 
@@ -193,6 +195,10 @@ def test_solidbody():
             s2 = b.evaluate.cauchy_stress(u)
             assert np.allclose(s1, s2)
 
+            t1 = b.evaluate.kirchhoff_stress()
+            t2 = b.evaluate.kirchhoff_stress(u)
+            assert np.allclose(t1, t2)
+
 
 def test_solidbody_axi():
 
@@ -234,6 +240,10 @@ def test_solidbody_axi():
             s2 = b.evaluate.cauchy_stress(u)
             assert np.allclose(s1, s2)
 
+            t1 = b.evaluate.kirchhoff_stress()
+            t2 = b.evaluate.kirchhoff_stress(u)
+            assert np.allclose(t1, t2)
+
 
 def test_solidbody_mixed():
 
@@ -280,6 +290,10 @@ def test_solidbody_mixed():
             s2 = b.evaluate.cauchy_stress(u)
             assert np.allclose(s1, s2)
 
+            t1 = b.evaluate.kirchhoff_stress()
+            t2 = b.evaluate.kirchhoff_stress(u)
+            assert np.allclose(t1, t2)
+
 
 if __name__ == "__main__":
     test_simple()