Old test for thermal evolution in SMP domain with analytical EM thermal source...

Old test for thermal evolution in SMP domain with analytical EM thermal source wEM computed using values from weak coupling test for same EM input parameters

dG3D/benchmarks/ThermoMechanicsInSMP/CMakeLists.txt

+# test file
+
+set(PYFILE ThermoMechanicsInSMP.py)
+
+set(FILES2DELETE 
+  *.msh
+  *.csv
+)
+
+add_cm3python_test(${PYFILE} "${FILES2DELETE}")

dG3D/benchmarks/ThermoMechanicsInSMP/ThermoMechanicsInSMP.py

+#coding-Utf-8-*-
+from gmshpy import *
+
+
+#from dG3DpyDebug import*
+from dG3Dpy import*
+
+#script to launch unit test to debug strong vs weak coupling in SMP
+
+# Here strong coupled EM-TM in SMP only with analytically computed EM thermal source
+# No solution of EM fields but directly EM thermal source is imported to 
+# resolve TM in SMP
+
+# Need to uncomment in mlawLinearEMTM
+# wAV = A² cos² (2 pi f t)
+# with A²/2 = mean computed over one period for considered f, I from weak coupling test
+
+import math
+
+# material law: SMP
+lawnum= 1 
+rho   = 270. # material density (Kg/m^3)
+tinitial = 273.+47. # (K)
+G=299.e6 # Shear modulus (Pa)
+Mu=0.29 # Poisson ratio
+alpha = beta=gamma=0. # parameter of anisotropy
+cp=10.*rho # specific heat capacity 
+Kx=Ky=Kz=237. # Thermal conductivity (W/m K)
+lx=ly=lz=1.0e4 # electrical conductivity (S/m) 
+seebeck=0.0 #21.e-6
+alphaTherm=0.0 #13.e-5 # thermal dilatation
+v0=0. # initial electric potential
+E_matrix = 771.e6 # (Pa)
+mag_mu_0 = 4.e-7 * math.pi # vacuum magnetic permeability
+mag_r_smp = 20.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
+Irms = 325. # (Ampere)
+freq = 1000. # (Hz)
+nTurnsCoil = 1000
+coilLx = coilLy = coilLz = 1. # m
+coilW = 3.e-2 # m
+
+useFluxT=True;
+evaluateCurlField = True;
+evaluateVoltage = True;
+
+# geometry
+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 = 160 # 8*EMncycles # number of step (used only if soltype=1)
+ftime = 5.e-2 # EMncycles/freq; # according to characteristic time
+tol=1.e-6   # 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.e0
+
+#  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,evaluateCurlField)
+
+#lawsmp.setUseBarF(True)
+
+# creation of ElasticField
+SMPfield = 1000 # number of the field (physical number of SMP)
+SurfSmp = 11110 # outer surface of conductor SMP
+SmpRef=11111 # top face of outer surface of SMP
+SmpRefBottom=11115 # bottom face of outer surface of SMP
+SmpRefPoint=11112 # single point on top face of SMP
+
+space1 = 0 # function space (Lagrange=0)
+
+SMP_field = ElecMagTherMechDG3DDomain(10,SMPfield,space1,lawnum,fullDg,eqRatio,3)
+
+SMP_field.setConstitutiveExtraDofDiffusionEqRatio(eqRatio)
+SMP_field.setConstitutiveCurlEqRatio(eqRatio)
+SMP_field.stabilityParameters(beta1)
+
+# flag to fix the reinitialization and 
+# reevaluation of field after restart
+# args: boolean flag and field index number
+SMP_field.setEvaluateCurlField(evaluateCurlField,0) # fix curl field and don't compute
+SMP_field.setEvaluateExtraDofDiffusionField(evaluateVoltage,1) # fix voltage and don't compute
+
+thermalSource=True
+mecaSource=False
+
+SMP_field.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
+
+#SMP_field.matrixByPerturbation(1,1,1,1.e-3)
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(SMP_field)
+mysolver.addMaterialLaw(lawsmp)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol,1.e-9)#forsmp
+mysolver.options("-ksp_type preonly")
+mysolver.options("-pc_type lu")
+mysolver.stepBetweenArchiving(nstepArch)
+mysolver.snlManageTimeStep(25, 10, 2, 10)
+
+# BC
+finalDisp = 0.01;
+cyclicFunctionDisp=cycleFunctionTime(0., 0., ftime , finalDisp);
+mysolver.displacementBC("Face",SmpRefBottom,0,0.)
+mysolver.displacementBC("Face",SmpRefBottom,1,0.)
+mysolver.displacementBC("Face",SmpRefBottom,2,0.)
+mysolver.displacementBC("Face",SmpRef,0,0.)
+mysolver.displacementBC("Face",SmpRef,1,cyclicFunctionDisp)
+mysolver.displacementBC("Face",SmpRef,2,0.)
+
+#thermal BC
+mysolver.initialBC("Volume","Position",SMPfield,3,tinitial)
+
+#electrical BC
+mysolver.initialBC("Volume","Position",SMPfield,4,0.0)
+mysolver.displacementBC("Node",SmpRefPoint,4,0.0)
+
+#magentic
+mysolver.curlDisplacementBC("Volume",SMPfield,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 = "11110" # input required as a string of comma separated ints
+PhysicalVolumes = "1000" 	# 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("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("eps_xx",IPField.STRAIN_XX, 1, 1)
+mysolver.internalPointBuildView("eps_yy",IPField.STRAIN_YY, 1, 1)
+mysolver.internalPointBuildView("eps_zz",IPField.STRAIN_ZZ, 1, 1)
+mysolver.internalPointBuildView("eps_xy",IPField.STRAIN_XY, 1, 1)
+mysolver.internalPointBuildView("eps_yz",IPField.STRAIN_YZ, 1, 1)
+mysolver.internalPointBuildView("eps_xz",IPField.STRAIN_XZ, 1, 1)
+mysolver.internalPointBuildView("temperature",IPField.TEMPERATURE, 1, 1)
+mysolver.internalPointBuildView("w_AV",IPField.EMFIELDSOURCE, 1, 1)
+mysolver.internalPointBuildView("w_T",IPField.THERMALSOURCE, 1, 1)
+mysolver.internalPointBuildView("jyx",IPField.THERMALFLUX_X, 1, 1)
+mysolver.internalPointBuildView("jyy",IPField.THERMALFLUX_Y, 1, 1)
+mysolver.internalPointBuildView("jyz",IPField.THERMALFLUX_Z, 1, 1)
+mysolver.internalPointBuildView("voltage",IPField.VOLTAGE, 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("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("jex",IPField.ELECTRICALFLUX_X, 1, 1)
+mysolver.internalPointBuildView("jey",IPField.ELECTRICALFLUX_Y, 1, 1)
+mysolver.internalPointBuildView("jez",IPField.ELECTRICALFLUX_Z, 1, 1)
+
+mysolver.archivingIPOnPhysicalGroup("Volume",SMPfield, IPField.TEMPERATURE,IPField.MEAN_VALUE,nstepArch);
+mysolver.archivingIPOnPhysicalGroup("Volume",SMPfield, IPField.VOLTAGE,IPField.MEAN_VALUE,nstepArch);
+mysolver.archivingIPOnPhysicalGroup("Volume",SMPfield, IPField.EMFIELDSOURCE,IPField.MEAN_VALUE,nstepArch);
+
+mysolver.archivingForceOnPhysicalGroup("Face", SurfSmp, 3, 1);
+mysolver.archivingForceOnPhysicalGroup("Face", SurfSmp, 4, 1);
+mysolver.archivingForceOnPhysicalGroup("Face", SurfSmp, 5, 1);
+
+mysolver.archivingNodeDisplacement(SmpRefPoint,0,1);
+mysolver.archivingNodeDisplacement(SmpRefPoint,1,1);
+mysolver.archivingNodeDisplacement(SmpRefPoint,2,1);
+mysolver.archivingNodeDisplacement(SmpRefPoint,3,1);
+mysolver.archivingNodeDisplacement(SmpRefPoint,4,1);
+
+#mysolver.setWriteDeformedMeshToFile(True);
+
+mysolver.solve()
+
+"""
+check = TestCheck()
+check.equal(0.000000e+00,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfSmp, 3),1.e-5)
+check.equal(0.000000e+00,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfSmp, 4),1.e-5)
+check.equal(1.427543e-01,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfSmp, 5),1.e-5)
+check.equal(0.0,mysolver.getArchivedNodalValue(SmpRefPoint,0,mysolver.displacement),1.e-5)
+check.equal(0.0,mysolver.getArchivedNodalValue(SmpRefPoint,1,mysolver.displacement),1.e-5)
+check.equal(0.0,mysolver.getArchivedNodalValue(SmpRefPoint,2,mysolver.displacement),1.e-5)
+check.equal(2.997321e+02,mysolver.getArchivedNodalValue(SmpRefPoint,3,mysolver.displacement),1.e-5)
+check.equal(3.548869e+03,mysolver.getArchivedNodalValue(SmpRefPoint,4,mysolver.displacement),1.e-5)
+"""
+

[Vinayak Gholap]

Note: When running this test, follow these comments.

[Vinayak Gholap]

Note: here for the test using weak coupling scheme, get magnitude needed for wEM from weak coupling. Make necessary changes in mlaw to use the analytical wEM instead of computing from EM problem.