Select Git revision
runTVP.py 8.35 KiB
#coding-Utf-8-*-
from gmshpy import *
from dG3Dpy import*
#from dG3DpyDebug import*
from math import*
import csv
import numpy as np
# Load the csv data for relaxation spectrum
with open('relaxationSpectrumBetter_N30.csv', newline='') as csvfile:
reader = csv.reader(csvfile, delimiter=';')
relSpec = np.zeros((1,))
i = 0
for row in reader:
if i == 0:
relSpec[i] = float(' '.join(row))
else:
relSpec = np.append(relSpec, float(' '.join(row)))
i += 1
# print(relSpec)
N = int(0.5*(i-1))
relSpec[0:N+1] *= 1e-6 # convert to MPa
#script to launch beam problem with a python script
# material law
lawnum = 11 # unique number of law
E = relSpec[0] #MPa
nu = 0.33
K = E/3./(1.-2.*nu) # Bulk mudulus
mu = E/2./(1.+nu) # Shear mudulus
rho = 905e-12 # Bulk mass 905 kg/m3 -> 905e-12 tonne/mm3
Cp = rho*1900e+6 # 1900 J/kg/K -> 1900e+6 Nmm/tonne/K
KThCon = 0.14 # 0.14 W/m/K -> 0.14 Nmm/s/mm/K
Alpha = 0.6e-6 # 1/K
Tinitial = 273.15 + 23 # K
C1 = 136.17201827
C2 = 256.58268053 # K
# Temp Settings
Tref = 273.15+20.
TemFuncOpt = 0 # 0-> Shift Factor based derivatives of Ki; 1-> Gibson/Huang based derivatives
# Default tempfuncs
E_G = 2.80717267e+09 * 1e-6
E_R = 5.52466562e+08 * 1e-6
Cp_G = rho*1700e+6
Cp_R = rho*2100e+6
Alpha_G = 0.2e-6
Alpha_R = 1.e-6
KThCon_G = 0.10
KThCon_R = 1.5*0.12
k_gibson = 3.36435338e-02
m_gibson = 0.0
Tg = Tref
E_tempfunc = GlassTransitionScalarFunction(1, E_R/E_G, Tg, k_gibson, m_gibson, 0) # can't be used
Alpha_tempfunc = GlassTransitionScalarFunction(1, Alpha_R/Alpha_G, Tg, k_gibson, m_gibson, 0)
Cp_tempfunc = GlassTransitionScalarFunction(1, Cp_R/Cp_G, Tg, k_gibson, m_gibson, 0)
KThCon_tempfunc = GlassTransitionScalarFunction(1, KThCon_R/KThCon_G, Tg, k_gibson, m_gibson, 0)
ShiftFactor_tempfunc = negExpWLFshiftFactor(C1,C2,Tref)
sy0c = 10./0.78 #MPa
hc = 30.
sy0t = sy0c*0.78
ht = 30.
hardenc = LinearExponentialJ2IsotropicHardening(1, sy0c, hc, 10., 100.)
hardent = LinearExponentialJ2IsotropicHardening(2, sy0t, ht, 10., 100.)
# hardenc.setTemperatureFunction_h1(ShiftFactor_tempfunc)
# hardenc.setTemperatureFunction_h2(ShiftFactor_tempfunc)
# hardenc.setTemperatureFunction_hexp(ShiftFactor_tempfunc)
# hardent.setTemperatureFunction_h1(ShiftFactor_tempfunc)
# hardent.setTemperatureFunction_h2(ShiftFactor_tempfunc)
# hardent.setTemperatureFunction_hexp(ShiftFactor_tempfunc)
hardenk = PolynomialKinematicHardening(3,1)
hardenk.setCoefficients(1,30.)
# hardenk.setTemperatureFunction_K(ShiftFactor_tempfunc)
law1 = NonLinearTVPDG3DMaterialLaw(lawnum,rho,E,nu,1e-6,Tinitial,Alpha,KThCon,Cp) #,True,1e-8)
law1.setCompressionHardening(hardenc)
law1.setTractionHardening(hardent)
law1.setKinematicHardening(hardenk)
law1.setStrainOrder(11)
law1.setYieldPowerFactor(3.5)
law1.setNonAssociatedFlow(True)
law1.nonAssociatedFlowRuleFactor(0.5)
law1.setViscoelasticMethod(0)
law1.setTestOption(3)
law1.setShiftFactorConstantsWLF(C1,C2)
law1.setReferenceTemperature(Tref)
law1.setTempFuncOption(0)
law1.setViscoElasticNumberOfElement(N)
if N > 0:
for i in range(1, N+1):
law1.setViscoElasticData(i, relSpec[i], relSpec[i+N])
law1.setTemperatureFunction_ThermalDilationCoefficient(Alpha_tempfunc)
law1.setTemperatureFunction_ThermalConductivity(KThCon_tempfunc)
# law1.setTemperatureFunction_SpecificHeat(Cp_tempfunc)
eta = constantViscosityLaw(1,3e4)
law1.setViscosityEffect(eta,0.21)
# eta.setTemperatureFunction(ShiftFactor_tempfunc)
law1.setIsotropicHardeningCoefficients(1.,1.,1.)
# geometry
geofile = "line.geo"
meshfile = "line.msh"
# solver
sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
nstep = 100 # number of step (used only if soltype=1)
ftime = 1 # Final time (used only if soltype=1)
tol=1.e-5 # relative tolerance for NR scheme (used only if soltype=1)
nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
fullDg = 0 #O = CG, 1 = DG
space1 = 0 # function space (Lagrange=0)
beta1 = 100
# creation of ElasticField
#matrix
#myfield1 = dG3D1DDomain(11,11,space1,lawnum,0,1) #for non local
# myfield1 = dG3D1DDomain(11,11,space1,lawnum, False) #for non local
#myfield1.matrixByPerturbation(1,1,1,1e-8)
ThermoMechanicsEqRatio = 1.e1 # setConstitutiveExtraDofDiffusionEqRatio(ThermoMechanicsEqRatio); WHAT is this parameter?
thermalSource = True
mecaSource = True
# myfield1 = ThermoMechanicsDG3DDomain(1000,11,space1,lawnum,fullDg,ThermoMechanicsEqRatio,1) # Added
myfield1 = dG3D1DDomain(11,11,space1,lawnum,0,0,1) #for non local
myfield1.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
myfield1.stabilityParameters(beta1)
myfield1.ThermoMechanicsStabilityParameters(beta1,bool(1)) # Added
myfield1.setThermoMechanicsEqRatio(ThermoMechanicsEqRatio)
myfield1.setState(1) #UniaxialStrain=0, UniaxialStress =1, ConstantTriaxiality=2
L = 1.
x0 = 0;
A0 = 0.99
xL = L
AL = 1
a = (AL - A0)/(xL-x0)
crossSection = linearScalarFunction(x0,A0,a)
myfield1.setCrossSectionDistribution(crossSection)
# creation of Solver
mysolver = nonLinearMechSolver(1000)
mysolver.createModel(geofile,meshfile,1,1)
mysolver.addDomain(myfield1)
mysolver.addMaterialLaw(law1)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol)
mysolver.stepBetweenArchiving(nstepArch)
mysolver.snlManageTimeStep(15,5,2.,20) # maximal 20 times
mysolver.displacementBC("Edge",11,1,0.)
mysolver.displacementBC("Edge",11,2,0.)
mysolver.displacementBC("Node",1,0,0.)
method = 0
withPF = False #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,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.)
if withPF:
mysolver.forceBC("Node",2,0,1e2)
else:
# mysolver.displacementBC("Node",2,0,7e-2*L)
mysolver.displacementBC("Node",2,0, 0.0834) # strain rate = 8.333e-2
# mysolver.displacementBC("Node",2,0,7e-2*L)
mysolver.initialBC("Volume","Position",11,3,Tinitial)
fT = LinearFunctionTime(0,Tinitial,ftime,Tinitial); # Linear Displacement with time
# mysolver.displacementBC("Volume",11,3,fT)
"""
mysolver.internalPointBuildView("P_xx",IPField.P_XX, 1, 1)
mysolver.internalPointBuildView("P_yy",IPField.P_YY, 1, 1)
mysolver.internalPointBuildView("P_zz",IPField.P_ZZ, 1, 1)
mysolver.internalPointBuildView("P_xy",IPField.P_XY, 1, 1)
mysolver.internalPointBuildView("P_yx",IPField.P_YX, 1, 1)
mysolver.internalPointBuildView("P_yz",IPField.P_YZ, 1, 1)
mysolver.internalPointBuildView("P_zy",IPField.P_ZY, 1, 1)
mysolver.internalPointBuildView("P_xz",IPField.P_XZ, 1, 1)
mysolver.internalPointBuildView("P_zx",IPField.P_ZX, 1, 1)
mysolver.internalPointBuildView("F_xx",IPField.F_XX, 1, 1)
mysolver.internalPointBuildView("F_yy",IPField.F_YY, 1, 1)
mysolver.internalPointBuildView("F_zz",IPField.F_ZZ, 1, 1)
mysolver.internalPointBuildView("F_xy",IPField.F_XY, 1, 1)
mysolver.internalPointBuildView("F_yx",IPField.F_YX, 1, 1)
mysolver.internalPointBuildView("F_yz",IPField.F_YZ, 1, 1)
mysolver.internalPointBuildView("F_zy",IPField.F_ZY, 1, 1)
mysolver.internalPointBuildView("F_xz",IPField.F_XZ, 1, 1)
mysolver.internalPointBuildView("F_zx",IPField.F_ZX, 1, 1)
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("epl",IPField.PLASTICSTRAIN, 1, 1)
"""
mysolver.archivingForceOnPhysicalGroup("Node", 1, 0)
mysolver.archivingNodeDisplacement(2,0)
# mysolver.archivingAverageValue(IPField.P_YY)
# mysolver.archivingAverageValue(IPField.F_YY)
mysolver.archivingAverageValue(IPField.P_XX)
mysolver.archivingAverageValue(IPField.F_XX)
mysolver.solve()
check = TestCheck()
check.equal(-3.089706e+01,mysolver.getArchivedForceOnPhysicalGroup("Node", 1, 0),1.e-4)