Select Git revision
Forked from
gmsh / gmsh
Source project has a limited visibility.
cyl_LinearEMT.py 9.11 KiB
#coding-Utf-8-*-
from gmshpy import *
#from dG3DpyDebug import*
from dG3Dpy import*
#script to launch ElecMag SMP problem with a python script
import math
Irms = 0. # Ampere
freq = 0. # Hz
nTurnsCoil = 0
coilLx = coilLy = coilLz = 0. # m
coilW = 0. # m
mag_mu_0 = 4.e-7 * math.pi # vacuum magnetic permeability
tinitial = 273.+25.
# material law: SMP
lawnum= 1
rho = 1020. # material density
tinitial = 273.+25.
G=156.e6 # Shear modulus
Mu=0.35 # Poisson ratio
alpha = beta=gamma=0. # parameter of anisotropy
cp=1710.*rho # specific heat capacity
Kx=Ky=Kz=0.5 # Thermal conductivity
lx=ly=lz=35. # electrical conductivity (S/m)
seebeck=21.e-6
alphaTherm=0.0 #13.e-5 # thermal dilatation
v0=0. # initial electric potential
E_matrix = 2.0 * G * (1.0 + Mu)
mag_mu_0 = 4.e-7 * math.pi # vacuum magnetic permeability
mag_r_smp = 5.0 # relative mag permeability of SMP
mag_mu_x = mag_mu_y = mag_mu_z = mag_r_smp * mag_mu_0 # mag permeability smp
magpot0_x = magpot0_y = magpot0_z = 0.0 # initial magnetic potential
# material law: vacuum/free space region
lawnumvac = 3
rhovac = 1.2
Gvac=156.e1 # Shear modulus (1.e-5 less than SMP)
Muvac = 0.35 # Poisson ratio (same as SMP)
Evac= 2.0 * Gvac * (1.0 + Muvac) #youngs modulus
alphavac = betavac = gammavac = 0. # parameter of anisotropy
cpvac= 1012.*1.2
Kxvac=Kyvac=Kzvac= 26.e-3 # thermal conductivity tensor components
lxvac=lyvac=lzvac= 1.e-12 # electrical conductivity
seebeckvac= 0.
alphaThermvac= 0.0 # thermal dilatation
v0vac = 0.
mag_r_vac = 1.0 # relative mag permeability of inductor
mag_mu_x_vac = mag_mu_y_vac = mag_mu_z_vac = mag_r_vac * mag_mu_0 # mag permeability inductor
magpot0_x_vac = magpot0_y_vac = magpot0_z_vac = 0.0 # initial magnetic potential
useFluxT=True;
evaluateMecaField = False;
# geometry
geofile="cylinder.geo" # name of mesh file
meshfile="cylinder.msh" # name of mesh file
# solver
sol = 2 # Gmm=0 (default) Taucs=1 PETsc=2
soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
nstep = 5 # number of step (used only if soltype=1)
ftime = 5. # Final time (used only if soltype=1)
tol=1.e-3 # 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
beta1 = 1000.
eqRatio =1.e8
# compute solution and BC (given directly to the solver
# creation of law
# material law for smp
lawsmp = LinearElecMagTherMechDG3DMaterialLaw(lawnum,rho,E_matrix,E_matrix,E_matrix,Mu,Mu,Mu,G,G,G,alpha,beta,gamma,
tinitial,Kx,Ky,Kz,alphaTherm,alphaTherm,alphaTherm,lx,ly,lz,seebeck,cp,v0, mag_mu_x, mag_mu_y, mag_mu_z, magpot0_x,
magpot0_y, magpot0_z, Irms, freq, nTurnsCoil, coilLx, coilLy, coilLz, coilW,useFluxT)
#material law for free space (vacuum)
lawvac = LinearElecMagTherMechDG3DMaterialLaw(lawnumvac,rhovac,Evac,Evac,Evac,Muvac,Muvac,Muvac,Gvac,Gvac,Gvac,alphavac,betavac,gammavac,tinitial,
Kxvac,Kyvac,Kzvac,alphaThermvac,alphaThermvac,alphaThermvac,lxvac,lyvac,lzvac,seebeckvac,cpvac,v0vac, mag_mu_x_vac, mag_mu_y_vac, mag_mu_z_vac, magpot0_x_vac,
magpot0_y_vac, magpot0_z_vac, Irms, freq, nTurnsCoil, coilLx, coilLy, coilLz, coilW,useFluxT)
# creation of ElasticField
SMPfield = 10 # number of the field (physical number of SMP region)
Vacfield = 20 # number of the field (physical number of Vacuum)
SurfSMPStart = 3 # delta voltage surface
SurfSMPEnd = 4 # ground voltage surface
space1 = 0 # function space (Lagrange=0)
SMP_field = ElecMagTherMechDG3DDomain(10,SMPfield,space1,lawnum,fullDg,eqRatio,3)
Vacuum_field = ElecMagTherMechDG3DDomain(10,Vacfield,space1,lawnumvac,fullDg,eqRatio,3)
SMP_field.setConstitutiveExtraDofDiffusionEqRatio(eqRatio)
Vacuum_field.setConstitutiveExtraDofDiffusionEqRatio(eqRatio)
SMP_field.setConstitutiveCurlEqRatio(eqRatio)
Vacuum_field.setConstitutiveCurlEqRatio(eqRatio)
SMP_field.stabilityParameters(beta1)
Vacuum_field.stabilityParameters(beta1)
thermalSource=True
mecaSource=False
Vacuum_field.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
SMP_field.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
#SMP_field.matrixByPerturbation(1,1,1,1.e-3)
#Vacuum_field.matrixByPerturbation(1,1,1,1.e-3)
# flag to fix the reinitialization and
# reevaluation of field after restart
# args: boolean flag and field index number
SMP_field.setEvaluateMecaField(evaluateMecaField) # fix displacements and don't compute
Vacuum_field.setEvaluateMecaField(evaluateMecaField)
# creation of Solver
mysolver = nonLinearMechSolver(1000)
mysolver.createModel(geofile,meshfile,3,1)
#mysolver.loadModel(meshfile)
mysolver.addDomain(SMP_field)
mysolver.addDomain(Vacuum_field)
mysolver.addMaterialLaw(lawsmp)
mysolver.addMaterialLaw(lawvac)
mysolver.Scheme(soltype)
mysolver.Solver(sol)
mysolver.snlData(nstep,ftime,tol,tol)#forsmp
mysolver.options("-ksp_type preonly")
mysolver.options("-pc_type lu")
mysolver.stepBetweenArchiving(nstepArch)
mysolver.snlManageTimeStep(25, 3, 2, 10)
# BC
mysolver.displacementBC("Volume",SMPfield,0,0.0)
mysolver.displacementBC("Volume",SMPfield,1,0.0)
mysolver.displacementBC("Volume",SMPfield,2,0.0)
mysolver.displacementBC("Volume",Vacfield,0,0.0)
mysolver.displacementBC("Volume",Vacfield,1,0.0)
mysolver.displacementBC("Volume",Vacfield,2,0.0)
#thermal BC
cyclicFunctionTemp1=cycleFunctionTime(0., tinitial,ftime,tinitial);
mysolver.initialBC("Volume","Position",SMPfield,3,tinitial)
mysolver.initialBC("Volume","Position",Vacfield,3,tinitial)
mysolver.displacementBC("Face",8,3,cyclicFunctionTemp1)
#mysolver.displacementBC("Face",5,3,cyclicFunctionTemp1)
#mysolver.displacementBC("Face",6,3,cyclicFunctionTemp1)
#mysolver.displacementBC("Face",SurfSMPStart,3,cyclicFunctionTemp1)
#mysolver.displacementBC("Face",SurfSMPEnd,3,cyclicFunctionTemp1)
#mysolver.displacementBC("Volume",SMPfield,3,cyclicFunctionTemp1)
#mysolver.displacementBC("Volume",Vacfield,3,cyclicFunctionTemp1)
#electrical BC
mysolver.initialBC("Volume","Position",SMPfield,4,0.)
mysolver.initialBC("Volume","Position",Vacfield,4,0.)
cyclicFunctionvolt1=cycleFunctionTime(0., 0., ftime,1.e2);
mysolver.displacementBC("Face",SurfSMPStart,4,cyclicFunctionvolt1)
cyclicFunctionvolt2=cycleFunctionTime(0., 0., ftime,0.);
mysolver.displacementBC("Face",SurfSMPEnd,4,cyclicFunctionvolt2)
#magentic
mysolver.curlDisplacementBC("Face",SurfSMPStart,5,0.0) # comp may also be 5
mysolver.curlDisplacementBC("Face",SurfSMPEnd,5,0.0) # comp may also be 5
mysolver.curlDisplacementBC("Face",5,5,0.0) # comp may also be 5
mysolver.curlDisplacementBC("Face",6,5,0.0) # comp may also be 5
mysolver.curlDisplacementBC("Face",8,5,0.0) # comp may also be 5
#Gauging the edges using tree-cotree method
PhysicalCurves = "" # input required as a string of comma separated ints
PhysicalSurfaces = "3,4,5,6,8" # input required as a string of comma separated ints
PhysicalVolumes = "10,20" # input required as a string of comma separated ints
OutputPhysical = 55 # input required as a int
mysolver.createTreeForBC(PhysicalCurves,PhysicalSurfaces,PhysicalVolumes,OutputPhysical)
mysolver.curlDisplacementBC("Edge",OutputPhysical,5,0.0)
mysolver.internalPointBuildView("temperature",IPField.TEMPERATURE, 1, 1)
mysolver.internalPointBuildView("w_T",IPField.THERMALSOURCE, 1, 1)
mysolver.internalPointBuildView("jy_x",IPField.THERMALFLUX_X, 1, 1)
mysolver.internalPointBuildView("jy_y",IPField.THERMALFLUX_Y, 1, 1)
mysolver.internalPointBuildView("jy_z",IPField.THERMALFLUX_Z, 1, 1)
mysolver.internalPointBuildView("voltage",IPField.VOLTAGE, 1, 1)
mysolver.internalPointBuildView("je_x",IPField.ELECTRICALFLUX_X, 1, 1)
mysolver.internalPointBuildView("je_y",IPField.ELECTRICALFLUX_Y, 1, 1)
mysolver.internalPointBuildView("je_z",IPField.ELECTRICALFLUX_Z, 1, 1)
mysolver.internalPointBuildView("ax",IPField.MAGNETICVECTORPOTENTIAL_X, 1, 1)
mysolver.internalPointBuildView("ay",IPField.MAGNETICVECTORPOTENTIAL_Y, 1, 1)
mysolver.internalPointBuildView("az",IPField.MAGNETICVECTORPOTENTIAL_Z, 1, 1)
mysolver.internalPointBuildView("bx",IPField.MAGNETICINDUCTION_X, 1, 1)
mysolver.internalPointBuildView("by",IPField.MAGNETICINDUCTION_Y, 1, 1)
mysolver.internalPointBuildView("bz",IPField.MAGNETICINDUCTION_Z, 1, 1)
mysolver.internalPointBuildView("js0_x",IPField.INDUCTORSOURCEVECTORFIELD_X, 1, 1)
mysolver.internalPointBuildView("js0_y",IPField.INDUCTORSOURCEVECTORFIELD_Y, 1, 1)
mysolver.internalPointBuildView("js0_z",IPField.INDUCTORSOURCEVECTORFIELD_Z, 1, 1)
mysolver.internalPointBuildView("wAV_x",IPField.EMSOURCEVECTORFIELD_X, 1, 1)
mysolver.internalPointBuildView("wAV_y",IPField.EMSOURCEVECTORFIELD_Y, 1, 1)
mysolver.internalPointBuildView("wAV_z",IPField.EMSOURCEVECTORFIELD_Z, 1, 1)
mysolver.archivingForceOnPhysicalGroup("Face", SurfSMPStart, 3, 1);
mysolver.archivingForceOnPhysicalGroup("Face", SurfSMPStart, 4, 1);
mysolver.archivingForceOnPhysicalGroup("Face", SurfSMPEnd, 4, 1);
mysolver.archivingForceOnPhysicalGroup("Face", 7, 5, 1);
mysolver.solve()
check = TestCheck()
check.equal(0.000000e+00,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfSMPStart, 3),1.e-5)
check.equal(-1.273695e+09,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfSMPStart, 4),1.e-5)
check.equal(1.273695e+09,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfSMPEnd, 4),1.e-5)
check.equal(1.387338e+08,mysolver.getArchivedForceOnPhysicalGroup("Face", 7, 5),1.e-5)