Added analytical model for hydrostatic flow dependence

Author: AHartmaier

src/pylabfea/basic.py

@@ -27,6 +27,16 @@
 """Tolerance: Plastic yielding if yield function > yf_tolerance"""
 
 
+def sig_hydro(sig):
+    if isinstance(sig, list):
+        sig = np.array(sig)
+    dim = len(sig.shape)  # determine whether sig is single stress tensor or array of tensors
+    if dim == 1:
+        return np.sum(sig[0:3]) / 3.
+    else:
+        return np.sum(sig[:, 0:3], axis=1) / 3.
+
+
 def sig_eq_j2(sig: np.ndarray):
     """Calculate sj2 equivalent stress from any stress tensor
 

src/pylabfea/data.py

@@ -460,7 +460,7 @@ def __init__(self, source, path_data='./',
             self.lc_data = source
             self.parse_data(epl_crit, epl_start, epl_max, depl)  # add data to mat_data
         elif isinstance(source, list) or isinstance(source, np.ndarray):
-            print('WARNING: This data type will be no longer supported.')
+            print('WARNING in initialization: Data of type list or ndarray will be no longer supported in future versions.')
             raw_data = np.array(source)
             self.convert_data(raw_data)  # add data to mat_data
         else:
@@ -526,7 +526,7 @@ def read_data(self, data_File: str):
                 except KeyError:
                     print("No texture_index found in this Data_Base.json -> Assign default value of 0")
                 if not self.mat_data['tx_data']:
-                    warnings.warn(f"WARNING: tx_data was set to false. I will just include qualitative texture info.")
+                    warnings.warn(f"WARNING in 'read_data': tx_data was set to false. I will just include qualitative texture info.")
                 else:
                     if 'GSH' in self.mat_data['tx_descriptor']:
                         # JS: Use GSH coefficients
@@ -614,7 +614,7 @@ def read_data(self, data_File: str):
                                              f"Data must be provided either im MPa or in GPa.")
                     else:
                         sfct = 1.
-                        print('Warning: No units for stresses are given. Assuming MPa.')
+                        print('Warning in "read_data": No units for stresses are given. Assuming MPa.')
                     Original_Stresses *= sfct
                     seq_full = FE.sig_eq_j2(Original_Stresses)
                     tens = [1]*6

src/pylabfea/material.py

@@ -643,7 +643,7 @@ def calc_seq(self, sig):
                 d0 = np.zeros(3)
             else:
                 hp = self.hill
-                if self.lhs:
+                if self.lhs is not None:
                     d0 = self.lhs
                 else:
                     d0 = np.ones(3) * self.drucker
@@ -658,7 +658,7 @@ def calc_seq(self, sig):
                      6. * hp[3] * np.square(sig[:, 3]) + \
                      6. * hp[4] * np.square(sig[:, 4]) + \
                      6. * hp[5] * np.square(sig[:, 5])
-                I2 *= 0.5
+                I2 *= 0.5  # corresponds to 3 * second invariant of stress deviator
                 self.msg['equiv'] = '6-parameter Hill, full Voigt stress'
                 # print('Full stress', np.sqrt(I2))
             else:
@@ -699,6 +699,10 @@ def calc_seqB(self, sv):
               np.abs(Stp1[1] - Stp2[0]) ** a + np.abs(Stp1[1] - Stp2[1]) ** a + np.abs(Stp1[1] - Stp2[2]) ** a + \
               np.abs(Stp1[2] - Stp2[0]) ** a + np.abs(Stp1[2] - Stp2[1]) ** a + np.abs(Stp1[2] - Stp2[2]) ** a
         seq = (0.25 * seq) ** (1. / a)
+        if self.lhs is not None:
+            d0 = self.lhs
+            I1 = (sv[0] * d0[0] + sv[1] * d0[1] + sv[2] * d0[2]) / 3.  # hydrostatic stress as 1st invariant
+            seq += I1
         return seq
 
     def calc_fgrad(self, sig, epl=None, seq=None,
@@ -827,7 +831,7 @@ def Jac(sig):
             h1 = self.hill[1]
             h2 = self.hill[2]
             if self.lhs is not None:
-                d3 = self.lhs
+                d3 = self.lhs / 3.
             else:
                 d3 = np.ones(3) * self.drucker / 3.
             if seq is None:
@@ -2523,9 +2527,6 @@ def plasticity(self, sy=None, sdim=6, drucker=0., khard=0.,
             self.sdim = sdim
         if hill is None and rv is None:
             hill = np.ones(self.sdim)
-            if lhs is not None:
-                raise ValueError('LHS parameters for anisotropic yield asymmetry provided, '
-                                 'but no anisotropy parameters for plastic yielding have been given.')
         elif hill is None:
             # compute Hill parameters from rv values
             hill = np.ones(self.sdim)
@@ -2592,12 +2593,17 @@ def plasticity(self, sy=None, sdim=6, drucker=0., khard=0.,
             self.barlat_exp = barlat_exp
         else:
             self.barlat = False
+        #if lhs is not None:
+        #    print(f'Using LHS parameters: {lhs}')
+        #    print(f'Using LHS with hill parameters: {hill}')
+        #    if self.barlat:
+        #        print(f'Using LHS with Barlat parameters: {barlat}')
 
     def from_data(self, param):
         """Define material properties from data sets generated in module `Data`:
         contains data on elastic and plastic behavior, including work hardening,
         for different crystallographic textures. Possible to extend to grain sizes,
-        grain shapes and porosities. Will invoke definition of elastic and plastic
+        grain shapes and porosities. Subroutine invokes definition of elastic and plastic
         parameters by calls to the methods `Material.elasticity` and `Material.plasticity`
         with the parameters provided in the data set.
         Also initializes current texture to first one in list and resets work hardening