kindred/create
1
1
Fork 0
Code Issues 53 Pull Requests Actions Packages Projects 9 Releases 2 Wiki Activity

[FEM] added critical strain ratio functionality (#7467)

Browse Source 
* added critical strain ratio functionality
This commit is contained in:
HarryvL
2022-09-14 01:30:17 +02:00
committed by GitHub
parent db7d615cb0
commit 9dc9d6735a
3 changed files with 90 additions and 30 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/Mod/Fem/App/FemVTKTools.cpp
View File

@@ -718,6 +718,8 @@ std::map<std::string, std::string> _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

7
src/Mod/Fem/femobjects/result_mechanical.py
View File

@@ -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

111
src/Mod/Fem/femresult/resulttools.py
View File

@@ -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)