Select Git revision
Options.cpp
Forked from
gmsh / gmsh
Source project has a limited visibility.
-
Christophe Geuzaine authoredChristophe Geuzaine authored
Plane_notch.py 7.56 KiB
#coding-Utf-8-*-
from gmshpy import*
#from dG3DpyDebug import*
from dG3Dpy import*
# Script for plane notched specimen with Gurson
# Material law creation
# ===================================================================================
BulkLawNum1 = 1
# material parameters - bulk law
rho = 7600. # kg/m3: density
young = 210.e3 # Pa: module de young
nu = 0.3 # -: poisson coefficient
sy0 = 354.41 # initial yield stress
h = 100.
p1 = 0.05
n1 = 0.2
p2 = 0.1
n2 = 0.15
# R = yield0 + H*p if p < p1
# R = syp1*(p/p1)^n1 p1 < p < p2
# R = syp2*(p/p2)^n2 p2 < p
harden_Law1 = LinearFollowedByMultiplePowerLawJ2IsotropicHardening(BulkLawNum1,sy0,h,p1,n1,p2,n2)
# - gurson
q1 = 1.5
q2 = 1.
q3 = 1.5
fVinitial = 1.5e-11 #cannot be 0
# - coalesence
fC = 0.0#0.1#0.005
ff = 1./q1#0.3#0.225
ffstar = 1./q1
# - nucleation
function1 = LinearNucleationFunctionLaw(0, 0.2)
criterion1 = PorosityBasedNucleationCriterionLaw(0, 0.005, 1.)
# - plastic instabilities
eps0=0.1
deltaY=0.2
instability_Law1 = ExponentialPlasticInstabilityFunctionLaw(1, eps0, deltaY)
# yield = (1-wp)*R
# wp = deltaY_*(1.-exp(-hatpm/eps0_))
# - non-local law
l = 0.25
cl = l*l
ClLawNum1 = 4
Cl_Law1 = IsotropicCLengthLaw(ClLawNum1, cl)
# material law creation
BulkLaw1 = NonLocalDamageGursonDG3DMaterialLaw(BulkLawNum1, rho, young, nu, q1,q2,q3,fVinitial, harden_Law1, Cl_Law1,1e-8,False,1e-8)
BulkLaw1.setOrderForLogExp(9)
#BulkLaw1.setNucleationLaw(function1, criterion1)
BulkLaw1.setPlasticInstabilityFunctionLaw(instability_Law1)
BulkLaw1.setSubStepping(True,3)
BulkLaw1.setNonLocalMethod(2)
BulkLaw1.clearCLengthLaw()
BulkLaw1.setCLengthLaw(Cl_Law1)
BulkLaw1.setCLengthLaw(Cl_Law1)
BulkLaw1.setCLengthLaw(Cl_Law1)
BulkLaw1.useI1J2J3Implementation(True)
# Solver parameters - shift
# ===================================================================================
ftime =1.0 # Final time
r_impl = 0.2 #1.0
# Solver parameters - implicit
# ===================================================================================
soltype = 1 # StaticLinear=0 (default, StaticNonLinear=1, Explicit=2, # Multi=3, Implicit=4, Eigen=5)
nstep = 100 # Number of step
tol=1.e-6 # Relative tolerance for NR scheme
tolAbs = 1.e-8 # Absolute tolerance for NR scheme
nstepArchIP=1 # Number of step between 2 archiving
nstepArchForce = 1 # Number of step between 2 force archiving
nstepArchEnergy = nstepArchForce # Number of step between 2 energy computation and archiving
MaxIter = 15 # Maximum number of iterations
StepIncrease = 2 # Number of successfull timestep before reincreasing it
StepReducFactor = 2.0 # Timestep reduction factor
NumberReduction = 10 # Maximum number of timespep reduction: max reduction = pow(StepReducFactor,this)
fullDg = bool(0) # O = CG, 1 = DG
dgnl = bool(0) # DG for non-local variables inside a domain (only if fullDg)
eqRatio = 1.0e3 # Ratio between "elastic" and non-local equations
space1 = 0 # Function space (Lagrange=0)
beta1 = 30.0 # Penality parameter for DG
# Domain creation
## ===================================================================================
numPhysVolume1 = 700 # Number of a physical volume of the model in .geo
numDomain1 = 1001 # Number of the domain
field1 = dG3DDomain(numDomain1,numPhysVolume1,space1,BulkLawNum1,fullDg,2,3)
#field1.stabilityParameters(beta1) # Adding stability parameters (for DG)
#field1.setNonLocalStabilityParameters(beta1,dgnl) # Adding stability parameters (for DG)
#field1.gaussIntegration(0,-1,-1)
#field1.matrixByPerturbation(1,1,1,1e-8) # Tangent computation analytically or by pertubation
#field1.strainSubstep(2,10)
field1.setNonLocalEqRatio(eqRatio)
# Solver creation
# ===================================================================================
mysolver = nonLinearMechSolver(1002) # Solver associated with numSolver1
meshfile= "plane_notch.msh" # name of mesh file
geofile= "plane_notch.geo" # name of geo file
#mysolver.loadModel(meshfile) # add mesh
mysolver.createModel(geofile,meshfile,3,2)
mysolver.addDomain(field1) # add domain
mysolver.addMaterialLaw(BulkLaw1) # add material law
# solver parametrisation
mysolver.Scheme(soltype) # solver scheme
mysolver.Solver(2) # Library solver: Gmm=0 (default) Taucs=1 PETsc=2
mysolver.snlData(nstep,ftime*r_impl,tol,tolAbs) # solver parameters for imp. solving (numb. of step, final time and tolerance)
mysolver.snlManageTimeStep(MaxIter,StepIncrease,StepReducFactor,NumberReduction)
# solver archiving
mysolver.stepBetweenArchiving(nstepArchIP) # archiving frequency
mysolver.energyComputation(nstepArchEnergy) # archiving frequency for energy
#mysolver.lineSearch(bool(0)) # lineSearch activation
# path following
method = 0
withPF = True
mysolver.pathFollowing(withPF,method)
if method == 0:
mysolver.setPathFollowingIncrementAdaptation(True,4)
mysolver.setPathFollowingControlType(0)
mysolver.setPathFollowingCorrectionMethod(0)
mysolver.setPathFollowingArcLengthStep(1e-1)
mysolver.setBoundsOfPathFollowingArcLengthSteps(0,1.);
else:
# time-step adaptation by number of NR iterations
mysolver.setPathFollowingIncrementAdaptation(True,4)
mysolver.setPathFollowingLocalSteps(1e-2,1.)
mysolver.setPathFollowingSwitchCriterion(0.)
mysolver.setPathFollowingLocalIncrementType(1);
mysolver.setBoundsOfPathFollowingLocalSteps(1.,3.)
# Boundary conditions
# ===============================
mysolver.displacementBC("Face",700,2,0.)
mysolver.displacementBC("Node",401,1,0.)
mysolver.displacementBC("Edge",1400,0,0.)
if withPF == True:
mysolver.sameDisplacementBC("Edge",400,401,0)
mysolver.forceBC("Node",401,0,1e4)
else:
mysolver.displacementBC("Edge",400,0,3.)
# Variable storage
# ===============================
mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
mysolver.internalPointBuildView("damage",IPField.DAMAGE,1,1)
mysolver.internalPointBuildView("local_fV",IPField.LOCAL_POROSITY,1,1)
mysolver.internalPointBuildView("local_fV_max",IPField.LOCAL_POROSITY,1,IPField.MAX_VALUE)
mysolver.internalPointBuildView("epl",IPField.PLASTICSTRAIN, 1, 1)
mysolver.internalPointBuildView("triaxiality",IPField.STRESS_TRIAXIALITY, 1, 1)
mysolver.internalPointBuildView("LOCAL_0",IPField.LOCAL_0, 1, 1)
mysolver.internalPointBuildView("NONLOCAL_0",IPField.NONLOCAL_0, 1, 1)
mysolver.internalPointBuildView("LOCAL_1",IPField.LOCAL_1, 1, 1)
mysolver.internalPointBuildView("NONLOCAL_1",IPField.NONLOCAL_1, 1, 1)
mysolver.internalPointBuildView("LOCAL_2",IPField.LOCAL_2, 1, 1)
mysolver.internalPointBuildView("NONLOCAL_2",IPField.NONLOCAL_2, 1, 1)
mysolver.internalPointBuildView("LOCAL_POROSITY",IPField.LOCAL_POROSITY,1,1)
mysolver.internalPointBuildView("CORRECTED_POROSITY",IPField.CORRECTED_POROSITY,1,1)
mysolver.internalPointBuildView("YIELD_POROSITY",IPField.YIELD_POROSITY,1,1)
mysolver.archivingForceOnPhysicalGroup("Edge", 1400, 0, nstepArchForce)
mysolver.archivingNodeDisplacement(401,0, nstepArchForce)
mysolver.archivingNodeDisplacementOnPhysical(0,100,1,nstepArchForce)
mysolver.archivingIPOnPhysicalGroup("Face", 700, IPField.DAMAGE,IPField.MAX_VALUE,nstepArchForce)
mysolver.archivingIPOnPhysicalGroup("Face", 700, IPField.LOCAL_POROSITY,IPField.MAX_VALUE,nstepArchForce)
mysolver.solve()
check = TestCheck()
check.equal(-1.957854e+03,mysolver.getArchivedForceOnPhysicalGroup("Edge", 1400, 0),1.e-6)