From 9dc9d6735ad809dc1f1a53ce0e0159c80802ad42 Mon Sep 17 00:00:00 2001 From: HarryvL <35259498+HarryvL@users.noreply.github.com> Date: Wed, 14 Sep 2022 01:30:17 +0200 Subject: [PATCH] [FEM] added critical strain ratio functionality (#7467) * added critical strain ratio functionality --- src/Mod/Fem/App/FemVTKTools.cpp | 2 + src/Mod/Fem/femobjects/result_mechanical.py | 7 ++ src/Mod/Fem/femresult/resulttools.py | 111 ++++++++++++++------ 3 files changed, 90 insertions(+), 30 deletions(-) diff --git a/src/Mod/Fem/App/FemVTKTools.cpp b/src/Mod/Fem/App/FemVTKTools.cpp index a56fd4434b..dd36bec364 100644 --- a/src/Mod/Fem/App/FemVTKTools.cpp +++ b/src/Mod/Fem/App/FemVTKTools.cpp @@ -718,6 +718,8 @@ std::map _getFreeCADMechResultScalarProperties() { resFCScalProp["NodeStrainXZ"] = "Strain xz component"; resFCScalProp["NodeStrainYZ"] = "Strain yz component"; resFCScalProp["Peeq"] = "Equivalent Plastic Strain"; + resFCScalProp["CriticalStrainRatio"] = "Critical Strain Ratio"; + // the following three are filled in all cases // https://forum.freecadweb.org/viewtopic.php?f=18&t=33106&start=70#p296317 // it might be these can be generated in paraview from stress tensor values as diff --git a/src/Mod/Fem/femobjects/result_mechanical.py b/src/Mod/Fem/femobjects/result_mechanical.py index a0a26ebf4b..a3faaf4f4c 100644 --- a/src/Mod/Fem/femobjects/result_mechanical.py +++ b/src/Mod/Fem/femobjects/result_mechanical.py @@ -292,6 +292,13 @@ class ResultMechanical(base_fempythonobject.BaseFemPythonObject): "", True ) + obj.addProperty( + "App::PropertyFloatList", + "CriticalStrainRatio", + "NodeData", + "", + True + ) # initialize the Stats with the appropriate count of items # see fill_femresult_stats in femresult/resulttools.py diff --git a/src/Mod/Fem/femresult/resulttools.py b/src/Mod/Fem/femresult/resulttools.py index af79f0fc7c..ed2572f966 100644 --- a/src/Mod/Fem/femresult/resulttools.py +++ b/src/Mod/Fem/femresult/resulttools.py @@ -411,13 +411,67 @@ def add_principal_stress_std(res_obj): res_obj.PrincipalMin = prinstress3 res_obj.MaxShear = shearstress FreeCAD.Console.PrintLog("Added standard principal stresses and max shear values.\n") + + # + # Add critical strain ratio using the Stress Modified Critical Strain (SMCS) criterion + # Forum Discussion: https://forum.freecadweb.org/viewtopic.php?f=18&t=35893#p303392 + # Background: https://www.vtt.fi/inf/julkaisut/muut/2017/VTT-R-01177-17.pdf + # + # critical strain ratio = peeq / critical_strain (>1.0 indicates ductile rupture) + # peeq = equivalent plastic strain + # critical strain = alpha * np.exp(-beta * T) + # alpha and beta are material parameters, where alpha can be related to unixial test data (user input) and + # beta is normally kept fixed at 1.5, unless available from extensive research experiments + # T = pressure / von Mises stress (stress triaxiality) + # + + MatMechNon = FreeCAD.ActiveDocument.getObject('MaterialMechanicalNonlinear') + if MatMechNon: + stress_strain = MatMechNon.YieldPoints + if stress_strain: + i = -1 + while stress_strain[i] == "": i -= 1 + critical_uniaxial_strain = float(stress_strain[i].split(",")[1]) + alpha = np.sqrt(np.e) * critical_uniaxial_strain # stress triaxiality T = 1/3 for uniaxial test + beta = 1.5 + if res_obj.Peeq: + res_obj.CriticalStrainRatio = calculate_csr(prinstress1, prinstress2, prinstress3, alpha, beta, + res_obj) + return res_obj -def get_concrete_nodes(res_obj): - +def calculate_csr(ps1, ps2, ps3, alpha, beta, res_obj): # - # HarryvL: determine concrete / non-concrete nodes + # calculate critical strain ratio + # Forum Discussion: https://forum.freecadweb.org/viewtopic.php?f=18&t=35893#p303392 + # Background: https://www.vtt.fi/inf/julkaisut/muut/2017/VTT-R-01177-17.pdf + # + # critical strain ratio = peeq / critical_strain (>1.0 indicates ductile rupture) + # peeq = equivalent plastic strain + # critical strain = alpha * np.exp(-beta * T) + # alpha and beta are material parameters, where alpha can be related to unixial test data (user input) and + # beta is normally kept fixed at 1.5, unless available from extensive research experiments + # T = pressure / von Mises stress (stress triaxiality) + # + # + csr = [] # critical strain ratio + nsr = len(ps1) # number of stress results + for i in range(nsr): + p = (ps1[i] + ps2[i] + ps3[i]) / 3.0 # pressure + svm = np.sqrt(1.5 * (ps1[i] - p) ** 2 + 1.5 * (ps2[i] - p) ** 2 + 1.5 * ( + ps3[i] - p) ** 2) # von Mises stress: https://en.wikipedia.org/wiki/Von_Mises_yield_criterion + if svm != 0.: + T = p / svm # stress triaxiality + else: + T = 0. + critical_strain = alpha * np.exp(-beta * T) # critical strain + csr.append(abs(res_obj.Peeq[i]) / critical_strain) # critical strain ratio + return csr + +def get_concrete_nodes(res_obj): + # + # determine concrete / non-concrete nodes # from femmesh.meshtools import get_femnodes_by_refshape @@ -460,9 +514,8 @@ def get_concrete_nodes(res_obj): def add_principal_stress_reinforced(res_obj): - # - # HarryvL: determine concrete / non-concrete nodes + # determine concrete / non-concrete nodes # ic = get_concrete_nodes(res_obj) @@ -481,7 +534,7 @@ def add_principal_stress_reinforced(res_obj): ps2v = [] ps3v = [] # - # HarryvL: additional arrays to hold reinforcement ratios + # additional arrays to hold reinforcement ratios # and mohr coulomb stress # rhx = [] @@ -522,7 +575,7 @@ def add_principal_stress_reinforced(res_obj): if ic[isv] == 1: # - # HarryvL: for concrete scxx etc. are affected by + # for concrete scxx etc. are affected by # reinforcement (see calculate_rho(stress_tensor)). for all other # materials scxx etc. are the original stresses # @@ -558,7 +611,7 @@ def add_principal_stress_reinforced(res_obj): res_obj.PrincipalMin = prinstress3 res_obj.MaxShear = shearstress # - # HarryvL: additional concrete and principal stress plot + # additional concrete and principal stress plot # results for use in _ViewProviderFemResultMechanical # res_obj.ReinforcementRatio_x = rhx @@ -610,13 +663,12 @@ def calculate_von_mises(stress_tensor): normal = stress_tensor[:3] shear = stress_tensor[3:] pressure = np.average(normal) - return np.sqrt(1.5 * np.linalg.norm(normal - pressure)**2 + 3.0 * np.linalg.norm(shear)**2) + return np.sqrt(1.5 * np.linalg.norm(normal - pressure) ** 2 + 3.0 * np.linalg.norm(shear) ** 2) def calculate_principal_stress_std( stress_tensor ): - # if NaN is inside the array, which can happen on Calculix frd result files return NaN # https://forum.freecadweb.org/viewtopic.php?f=22&t=33911&start=10#p284229 # https://forum.freecadweb.org/viewtopic.php?f=18&t=32649#p274291 @@ -645,7 +697,7 @@ def calculate_principal_stress_std( def calculate_principal_stress_reinforced(stress_tensor): # - # HarryvL - calculate principal stress vectors and values + # - calculate principal stress vectors and values # - for total stresses use stress_tensor[0], stress_tensor[1], stress_tensor[2] # on the diagonal of the stress tensor # @@ -668,7 +720,7 @@ def calculate_principal_stress_reinforced(stress_tensor): eigenvalues, eigenvectors = np.linalg.eig(sigma) # - # HarryvL: suppress complex eigenvalue and vectors that may occur for + # suppress complex eigenvalue and vectors that may occur for # near-zero (numerical noise) stress fields # @@ -690,9 +742,8 @@ def calculate_principal_stress_reinforced(stress_tensor): def calculate_rho(stress_tensor, fy): - # - # HarryvL - Calculation of Reinforcement Ratios and + # Calculation of Reinforcement Ratios and # Concrete Stresses according to http://heronjournal.nl/53-4/3.pdf # - See post: # https://forum.freecadweb.org/viewtopic.php?f=18&t=28821 @@ -715,29 +766,29 @@ def calculate_rho(stress_tensor, fy): # i1=sxx+syy+szz NOT USED # i2=sxx*syy+syy*szz+szz*sxx-sxy**2-sxz**2-syz**2 NOT USED - i3 = (sxx * syy * szz + 2 * sxy * sxz * syz - sxx * syz**2 - - syy * sxz**2 - szz * sxy**2) + i3 = (sxx * syy * szz + 2 * sxy * sxz * syz - sxx * syz ** 2 + - syy * sxz ** 2 - szz * sxy ** 2) # Solution (5) - d = (sxx * syy - sxy**2) + d = (sxx * syy - sxy ** 2) if d != 0.: rhoz[0] = i3 / d / fy # Solution (6) - d = (sxx * szz - sxz**2) + d = (sxx * szz - sxz ** 2) if d != 0.: rhoy[1] = i3 / d / fy # Solution (7) - d = (syy * szz - syz**2) + d = (syy * szz - syz ** 2) if d != 0.: rhox[2] = i3 / d / fy # Solution (9) if sxx != 0.: fc = sxz * sxy / sxx - syz - fxy = sxy**2 / sxx - fxz = sxz**2 / sxx + fxy = sxy ** 2 / sxx + fxz = sxz ** 2 / sxx # Solution (9+) rhoy[3] = syy - fxy + fc @@ -754,8 +805,8 @@ def calculate_rho(stress_tensor, fy): # Solution (10) if syy != 0.: fc = syz * sxy / syy - sxz - fxy = sxy**2 / syy - fyz = syz**2 / syy + fxy = sxy ** 2 / syy + fyz = syz ** 2 / syy # Solution (10+) rhox[5] = sxx - fxy + fc @@ -773,8 +824,8 @@ def calculate_rho(stress_tensor, fy): # Solution (11) if szz != 0.: fc = sxz * syz / szz - sxy - fxz = sxz**2 / szz - fyz = syz**2 / szz + fxz = sxz ** 2 / szz + fyz = syz ** 2 / szz # Solution (11+) rhox[7] = sxx - fxz + fc @@ -825,10 +876,10 @@ def calculate_rho(stress_tensor, fy): scyy = syy - rhoy[ir] * fy sczz = szz - rhoz[ir] * fy ic1 = (scxx + scyy + sczz) - ic2 = (scxx * scyy + scyy * sczz + sczz * scxx - sxy**2 - - sxz**2 - syz**2) - ic3 = (scxx * scyy * sczz + 2 * sxy * sxz * syz - scxx * syz**2 - - scyy * sxz**2 - sczz * sxy**2) + ic2 = (scxx * scyy + scyy * sczz + sczz * scxx - sxy ** 2 + - sxz ** 2 - syz ** 2) + ic3 = (scxx * scyy * sczz + 2 * sxy * sxz * syz - scxx * syz ** 2 + - scyy * sxz ** 2 - sczz * sxy ** 2) if ic1 <= 1.e-6 and ic2 >= -1.e-6 and ic3 <= 1.0e-6: @@ -843,7 +894,7 @@ def calculate_rho(stress_tensor, fy): def calculate_mohr_coulomb(prin1, prin3, phi, fck): # - # HarryvL - Calculation of Mohr Coulomb yield criterion to judge + # Calculation of Mohr Coulomb yield criterion to judge # concrete curshing and shear failure # phi: angle of internal friction # fck: factored compressive strength of the matrix material (usually concrete)