[Getfem-commits] [getfem-commits] branch master updated: Add quasi-linear viscoelasticity example (finite strain)

[Konstantinos Poulios]

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 2b5bbe3  Add quasi-linear viscoelasticity example (finite strain)
2b5bbe3 is described below

commit 2b5bbe33723168066f400170ca795b4c2e9894db
Author: Konstantinos Poulios <logari81@gmail.com>
AuthorDate: Sun Aug 15 17:51:49 2021 +0200

    Add quasi-linear viscoelasticity example (finite strain)
    
    (de Pascalis, Abrahams, Parnell 2014)
---
 contrib/continuum_mechanics/QLV_viscoelasticity.py | 196 +++++++++++++++++++++
 1 file changed, 196 insertions(+)

diff --git a/contrib/continuum_mechanics/QLV_viscoelasticity.py 
b/contrib/continuum_mechanics/QLV_viscoelasticity.py
new file mode 100644
index 0000000..6c22f0d
--- /dev/null
+++ b/contrib/continuum_mechanics/QLV_viscoelasticity.py
@@ -0,0 +1,196 @@
+#!/usr/bin/env python3
+# -*- coding: UTF8 -*-
+# Python GetFEM interface
+#
+# Copyright (C) 2021-2021 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 provides an implementation of large strain quasilinear
+viscoelasticity according to (de Pascalis, Abrahams, Parnell 2014). It
+implements the compressible case with the Horgan-Murphy law from the paper,
+as well as a neo-Hookean law from (Simo, Taylor, Pister 1985). The
+implementation only covers Prony stress relaxation functions, which allow
+a much simpler time integration, requiring storage of data at only one
+previous time step.
+"""
+
+import getfem as gf
+
+gf.util_trace_level(1)
+gf.util_warning_level(1)
+
+# Material parameters under high strain rate
+E =   1e2                    # Young's modulus [Pa]
+nu =  0.49                   # Poisson's ratio [-]
+kappa = E/(3*(1-2*nu))
+mu = E/(2*(1+nu))
+MAT_LAW = ('neohookean-Simo','Horgan-Murphy')[1]
+gamma = 1./6.  # only used for Horgan-Murphy
+
+# Viscoelastic material constants
+tauH = 1e-3
+tauD = 1e-3
+rH = 0.5 # the low strain rate bulk modulus is rH*kappa
+rD = 0.5 # the low strain rate shear modulus is rD*mu
+
+# Dimensions [mm]
+W = 0.5     # Block width
+H = 2       # Block height
+T = 0.5     # Block thikness
+
+# Time
+t_max = 28e-3
+
+# Numerical parameters
+N_W = 1
+N_H = 1
+N_T = 1
+steps = 112
+fem_order = 2
+
+# Mesh
+m = gf.Mesh("import",
+            "structured","GT='GT_QK(3,2)'; ORG=[0,0,0]; SIZES=[%f,%f,%f]; 
NSUBDIV=[%i,%i,%i]"
+            % (W, H, T, N_W, N_H, N_T))
+
+# Boundaries
+RG_LEFT = 11
+RG_BOTTOM = 12
+RG_BACK = 13
+RG_TOP = 14
+m.set_region(RG_LEFT, m.outer_faces_in_box([-1e-5,-1e-5,-1e-5], 
[1e-5,H+1e-5,T+1e-5]))
+m.set_region(RG_BOTTOM, m.outer_faces_in_box([-1e-5,-1e-5,-1e-5], 
[W+1e-5,1e-5,T+1e-5]))
+m.set_region(RG_BACK, m.outer_faces_in_box([-1e-5,-1e-5,-1e-5], 
[W+1e-5,H+1e-5,1e-5]))
+m.set_region(RG_TOP, m.outer_faces_in_box([-1e-5,H-1e-5,-1e-5], 
[W+1e-5,H+1e-5,T+1e-5]))
+
+
+mfu_ = gf.MeshFem(m, 3)
+mfu_.set_classical_fem(fem_order)
+
+kept_dofs = list(range(mfu_.nbdof()))
+# remove x-dofs on RG_LEFT, y-dofs on RG_BOTTOM, and z-dofs on RG_BACK
+for skipped_dof, RG in enumerate((RG_LEFT,RG_BOTTOM,RG_BACK)):
+  for dof in mfu_.basic_dof_on_region(RG):
+    if dof % 3 == skipped_dof:
+      kept_dofs.remove(dof)
+mfu = gf.MeshFem('partial', mfu_, kept_dofs)
+
+mfmult = gf.MeshFem(m, 1)
+mfmult.set_classical_fem(fem_order)
+
+mfout = gf.MeshFem(m, 1)
+mfout.set_classical_discontinuous_fem(2)
+
+mim = gf.MeshIm(m, 5)
+
+# Model
+md = gf.Model('real')
+md.add_fem_variable('u', mfu)
+
+mimd3x3 = gf.MeshImData(mim, -1, [3,3])
+
+md.add_initialized_data('epsYY', 0)
+
+md.add_initialized_data('kappa', kappa)
+md.add_initialized_data('mu', mu)
+
+md.add_initialized_data('tauH', tauH)
+md.add_initialized_data('tauD', tauD)
+md.add_initialized_data('rH', rH)
+md.add_initialized_data('rD', rD)
+
+dt = t_max/steps
+md.add_initialized_data('dt', dt)
+
+md.add_im_data("SeHprev", mimd3x3)   # Elastic hydrostatic stress at previous 
timestep
+md.add_im_data("SeDprev", mimd3x3)   # Elastic deviatoric stress at previous 
timestep
+md.add_im_data("SvHprev", mimd3x3)   # Viscous hydrostatic stress at previous 
timestep
+md.add_im_data("SvDprev", mimd3x3)   # Viscous deviatoric stress at previous 
timestep
+
+md.add_macro("F", "Id(3)+Grad_u")
+md.add_macro("Fprev", "Id(3)+Grad_Previous_u")
+
+if MAT_LAW == 'neohookean-Simo':
+  md.add_macro("SeH", "kappa*log(Det(F))*Inv(Right_Cauchy_Green(F))")
+  md.add_macro("SeD", 
"mu*pow(Det(F),-2/3)*Inv(F)*Deviator(Left_Cauchy_Green(F))*Inv(F)'")
+elif MAT_LAW == 'Horgan-Murphy':
+  md.add_initialized_data('gamma', gamma)
+  md.add_macro("Sqr(a)", "a*a");
+  md.add_macro("WW1", "mu*(0.5+gamma)")
+  md.add_macro("WW2", "mu*(0.5-gamma)")
+  md.add_macro("WW3", "kappa*(1-1/Det(F)) - 2*WW2*pow(Det(F),-2/3) - 
WW1*pow(Det(F),-4/3)")
+  md.add_macro("I1", "Norm_sqr(F)")
+  md.add_macro("I2", "Sqr(Det(F))*Norm_sqr(Inv(F))")
+  md.add_macro("I3", "Sqr(Det(F))")
+  md.add_macro("SeH", "(1/3*WW1*I1 + 2/3*WW2*I2 + 
WW3*I3)*Inv(Right_Cauchy_Green(F))")
+  md.add_macro("SeD", "WW1*Id(3) + 
1/3*(WW2*I2-WW1*I1)*Inv(Right_Cauchy_Green(F))"
+                                "- WW2*I3*Sqr(Inv(Right_Cauchy_Green(F)))")
+
+md.add_macro("SvH", "exp(-dt/tauH)*SvHprev"
+                    "-(1-rH)*(tauH/dt-(dt+tauH)/dt*exp(-dt/tauH))*SeHprev"
+                    "-(1-rH)*((dt-tauH)/dt+tauH/dt*exp(-dt/tauH))*SeH")
+#                    "-(1-rH)*(1-exp(-dt/tauH))*SeHprev"
+#                    "-(1-rH)*(1-tauH/dt*(1-exp(-dt/tauH)))*(SeH-SeHprev)")
+md.add_macro("SvD", "exp(-dt/tauD)*SvDprev"
+                    "-(1-rD)*(tauD/dt-(dt+tauD)/dt*exp(-dt/tauD))*SeDprev"
+                    "-(1-rD)*((dt-tauD)/dt+tauD/dt*exp(-dt/tauD))*SeD")
+#                    "-(1-rD)*(1-exp(-dt/tauD))*SeDprev"
+#                    "-(1-rD)*(1-tauD/dt*(1-exp(-dt/tauD)))*(SeD-SeDprev)")
+
+# Virtual work expression
+md.add_nonlinear_generic_assembly_brick(mim, 
"(F*(SeH+SvH+SeD+SvD)):Grad_Test_u")
+
+# Cauchy and von Mises definitions
+md.add_macro("Cauchy", "(F*(SeH+SvH+SeD+SvD)*F')/Det(F)")
+md.add_macro("VM", "sqrt(1.5)*Norm(Deviator(F*(SeH+SvH+SeD+SvD)*F'))/Det(F)")
+
+# Dirichlet condition
+md.add_filtered_fem_variable("dirmult", mfmult, RG_TOP)
+md.add_linear_generic_assembly_brick(mim, "(epsYY*X(2)-u(2))*dirmult", RG_TOP)
+
+Vol = gf.asm_generic(mim, 0, "1", -1, md)
+with open("QLV_viscoelasticity_results.dat", "w") as f:
+  for i in range(steps+1):
+    t = i*dt
+    epsYY = 0.3*max(0,min(t/(2.*t_max/7.),((6.*t_max/7.)-t)/(2.*t_max/7.),1.))
+    md.set_variable("epsYY", epsYY)
+
+    md.solve('noisy', 'lsolver', 'mumps', 'max_res', 1E-7,'max_iter',20,
+             'lsearch', 'simplest', 'alpha max ratio', 1., 'alpha min', 0.2, 
'alpha mult', 0.6)
+
+    sigmayy = gf.asm_generic(mim, 0, "Cauchy(2,2)", -1, md) / Vol
+    f.write("t=%.10g sigma_yy=%.10g\n" % (t,sigmayy))
+
+    U = md.variable('u')
+    VM = md.local_projection(mim, "VM", mfout)
+    SIGMA11 = md.local_projection(mim, "Cauchy(1,1)", mfout)
+    SIGMA22 = md.local_projection(mim, "Cauchy(2,2)", mfout)
+    SIGMA33 = md.local_projection(mim, "Cauchy(3,3)", mfout)
+    mfout.export_to_vtk('QLV_viscoelasticity%i.vtk' % i, mfu, U, 
'Displacements',
+                                                     mfout, VM,  'Von Mises 
Stresses',
+                                                     mfout, SIGMA11,  'Sigma 
11',
+                                                     mfout, SIGMA22,  'Sigma 
22',
+                                                     mfout, SIGMA33,  'Sigma 
33')
+
+    md.set_variable("SeHprev", md.interpolation("SeH", mimd3x3, -1))
+    md.set_variable("SeDprev", md.interpolation("SeD", mimd3x3, -1))
+    md.set_variable("SvHprev", md.interpolation("SvH", mimd3x3, -1))
+    md.set_variable("SvDprev", md.interpolation("SvD", mimd3x3, -1))
+