[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] [getfem-commits] branch master updated: Add finite stra
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] [getfem-commits] branch master updated: Add finite strain plasticity examples with linear isotropic hardening |
|
Date: |
Fri, 25 Dec 2020 19:04:14 -0500 |
This is an automated email from the git hooks/post-receive script.
logari81 pushed a commit to branch master
in repository getfem.
The following commit(s) were added to refs/heads/master by this push:
new 4ff15ca Add finite strain plasticity examples with linear isotropic
hardening
4ff15ca is described below
commit 4ff15ca9bd431068e9318a645a320edf383a4e13
Author: Konstantinos Poulios <logari81@gmail.com>
AuthorDate: Sat Dec 26 01:03:59 2020 +0100
Add finite strain plasticity examples with linear isotropic hardening
- 3D
- plane strain (2D)
- axisymmetric (2D)
---
configure.ac | 1 +
contrib/Makefile.am | 3 +-
contrib/{ => continuum_mechanics}/Makefile.am | 20 +-
...ty_finite_strain_linear_hardening_tension_3D.py | 250 +++++++++++++++++++++
...strain_linear_hardening_tension_axisymmetric.py | 239 ++++++++++++++++++++
...strain_linear_hardening_tension_plane_strain.py | 237 +++++++++++++++++++
6 files changed, 746 insertions(+), 4 deletions(-)
diff --git a/configure.ac b/configure.ac
index a19ba7e..14fb370 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1192,6 +1192,7 @@ contrib/level_set_contact/Makefile
\
contrib/static_contact_gears/Makefile \
contrib/test_plasticity/Makefile \
contrib/opt_assembly/Makefile \
+contrib/continuum_mechanics/Makefile \
bin/Makefile \
interface/Makefile \
interface/src/Makefile \
diff --git a/contrib/Makefile.am b/contrib/Makefile.am
index 02a044e..5d1d655 100644
--- a/contrib/Makefile.am
+++ b/contrib/Makefile.am
@@ -17,4 +17,5 @@
SUBDIRS = icare delaminated_crack aposteriori xfem_stab_unilat_contact \
bimaterial_crack_test mixed_elastostatic xfem_contact crack_plate \
- static_contact_gears level_set_contact test_plasticity opt_assembly
+ static_contact_gears level_set_contact test_plasticity opt_assembly \
+ continuum_mechanics
diff --git a/contrib/Makefile.am b/contrib/continuum_mechanics/Makefile.am
similarity index 65%
copy from contrib/Makefile.am
copy to contrib/continuum_mechanics/Makefile.am
index 02a044e..fddd257 100644
--- a/contrib/Makefile.am
+++ b/contrib/continuum_mechanics/Makefile.am
@@ -15,6 +15,20 @@
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
-SUBDIRS = icare delaminated_crack aposteriori xfem_stab_unilat_contact \
- bimaterial_crack_test mixed_elastostatic xfem_contact crack_plate \
- static_contact_gears level_set_contact test_plasticity opt_assembly
+check_PROGRAMS =
+
+CLEANFILES =
+
+if BUILDPYTHON
+TESTS = plasticity_finite_strain_linear_hardening_tension_plane_strain.py
+
+AM_TESTS_ENVIRONMENT = \
+ export PYTHONPATH=$(top_builddir)/interface/src/python; \
+ export LD_LIBRARY_PATH=$(LD_LIBRARY_PATH):$(top_builddir)/src/.libs;
+LOG_COMPILER = $(PYTHON)
+endif
+
+EXTRA_DIST = \
+ plasticity_finite_strain_linear_hardening_tension_3D.py \
+ plasticity_finite_strain_linear_hardening_tension_axisymmetric.py \
+ plasticity_finite_strain_linear_hardening_tension_plane_strain.py
diff --git
a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py
b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py
new file mode 100644
index 0000000..f041500
--- /dev/null
+++
b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_3D.py
@@ -0,0 +1,250 @@
+#!/usr/bin/env python3
+# -*- coding: UTF8 -*-
+# Python GetFEM interface
+#
+# Copyright (C) 2020-2020 Konstantinos Poulios.
+#
+# This file is a part of GetFEM
+#
+# GetFEM is free software; you can redistribute it and/or modify it
+# under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version along with the GCC Runtime Library
+# Exception either version 3.1 or (at your option) any later version.
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+# License and GCC Runtime Library Exception for more details.
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+############################################################################
+""" This example simulates necking during uniaxial tension of an
+ axisymmetric rod, with a hyperelastoplastic constitutive law with
+ linear isotropic hardening. Only one eighth of the rod is actually
+ modeled, using hexahedral 3D elements.
+
+ This is a reference implementation of finite strain plasticity with
+ linear hardening in 3D using the generic weak form language of GetFEM.
+"""
+import getfem as gf
+import numpy as np
+import os, sys, subprocess
+import time
+import shutil
+gf.util_trace_level(1)
+gf.util_warning_level(1)
+
+# Input data
+L = 2*26.667 # block length
+H = 2*6.413 # block height/diameter
+dH = 0.018*H # height reduction at the center of the block
+
+N_L = 16 # number of elements in block length direction
+N_R1 = 4 # number of elements in radial direction (core)
+N_R2 = 2 # number of elements in radial direction (peel)
+
+E = 210e3 # Young's modulus
+nu = 0.3 # Poisson's ratio
+pl_sigma_0 = 5e2 # Initial yield stress
+pl_H = 21e1 # Plastic modulus (0.1% of E)
+
+disp = 8. # maximum displacement
+steps_t = 200 # number of load steps for increasing the load
+
+#------------------------------------
+geotrans2d = "GT_QK(2,2)" # geometric transformation
+
+disp_fem_order = 2 # displacements finite element order
+mult_fem_order = 2 # dirichlet multipliers finite element order
+
+#integration_degree = 3 # 4 gauss points per quad
+integration_degree = 5 # 9 gauss points per quad
+
+#------------------------------------
+resultspath = "./results"
+if not os.path.exists(resultspath):
+ os.makedirs(resultspath)
+
+tee = subprocess.Popen(["tee", "%s/tension_3D.log" % resultspath],
+ stdin=subprocess.PIPE)
+sys.stdout.flush()
+os.dup2(tee.stdin.fileno(), sys.stdout.fileno())
+sys.stderr.flush()
+os.dup2(tee.stdin.fileno(), sys.stderr.fileno())
+
+# auxiliary constants
+XM_RG = 1
+XP_RG = 2
+YM_RG = 3
+ZM_RG = 4
+
+# Mesh
+mesh2d = gf.Mesh("import", "structured_ball",
+ "GT='%s';ORG=[0,0];SIZES=[%f];NSUBDIV=[%i,%i];SYMMETRIES=2"
+ % (geotrans2d, H/2., N_R1, N_R2))
+mesh = gf.Mesh("prismatic", mesh2d, N_L, disp_fem_order)
+
+trsf = np.zeros([3,3])
+trsf[0,2] = L/2.
+trsf[2,1] = 1
+trsf[1,0] = 1
+mesh.transform(trsf)
+
+N = mesh.dim()
+
+outer_faces = mesh.outer_faces()
+outer_normals = mesh.normal_of_faces(outer_faces)
+mesh.set_region(XM_RG,
+ outer_faces[:,np.nonzero(outer_normals[0] < -0.95)[0]])
+mesh.set_region(XP_RG,
+ outer_faces[:,np.nonzero(outer_normals[0] > 0.95)[0]])
+mesh.set_region(YM_RG,
+ outer_faces[:,np.nonzero(outer_normals[1] < -0.95)[0]])
+mesh.set_region(ZM_RG,
+ outer_faces[:,np.nonzero(outer_normals[2] < -0.95)[0]])
+
+if dH > 0:
+ pts = mesh.pts()
+ dr = 0.2 # fine element size ratio w.r.t. uniform mesh size
+ r0 = 0.04 # ratio of the finer meshed area w.r.t. the total length
+ x0 = r0/dr
+ b2 = (1-dr)/(x0-1)**2
+ b1 = dr - 2*b2*x0
+ b0 = 1 - b1 -b2
+ for i in range(pts.shape[1]):
+ x = pts[0,i]
+ y = pts[1,i]
+ z = pts[2,i]
+ r = np.sqrt(y**2+z**2)
+ t = 2*abs(x)/L
+ if t < x0:
+ x *= dr;
+ else:
+ x = (b0 + b1*t + b2*t**2) * np.sign(x)*L/2
+ pts[0,i] = x
+ pts[1,i] -= (y*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L))
+ pts[2,i] -= (z*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L))
+ mesh.set_pts(pts)
+
+mesh.export_to_vtk("%s/mesh.vtk" % resultspath)
+
+# FEM
+mfN = gf.MeshFem(mesh, N)
+mfN.set_classical_fem(disp_fem_order)
+#if disp_fem_order == 2:
+# mfN.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(3)"))
+#else:
+# mfN.set_classical_fem(disp_fem_order)
+keptdofs = np.arange(mfN.nbdof())
+keptdofs = np.setdiff1d(keptdofs, mfN.basic_dof_on_region(XM_RG)[0::N])
+keptdofs = np.setdiff1d(keptdofs, mfN.basic_dof_on_region(YM_RG)[1::N])
+keptdofs = np.setdiff1d(keptdofs, mfN.basic_dof_on_region(ZM_RG)[2::N])
+mfu = gf.MeshFem("partial", mfN, keptdofs)
+
+mfmult = gf.MeshFem(mesh, 1)
+mfmult.set_classical_fem(mult_fem_order)
+
+mfout1 = gf.MeshFem(mesh)
+mfout1.set_classical_discontinuous_fem(disp_fem_order-1)
+mfout2 = gf.MeshFem(mesh)
+mfout2.set_classical_discontinuous_fem(disp_fem_order)
+
+mim = gf.MeshIm(mesh, integration_degree)
+
+mimd1 = gf.MeshImData(mim)
+mimd6 = gf.MeshImData(mim, -1, 6)
+
+# Model
+md = gf.Model("real")
+
+md.add_fem_variable("u", mfu)
+
+# Vertical displacement
+md.add_initialized_data("disp", [0.])
+
+md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
+md.add_initialized_data("mu", E/(2*(1+nu))) # Shear modulus
+md.add_initialized_data("Y", np.sqrt(2./3.)*pl_sigma_0) # Initial yield limit
+md.add_macro("F", "Id(3)+Grad_u")
+md.add_macro("J", "Det(F)")
+md.add_macro("tauH", "K*log(J)")
+md.add_im_data("gamma0", mimd1) # accumulated plastic
strain at previous time step
+md.add_im_data("invCp0vec", mimd6) # Components 11, 22,
33, 12, 23 and 31 of the plastic
+md.set_variable("invCp0vec", # part ofthe inverse
right Cauchy Green tensor at the
+ np.tile([1,1,1,0,0,0], mimd6.nbpts())) # previous step.
Symmetric components are omitted.
+md.add_macro("invCp0", "[[[1,0,0],[0,0,0],[0,0,0]],"+\
+ " [[0,0,0],[0,1,0],[0,0,0]],"+\
+ " [[0,0,0],[0,0,0],[0,0,1]],"+\
+ " [[0,1,0],[1,0,0],[0,0,0]],"+\
+ " [[0,0,0],[0,0,1],[0,1,0]],"+\
+ " [[0,0,1],[0,0,0],[1,0,0]]].invCp0vec") #Vec6ToMat3x3
+md.add_macro("devlogbetr", "Deviator(Logm(F*invCp0*F'))")
+md.add_macro("Y0","{A}+{B}*gamma0".format(A=np.sqrt(2./3.)*pl_sigma_0,
B=2./3.*pl_H))
+md.add_macro("ksi",
+
"(1-1/max(1,mu*pow(J,-5/3)*Norm(devlogbetr)/Y0))/(2+{B}/(mu*pow(J,-5/3)))"
+ .format(B=2./3.*pl_H))
+md.add_macro("gamma", "gamma0+ksi*Norm(devlogbetr)")
+md.add_macro("devlogbe", "(1-2*ksi)*devlogbetr")
+md.add_macro("tauD", "mu*pow(J,-2/3)*devlogbe")
+
+md.add_nonlinear_generic_assembly_brick\
+ (mim, "((tauH*Id(3)+tauD)*Inv(F')):Grad_Test_u")
+
+md.add_macro("sigmaD", "(mu*pow(J,-5/3)*devlogbe)")
+md.add_macro("sigma", "tauH/J*Id(3)+mu*pow(J,-5/3)*devlogbe")
+md.add_macro("invCp", "(Inv(F)*Expm(-ksi*devlogbetr)*(F))*invCp0"\
+ "*(Inv(F)*Expm(-ksi*devlogbetr)*(F))'")
+
+# Dirichlet condition
+md.add_filtered_fem_variable("dirmult", mfmult, XP_RG)
+md.add_nonlinear_generic_assembly_brick(mim, "(disp-u(1))*dirmult", XP_RG)
+
+print("Model dofs: %i" % md.nbdof())
+print("Displacement fem dofs: %i" % mfu.nbdof())
+print("Dirichlet mult dofs: %i" % md.mesh_fem_of_variable("dirmult").nbdof())
+
+shutil.copyfile(os.path.abspath(sys.argv[0]),resultspath+"/"+sys.argv[0])
+starttime_overall = time.process_time()
+with open("%s/tension_3D.dat" % resultspath, "w") as f1:
+ for step in range(steps_t+1):
+ md.set_variable("disp", disp*step/float(steps_t))
+ print('STEP %i: Solving with disp = %g' % (step, md.variable("disp")))
+
+ starttime = time.process_time()
+ md.solve('noisy', 'max_iter', 25, 'max_res', 1e-10,
+ 'lsearch', 'simplest', 'alpha max ratio', 100, 'alpha min',
0.2, 'alpha mult', 0.6,
+ 'alpha threshold res', 1e9)
+ print('STEP %i COMPLETED IN %f SEC' % (step,
time.process_time()-starttime))
+
+ F = gf.asm_generic(mim, 0, "dirmult", XP_RG, md)
+ print("Displacement %g, total force %g" % (md.variable("disp"), F))
+
+ A = gf.asm_generic(mim, 0, "Norm(J*Inv(F')*[1;0;0])", XP_RG, md)
+ V = gf.asm_generic(mim, 0, "1", -1, md)
+ sigma11 = gf.asm_generic(mim, 0, "sigma(1,1)", -1, md)/V
+ gamma = gf.asm_generic(mim, 0, "gamma", -1, md)/V
+ f1.write("%.10g %.10g %.10g %.10g %10g %10g\n"
+ % (md.variable("disp"), F, A, F/A, sigma11, gamma))
+ f1.flush()
+
+ output = (mfout1, md.local_projection(mim, "sqrt(1.5)*Norm(sigmaD)",
mfout1), "Von Mises Stress",
+ mfout1, md.local_projection(mim, "J", mfout1), "J",
+ mfout1, md.local_projection(mim, "sigma(1,1)", mfout1),
"Cauchy stress 11",
+ mfout1, md.local_projection(mim, "sigma(2,2)", mfout1),
"Cauchy stress 22",
+ mfout1, md.local_projection(mim, "sigma(1,2)", mfout1),
"Cauchy stress 12",
+ mfout1, md.local_projection(mim, "sigma(3,3)", mfout1),
"Cauchy stress 33",
+ mfu, md.variable("u"), "Displacements",
+ mfout2, md.interpolation("dirmult", mfout2, XP_RG), "Nominal
reaction traction",
+ mfout2, md.local_projection(mim, "gamma", mfout2), "plastic
strain")
+ mfout2.export_to_vtk("%s/tension_3D_%i.vtk" % (resultspath, step),
*output)
+
+ md.set_variable("gamma0", md.interpolation("gamma", mimd1, -1))
+ md.set_variable("invCp0vec",
+ md.interpolation("[[[1,0,0,0,0,0]
,[0,0,0,0.5,0,0],[0,0,0,0,0,0.5]],"+\
+ " [[0,0,0,0.5,0,0],[0,1,0,0,0,0]
,[0,0,0,0,0.5,0]],"+\
+ "
[[0,0,0,0,0,0.5],[0,0,0,0,0.5,0],[0,0,1,0,0,0]]]:invCp", mimd6, -1))
+
+print('OVERALL SOLUTION TIME IS %f SEC' %
(time.process_time()-starttime_overall))
+
diff --git
a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_axisymmetric.py
b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_axisymmetric.py
new file mode 100644
index 0000000..8fc92ae
--- /dev/null
+++
b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_axisymmetric.py
@@ -0,0 +1,239 @@
+#!/usr/bin/env python3
+# -*- coding: UTF8 -*-
+# Python GetFEM interface
+#
+# Copyright (C) 2020-2020 Konstantinos Poulios.
+#
+# This file is a part of GetFEM
+#
+# GetFEM is free software; you can redistribute it and/or modify it
+# under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version along with the GCC Runtime Library
+# Exception either version 3.1 or (at your option) any later version.
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+# License and GCC Runtime Library Exception for more details.
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+############################################################################
+""" This example simulates necking during uniaxial tension of an
+ axisymmetric rod, with a hyperelastoplastic constitutive law with
+ linear isotropic hardening. Only one quarter of the longitudinal
+ section of the rod is actually modeled, using quadrilateral 2D
+ elements in an axisymmetric formulation.
+
+ This is a reference implementation of axisymmetric (i.e. 2D) finite
+ strain plasticity with linear hardening using the generic weak form
+ language of GetFEM.
+"""
+import getfem as gf
+import numpy as np
+import os, sys, subprocess
+import time
+import shutil
+gf.util_trace_level(1)
+gf.util_warning_level(1)
+
+# Input data
+L = 2*26.667 # block length
+H = 2*6.413 # block height
+dH = 0.018*H # height reduction at the center of the block
+
+N_L = 20 # number of elements in block length direction
+N_H = 10 # number of elements in block height direction
+
+E = 210e3 # Young's modulus
+nu = 0.3 # Poisson's ratio
+pl_sigma_0 = 5e2 # Initial yield stress
+pl_H = 21e1 # Plastic modulus (0.1% of E)
+
+disp = 8. # maximum displacement
+steps_t = 200 # number of load steps for increasing the load
+
+#------------------------------------
+geotrans = "GT_QK(2,2)" # geometric transformation
+
+disp_fem_order = 2 # displacements finite element order
+mult_fem_order = 2 # dirichlet multipliers finite element order
+
+#integration_degree = 3 # 4 gauss points per quad
+integration_degree = 5 # 9 gauss points per quad
+
+#------------------------------------
+resultspath = "./results"
+if not os.path.exists(resultspath):
+ os.makedirs(resultspath)
+
+tee = subprocess.Popen(["tee", "%s/tension_axisymmetric.log" % resultspath],
+ stdin=subprocess.PIPE)
+sys.stdout.flush()
+os.dup2(tee.stdin.fileno(), sys.stdout.fileno())
+sys.stderr.flush()
+os.dup2(tee.stdin.fileno(), sys.stderr.fileno())
+
+# auxiliary constants
+XM_RG = 1
+XP_RG = 2
+YM_RG = 3
+
+# Mesh
+xmin = 0.
+ymin = 0.
+dx = L/2
+dy = H/2
+mesh = gf.Mesh("import", "structured",
+ "GT='%s';ORG=[%f,%f];SIZES=[%f,%f];NSUBDIV=[%i,%i]"
+ % (geotrans, xmin, ymin, dx, dy, N_L, N_H))
+N = mesh.dim()
+
+outer_faces = mesh.outer_faces()
+outer_normals = mesh.normal_of_faces(outer_faces)
+mesh.set_region(XM_RG,
+ outer_faces[:,np.nonzero(outer_normals[0] < -0.95)[0]])
+mesh.set_region(XP_RG,
+ outer_faces[:,np.nonzero(outer_normals[0] > 0.95)[0]])
+mesh.set_region(YM_RG,
+ outer_faces[:,np.nonzero(outer_normals[1] < -0.95)[0]])
+
+if dH > 0:
+ pts = mesh.pts()
+ dy = 0.2 # fine element size ratio w.r.t. uniform mesh size
+ y0 = 0.04 # ratio of the finer meshed area w.r.t. the total length
+ x0 = y0/dy
+ b2 = (1-dy)/(x0-1)**2
+ b1 = dy - 2*b2*x0
+ b0 = 1 - b1 -b2
+ for i in range(pts.shape[1]):
+ x = pts[0,i]
+ y = pts[1,i]
+ t = 2*abs(x)/L
+ if t < x0:
+ x *= dy;
+ else:
+ x = (b0 + b1*t + b2*t**2) * np.sign(x)*L/2
+ pts[0,i] = x
+ pts[1,i] -= (y*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L))
+ mesh.set_pts(pts)
+
+mesh.export_to_vtk("%s/mesh.vtk" % resultspath)
+
+# FEM
+mfN = gf.MeshFem(mesh, N)
+mfN.set_classical_fem(disp_fem_order)
+#if disp_fem_order == 2:
+# mfN.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(2)"))
+#else:
+# mfN.set_classical_fem(disp_fem_order)
+keptdofs = np.arange(mfN.nbdof())
+keptdofs = np.setdiff1d(keptdofs, mfN.basic_dof_on_region(XM_RG)[0::N])
+keptdofs = np.setdiff1d(keptdofs, mfN.basic_dof_on_region(YM_RG)[1::N])
+mfu = gf.MeshFem("partial", mfN, keptdofs)
+
+mfmult = gf.MeshFem(mesh, 1)
+mfmult.set_classical_fem(mult_fem_order)
+
+mfout1 = gf.MeshFem(mesh)
+mfout1.set_classical_discontinuous_fem(disp_fem_order-1)
+mfout2 = gf.MeshFem(mesh)
+mfout2.set_classical_discontinuous_fem(disp_fem_order)
+
+mim = gf.MeshIm(mesh, integration_degree)
+
+mimd1 = gf.MeshImData(mim)
+mimd4 = gf.MeshImData(mim, -1, 4)
+
+# Model
+md = gf.Model("real")
+
+md.add_fem_variable("u", mfu)
+
+# Vertical displacement
+md.add_initialized_data("disp", [0.])
+
+md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
+md.add_initialized_data("mu", E/(2*(1+nu))) # Shear modulus
+md.add_initialized_data("Y", np.sqrt(2./3.)*pl_sigma_0) # Initial yield limit
+md.add_macro("F", "Id(2)+Grad_u")
+md.add_macro("F3d",
"Id(3)+[0,0,0;0,0,0;0,0,1/X(2)]*u(2)+[1,0;0,1;0,0]*Grad_u*[1,0,0;0,1,0]")
+md.add_macro("J", "Det(F)*(1+u(2)/X(2))")
+md.add_macro("tauH", "K*log(J)")
+md.add_im_data("gamma0", mimd1) # accumulated plastic
strain at previous time step
+md.add_im_data("invCp0vec", mimd4) # Components 11, 22, 33
and 12 of the plastic part of
+md.set_variable("invCp0vec", # the inverse right
Cauchy Green tensor at the previous
+ np.tile([1,1,1,0], mimd4.nbpts())) # step. Symmetric
components are omitted.
+md.add_macro("invCp0", "[[[1,0,0],[0,0,0],[0,0,0]],"+\
+ " [[0,0,0],[0,1,0],[0,0,0]],"+\
+ " [[0,0,0],[0,0,0],[0,0,1]],"+\
+ " [[0,1,0],[1,0,0],[0,0,0]]].invCp0vec") #Vec4ToMat3x3
+md.add_macro("devlogbetr", "Deviator(Logm(F3d*invCp0*F3d'))")
+md.add_macro("Y0","{A}+{B}*gamma0".format(A=np.sqrt(2./3.)*pl_sigma_0,
B=2./3.*pl_H))
+md.add_macro("ksi",
+
"(1-1/max(1,mu*pow(J,-5/3)*Norm(devlogbetr)/Y0))/(2+{B}/(mu*pow(J,-5/3)))"
+ .format(B=2./3.*pl_H))
+md.add_macro("gamma", "gamma0+ksi*Norm(devlogbetr)")
+md.add_macro("devlogbe", "(1-2*ksi)*devlogbetr")
+md.add_macro("tauD2d", "mu*pow(J,-2/3)*[1,0,0;0,1,0]*devlogbe*[1,0;0,1;0,0]")
+md.add_macro("tauD33", "mu*pow(J,-2/3)*devlogbe(3,3)")
+
+md.add_nonlinear_generic_assembly_brick\
+ (mim,
"2*pi*X(2)*(((tauH*Id(2)+tauD2d)*Inv(F')):Grad_Test_u+(tauH+tauD33)/(X(2)+u(2))*Test_u(2))")
+
+md.add_macro("sigmaD", "(mu*pow(J,-5/3)*devlogbe)")
+md.add_macro("sigma", "tauH/J*Id(3)+mu*pow(J,-5/3)*devlogbe")
+md.add_macro("invCp", "(Inv(F3d)*Expm(-ksi*devlogbetr)*(F3d))*invCp0"\
+ "*(Inv(F3d)*Expm(-ksi*devlogbetr)*(F3d))'")
+
+# Dirichlet condition
+md.add_filtered_fem_variable("dirmult", mfmult, XP_RG)
+md.add_nonlinear_generic_assembly_brick(mim, "2*pi*X(2)*(disp-u(1))*dirmult",
XP_RG)
+
+print("Model dofs: %i" % md.nbdof())
+print("Displacement fem dofs: %i" % mfu.nbdof())
+print("Dirichlet mult dofs: %i" % md.mesh_fem_of_variable("dirmult").nbdof())
+
+shutil.copyfile(os.path.abspath(sys.argv[0]),resultspath+"/"+sys.argv[0])
+starttime_overall = time.process_time()
+with open("%s/tension_axisymmetric.dat" % resultspath, "w") as f1:
+ for step in range(steps_t+1):
+ md.set_variable("disp", disp*step/float(steps_t))
+ print('STEP %i: Solving with disp = %g' % (step, md.variable("disp")))
+
+ starttime = time.process_time()
+ md.solve('noisy', 'max_iter', 25, 'max_res', 1e-10,
+ 'lsearch', 'simplest', 'alpha max ratio', 100, 'alpha min',
0.2, 'alpha mult', 0.6,
+ 'alpha threshold res', 1e9)
+ print('STEP %i COMPLETED IN %f SEC' % (step,
time.process_time()-starttime))
+
+ F = gf.asm_generic(mim, 0, "2*pi*X(2)*dirmult", XP_RG, md)
+ print("Displacement %g, total force %g" % (md.variable("disp"), F))
+ A = gf.asm_generic(mim, 0, "2*pi*X(2)*Norm(J*Inv(F')*[1;0])", XP_RG, md)
+ V = gf.asm_generic(mim, 0, "2*pi*X(2)", -1, md)
+ sigma11 = gf.asm_generic(mim, 0, "sigma(1,1)", -1, md)/V
+ gamma = gf.asm_generic(mim, 0, "gamma", -1, md)/V
+ f1.write("%.10g %.10g %.10g %.10g %10g %10g\n"
+ % (md.variable("disp"), F, A, F/A, sigma11, gamma))
+ f1.flush()
+
+ output = (mfout1, md.local_projection(mim, "sqrt(1.5)*Norm(sigmaD)",
mfout1), "Von Mises Stress",
+ mfout1, md.local_projection(mim, "J", mfout1), "J",
+ mfout1, md.local_projection(mim, "sigma(1,1)", mfout1),
"Cauchy stress 11",
+ mfout1, md.local_projection(mim, "sigma(2,2)", mfout1),
"Cauchy stress 22",
+ mfout1, md.local_projection(mim, "sigma(1,2)", mfout1),
"Cauchy stress 12",
+ mfout1, md.local_projection(mim, "sigma(3,3)", mfout1),
"Cauchy stress 33",
+ mfu, md.variable("u"), "Displacements",
+ mfout2, md.interpolation("dirmult", mfout2, XP_RG), "Nominal
reaction traction",
+ mfout2, md.local_projection(mim, "gamma", mfout2), "plastic
strain")
+ mfout2.export_to_vtk("%s/tension_axisymmetric_%i.vtk" % (resultspath,
step), *output)
+
+ md.set_variable("gamma0", md.interpolation("gamma", mimd1, -1))
+ md.set_variable("invCp0vec",
+ md.interpolation("[[[1,0,0,0]
,[0,0,0,0.5],[0,0,0,0]],"+\
+ " [[0,0,0,0.5],[0,1,0,0]
,[0,0,0,0]],"+\
+ " [[0,0,0,0] ,[0,0,0,0]
,[0,0,1,0]]]:invCp", mimd4, -1))
+
+print('OVERALL SOLUTION TIME IS %f SEC' %
(time.process_time()-starttime_overall))
+
diff --git
a/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_plane_strain.py
b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_plane_strain.py
new file mode 100644
index 0000000..422efe1
--- /dev/null
+++
b/contrib/continuum_mechanics/plasticity_finite_strain_linear_hardening_tension_plane_strain.py
@@ -0,0 +1,237 @@
+#!/usr/bin/env python3
+# -*- coding: UTF8 -*-
+# Python GetFEM interface
+#
+# Copyright (C) 2020-2020 Konstantinos Poulios.
+#
+# This file is a part of GetFEM
+#
+# GetFEM is free software; you can redistribute it and/or modify it
+# under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version along with the GCC Runtime Library
+# Exception either version 3.1 or (at your option) any later version.
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+# License and GCC Runtime Library Exception for more details.
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+############################################################################
+""" This example simulates necking during uniaxial tension of a plate
+ under plane strain, with a hyperelastoplastic constitutive law with
+ linear isotropic hardening. Only the cross section of the plate is
+ modeled, using quadrilateral 2D elements.
+
+ This is a reference implementation of finite strain plasticity with
+ linear hardening in 2D (plane strain) using the generic weak form
+ language of GetFEM.
+"""
+import getfem as gf
+import numpy as np
+import os, sys, subprocess
+import time
+import shutil
+gf.util_trace_level(1)
+gf.util_warning_level(1)
+
+# Input data
+L = 2*26.667 # block length
+H = 2*6.413 # block height
+dH = 0.018*H # height reduction at the center of the block
+
+N_L = 20 # number of elements in block length direction
+N_H = 10 # number of elements in block height direction
+
+E = 210e3 # Young's modulus
+nu = 0.3 # Poisson's ratio
+pl_sigma_0 = 5e2 # Initial yield stress
+pl_H = 21e1 # Plastic modulus (0.1% of E)
+
+disp = 8. # maximum displacement
+steps_t = 200 # number of load steps for increasing the load
+
+#------------------------------------
+geotrans = "GT_QK(2,2)" # geometric transformation
+
+disp_fem_order = 2 # displacements finite element order
+mult_fem_order = 2 # dirichlet multipliers finite element order
+
+#integration_degree = 3 # 4 gauss points per quad
+integration_degree = 5 # 9 gauss points per quad
+
+#------------------------------------
+resultspath = "./results"
+if not os.path.exists(resultspath):
+ os.makedirs(resultspath)
+
+tee = subprocess.Popen(["tee", "%s/tension_plane_strain.log" % resultspath],
+ stdin=subprocess.PIPE)
+sys.stdout.flush()
+os.dup2(tee.stdin.fileno(), sys.stdout.fileno())
+sys.stderr.flush()
+os.dup2(tee.stdin.fileno(), sys.stderr.fileno())
+
+# auxiliary constants
+XM_RG = 1
+XP_RG = 2
+YM_RG = 3
+
+# Mesh
+xmin = 0.
+ymin = 0.
+dx = L/2
+dy = H/2
+mesh = gf.Mesh("import", "structured",
+ "GT='%s';ORG=[%f,%f];SIZES=[%f,%f];NSUBDIV=[%i,%i]"
+ % (geotrans, xmin, ymin, dx, dy, N_L, N_H))
+N = mesh.dim()
+
+outer_faces = mesh.outer_faces()
+outer_normals = mesh.normal_of_faces(outer_faces)
+mesh.set_region(XM_RG,
+ outer_faces[:,np.nonzero(outer_normals[0] < -0.95)[0]])
+mesh.set_region(XP_RG,
+ outer_faces[:,np.nonzero(outer_normals[0] > 0.95)[0]])
+mesh.set_region(YM_RG,
+ outer_faces[:,np.nonzero(outer_normals[1] < -0.95)[0]])
+
+if dH > 0:
+ pts = mesh.pts()
+ dy = 0.2 # fine element size ratio w.r.t. uniform mesh size
+ y0 = 0.04 # ratio of the finer meshed area w.r.t. the total length
+ x0 = y0/dy
+ b2 = (1-dy)/(x0-1)**2
+ b1 = dy - 2*b2*x0
+ b0 = 1 - b1 -b2
+ for i in range(pts.shape[1]):
+ x = pts[0,i]
+ y = pts[1,i]
+ t = 2*abs(x)/L
+ if t < x0:
+ x *= dy;
+ else:
+ x = (b0 + b1*t + b2*t**2) * np.sign(x)*L/2
+ pts[0,i] = x
+ pts[1,i] -= (y*dH)/(2*H) * (1 + np.cos(2.*np.pi*x/L))
+ mesh.set_pts(pts)
+
+mesh.export_to_vtk("%s/mesh.vtk" % resultspath)
+
+# FEM
+mfN = gf.MeshFem(mesh, N)
+mfN.set_classical_fem(disp_fem_order)
+#if disp_fem_order == 2:
+# mfN.set_fem(gf.Fem("FEM_Q2_INCOMPLETE(2)"))
+#else:
+# mfN.set_classical_fem(disp_fem_order)
+keptdofs = np.arange(mfN.nbdof())
+keptdofs = np.setdiff1d(keptdofs, mfN.basic_dof_on_region(XM_RG)[0::N])
+keptdofs = np.setdiff1d(keptdofs, mfN.basic_dof_on_region(YM_RG)[1::N])
+mfu = gf.MeshFem("partial", mfN, keptdofs)
+
+mfmult = gf.MeshFem(mesh, 1)
+mfmult.set_classical_fem(mult_fem_order)
+
+mfout1 = gf.MeshFem(mesh)
+mfout1.set_classical_discontinuous_fem(disp_fem_order-1)
+mfout2 = gf.MeshFem(mesh)
+mfout2.set_classical_discontinuous_fem(disp_fem_order)
+
+mim = gf.MeshIm(mesh, integration_degree)
+
+mimd1 = gf.MeshImData(mim)
+mimd4 = gf.MeshImData(mim, -1, 4)
+
+# Model
+md = gf.Model("real")
+
+md.add_fem_variable("u", mfu)
+
+# Vertical displacement
+md.add_initialized_data("disp", [0.])
+
+md.add_initialized_data("K", E/(3.*(1.-2.*nu))) # Bulk modulus
+md.add_initialized_data("mu", E/(2*(1+nu))) # Shear modulus
+md.add_initialized_data("Y", np.sqrt(2./3.)*pl_sigma_0) # Initial yield limit
+md.add_macro("F", "Id(2)+Grad_u")
+md.add_macro("F3d", "Id(3)+[1,0;0,1;0,0]*Grad_u*[1,0,0;0,1,0]")
+md.add_macro("J", "Det(F)")
+md.add_macro("tauH", "K*log(J)")
+md.add_im_data("gamma0", mimd1) # accumulated plastic
strain at previous time step
+md.add_im_data("invCp0vec", mimd4) # Components 11, 22, 33
and 12 of the plastic part of
+md.set_variable("invCp0vec", # the inverse right
Cauchy Green tensor at the previous
+ np.tile([1,1,1,0], mimd4.nbpts())) # step. Symmetric
components are omitted.
+md.add_macro("invCp0", "[[[1,0,0],[0,0,0],[0,0,0]],"+\
+ " [[0,0,0],[0,1,0],[0,0,0]],"+\
+ " [[0,0,0],[0,0,0],[0,0,1]],"+\
+ " [[0,1,0],[1,0,0],[0,0,0]]].invCp0vec") #Vec4ToMat3x3
+md.add_macro("devlogbetr", "Deviator(Logm(F3d*invCp0*F3d'))")
+md.add_macro("Y0","{A}+{B}*gamma0".format(A=np.sqrt(2./3.)*pl_sigma_0,
B=2./3.*pl_H))
+md.add_macro("ksi",
+
"(1-1/max(1,mu*pow(J,-5/3)*Norm(devlogbetr)/Y0))/(2+{B}/(mu*pow(J,-5/3)))"
+ .format(B=2./3.*pl_H))
+md.add_macro("gamma", "gamma0+ksi*Norm(devlogbetr)")
+md.add_macro("devlogbe", "(1-2*ksi)*devlogbetr")
+md.add_macro("tauD2d", "mu*pow(J,-2/3)*[1,0,0;0,1,0]*devlogbe*[1,0;0,1;0,0]")
+
+md.add_nonlinear_generic_assembly_brick\
+ (mim, "((tauH*Id(2)+tauD2d)*Inv(F')):Grad_Test_u")
+
+md.add_macro("sigmaD", "(mu*pow(J,-5/3)*devlogbe)")
+md.add_macro("sigma", "tauH/J*Id(3)+mu*pow(J,-5/3)*devlogbe")
+md.add_macro("invCp", "(Inv(F3d)*Expm(-ksi*devlogbetr)*(F3d))*invCp0"\
+ "*(Inv(F3d)*Expm(-ksi*devlogbetr)*(F3d))'")
+
+# Dirichlet condition
+md.add_filtered_fem_variable("dirmult", mfmult, XP_RG)
+md.add_nonlinear_generic_assembly_brick(mim, "(disp-u(1))*dirmult", XP_RG)
+
+print("Model dofs: %i" % md.nbdof())
+print("Displacement fem dofs: %i" % mfu.nbdof())
+print("Dirichlet mult dofs: %i" % md.mesh_fem_of_variable("dirmult").nbdof())
+
+shutil.copyfile(os.path.abspath(sys.argv[0]),resultspath+"/"+sys.argv[0])
+starttime_overall = time.process_time()
+with open("%s/tension_plane_strain.dat" % resultspath, "w") as f1:
+ for step in range(steps_t+1):
+ md.set_variable("disp", disp*step/float(steps_t))
+ print('STEP %i: Solving with disp = %g' % (step, md.variable("disp")))
+
+ starttime = time.process_time()
+ md.solve('noisy', 'max_iter', 25, 'max_res', 1e-10,
+ 'lsearch', 'simplest', 'alpha max ratio', 100, 'alpha min', 1,
'alpha mult', 0.6,
+ 'alpha threshold res', 1e9)
+ print('STEP %i COMPLETED IN %f SEC' % (step,
time.process_time()-starttime))
+
+ F = gf.asm_generic(mim, 0, "dirmult", XP_RG, md)
+ print("Displacement %g, total force %g" % (md.variable("disp"), F))
+ A = gf.asm_generic(mim, 0, "Norm(J*Inv(F')*[1;0])", XP_RG, md)
+ V = gf.asm_generic(mim, 0, "1", -1, md)
+ sigma11 = gf.asm_generic(mim, 0, "sigma(1,1)", -1, md)/V
+ gamma = gf.asm_generic(mim, 0, "gamma", -1, md)/V
+ f1.write("%.10g %.10g %.10g %.10g %10g %10g\n"
+ % (md.variable("disp"), F, A, F/A, sigma11, gamma))
+ f1.flush()
+
+ output = (mfout1, md.local_projection(mim, "sqrt(1.5)*Norm(sigmaD)",
mfout1), "Von Mises Stress",
+ mfout1, md.local_projection(mim, "J", mfout1), "J",
+ mfout1, md.local_projection(mim, "sigma(1,1)", mfout1),
"Cauchy stress 11",
+ mfout1, md.local_projection(mim, "sigma(2,2)", mfout1),
"Cauchy stress 22",
+ mfout1, md.local_projection(mim, "sigma(1,2)", mfout1),
"Cauchy stress 12",
+ mfout1, md.local_projection(mim, "sigma(3,3)", mfout1),
"Cauchy stress 33",
+ mfu, md.variable("u"), "Displacements",
+ mfout2, md.interpolation("dirmult", mfout2, XP_RG), "Nominal
reaction traction",
+ mfout2, md.local_projection(mim, "gamma", mfout2), "plastic
strain")
+ mfout2.export_to_vtk("%s/tension_plane_strain_%i.vtk" % (resultspath,
step), *output)
+
+ md.set_variable("gamma0", md.interpolation("gamma", mimd1, -1))
+ md.set_variable("invCp0vec",
+ md.interpolation("[[[1,0,0,0]
,[0,0,0,0.5],[0,0,0,0]],"+\
+ " [[0,0,0,0.5],[0,1,0,0]
,[0,0,0,0]],"+\
+ " [[0,0,0,0] ,[0,0,0,0]
,[0,0,1,0]]]:invCp", mimd4, -1))
+
+print('OVERALL SOLUTION TIME IS %f SEC' %
(time.process_time()-starttime_overall))
+
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] [getfem-commits] branch master updated: Add finite strain plasticity examples with linear isotropic hardening,
Konstantinos Poulios <=