Skip to content
Snippets Groups Projects
Commit 99064080 authored by Ludovic Noels's avatar Ludovic Noels
Browse files

Merge branch 'vinayak' into 'master'

Updates in strong coupling tests (change cpvac, Kvac to get right delatT in...

See merge request !419
parents e413179d cf38f12d
No related branches found
No related tags found
1 merge request!419Updates in strong coupling tests (change cpvac, Kvac to get right delatT in...
Showing
with 179076 additions and 55264 deletions
......@@ -290,7 +290,6 @@ add_subdirectory(Cube_EM_SMPPhenomenological_Mech)
add_subdirectory(Cube_EM_SMPPhenomenological_Thermo)
add_subdirectory(EMTM_SMPPhenomenological_Mech)
add_subdirectory(EM_heating_SMPPheno_withAir)
add_subdirectory(strongCoupling_SMPPhenoMech_withAir)
add_subdirectory(steinmann_MTM_Benchmarks)
add_subdirectory(steinmann_ETM_Benchmarks)
add_subdirectory(air_wire_smp)
......
......@@ -10,20 +10,26 @@ from dG3Dpy import*
import math
#==========Electric & magnetic parameters SMPC========================
rho = 1200.
lx=ly=lz=1.0e5 # electrical conductivity (S/m)
seebeck=0.0 #21.e-6
v0=0. # initial electric potential
epsilon_0 = 8.854 * 1.e-12 # vacuum electric permittivity
mag_mu_0 = 4.e-7 * math.pi # vacuum magnetic permeability
epsilon_r_smp = 1.0 # relative elec permittivity of SMP
mag_r_smp = 20.0 # relative mag permeability of SMP
mag_r_smp = 5.0 # relative mag permeability of SMP
# use of electric permittivity constants to formulate polarization
epsilon_4 = 0.0
epsilon_6 = 0.0
epsilon_5 = 0.0 #(epsilon_0 * (1.0-epsilon_r_smp)/rho)
# use of magnetic permeability constants to formulate magnetization
mu_7 = 0.0 # magnetic permeability constants of SMP
mu_8 = -1.58733e-3 # magnetic permeability constants of SMP (rho * mag_mu_0)/((1.0/mag_r_smp)-1.0)
if(mag_r_smp == 1.0):
mu_8 = 0.0
else:
mu_8 = (rho * mag_mu_0)/((1.0/mag_r_smp)-1.0) # magnetic permeability constants of SMP
#print("mu_8 = ",mu_8)
mu_9 = 0.0 # magnetic permeability constants 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
......@@ -43,7 +49,7 @@ lawnum= 1
#=====SMP Phenomenological Thermo-Mech=============================================
#crystalline phase#
alpha=beta=gamma= 0.
rho = 1200.
#rho = 1200.
Eg = 1250.e6
Er = 9.7e6 #3.2e6
nug = 0.26 #0.4
......@@ -221,7 +227,7 @@ lawSMPPheno= PhenomenologicalSMPDG3DMaterialLaw(lawnum+7, rho, alpha, beta, gamm
TaylorQuineyG,TaylorQuineyR,TaylorQuineyAM, zAM)
lawSMPPheno.setFunctionCrystallizationVolumeRate(crystallizationRate)
lawSMPPheno.setStrainOrder(5)
lawSMPPheno.setStrainOrder(-1)
lawSMPPheno.setJ2IsotropicHardeningTractionG(_j2IHGt)
lawSMPPheno.setJ2IsotropicHardeningCompressionG(_j2IHGc)
lawSMPPheno.setKinematicHardeningG(hardenkG)
......@@ -307,8 +313,8 @@ Gvac=156.e1 # Shear modulus
Muvac = 0.35 # Poisson ratio
Evac= 2.0 * Gvac * (1.0 + Muvac) #youngs modulus
alphavac = betavac = gammavac = 0. # parameter of anisotropy
cpvac= 1012.*rhovac
Kxvac=Kyvac=Kzvac= 1.e-12 #26.e-3 # thermal conductivity tensor components
cpvac= 1.e-12 #1012.*rhovac
Kxvac=Kyvac=Kzvac= 0.0 #1.e-12 #26.e-3 # thermal conductivity tensor components
lxvac=lyvac=lzvac= 1.e-12 # electrical conductivity
seebeckvac= 0.
alphaThermvac= 0.0 # thermal dilatation
......
import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
Forcefilename = "/force11111comp1.csv"
Displacementfilename = "/NodalDisplacementPhysical11112Num1comp1.csv"
allDir = os.listdir("./TMresults/")
for d in allDir[:]:
temppath = "./TMresults/"+d
if os.path.isdir(temppath) == False:
allDir.remove(d)
if os.path.isfile(d) == True:
allDir.remove(d)
allDir.sort()
print(allDir)
force = "./TMresults/"+allDir[0]+Forcefilename
displ = "./TMresults/"+allDir[0]+Displacementfilename
newdataframeForce = pd.read_csv(force,sep=";",header=None).values
print(newdataframeForce.shape)
print()
newdataframeDispl = pd.read_csv(displ,sep=";",header=None).values
print(newdataframeDispl.shape)
print()
for d in allDir[1:]:
force = "./TMresults/"+d+Forcefilename
tempframeForce = pd.read_csv(force,sep=";",header=None).values
print(tempframeForce.shape)
print()
newdataframeForce = np.concatenate((newdataframeForce, tempframeForce), axis=0)
displ = "./TMresults/"+d+Displacementfilename
tempframeDispl = pd.read_csv(displ,sep=";",header=None).values
print(tempframeDispl.shape)
print()
newdataframeDispl = np.concatenate((newdataframeDispl, tempframeDispl), axis=0)
print(newdataframeForce.shape)
print(newdataframeDispl.shape)
# convert array into dataframe
DF = pd.DataFrame({'displ':newdataframeDispl[:,1],'force':newdataframeForce[:,1]})
# save the dataframe as a csv file
DF.to_csv("./TMresults/DisplvsForce.csv",sep=";",header=False,index=False)
plt.figure()
plt.title('Nodal displacement vs Nodal force (Y component) on top face in SMP', size=12)
plt.xlabel('Nodal displacement [m]', size=14)
plt.ylabel('Nodal force [N]', size=14)
plt.plot(newdataframeDispl[:,1],newdataframeForce[:,1])
plt.savefig('./TMresults/DisplvsForce.png', format='png')
plt.show()
......@@ -26,7 +26,7 @@ E_matrix = 771.e6 # (Pa)
epsilon_0 = 8.854 * 1.e-12 # vacuum electric permittivity
mag_mu_0 = 4.e-7 * math.pi # vacuum magnetic permeability
epsilon_r_smp = 1.0 # relative elec permittivity of SMP
mag_r_smp = 20.0 # relative mag permeability of SMP
mag_r_smp = 5.0 # relative mag permeability of SMP
# use of electric permittivity constants to formulate polarization
epsilon_4 = 0.0
epsilon_6 = 0.0
......@@ -76,8 +76,8 @@ Gvac=156.e1 # Shear modulus
Muvac = 0.35 # Poisson ratio
Evac= 2.0 * Gvac * (1.0 + Muvac) #youngs modulus
alphavac = betavac = gammavac = 0. # parameter of anisotropy
cpvac= 1012.*rhovac
Kxvac=Kyvac=Kzvac= 26.e-3 # thermal conductivity tensor components
cpvac= 1.e-12 #1012.*rhovac
Kxvac=Kyvac=Kzvac= 0.0 #26.e-3 # thermal conductivity tensor components
lxvac=lyvac=lzvac= 1.e-12 # electrical conductivity
seebeckvac= 0.
alphaThermvac= 0.0 # thermal dilatation
......
import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
Forcefilename = "force11111comp1.csv"
Displacementfilename = "NodalDisplacementPhysical11112Num1comp1.csv"
newdataframeForce = pd.read_csv(Forcefilename,sep=";",header=None).values
print(newdataframeForce.shape)
print()
newdataframeDispl = pd.read_csv(Displacementfilename,sep=";",header=None).values
print(newdataframeDispl.shape)
print()
# convert array into dataframe
DF = pd.DataFrame({'displ':newdataframeDispl[:,1],'force':newdataframeForce[:,1]})
# save the dataframe as a csv file
DF.to_csv("./DisplvsForce.csv",sep=";",header=False,index=False)
plt.figure()
plt.title('Nodal displacement vs Nodal force (Y component) on top face in SMP', size=12)
plt.xlabel('Nodal displacement [m]', size=14)
plt.ylabel('Nodal force [N]', size=14)
plt.plot(newdataframeDispl[:,1],newdataframeForce[:,1])
plt.savefig('./DisplvsForce.png', format='png')
plt.show()
# test file
set(PYFILE strongCoupling_SMPPhenoMech_withAir.py)
set(FILES2DELETE
*.msh
*.csv
)
add_cm3python_test(${PYFILE} "${FILES2DELETE}")
SetFactory("OpenCASCADE");
Mesh.Optimize = 1;
cm = 1e-2; // Unit
pp = "Input/10Geometric dimensions/0";
close_menu = 1;
colorpp = "Ivory";
// Smp
RintSMP = 0.95*cm;
RextSMP = 1*cm;
H = 2*cm;
DefineConstant[
RintSMP = {0.95*cm, Name StrCat[pp, "4Sample inner radius [m]"], Highlight Str[colorpp]}
RextSMP = {1*cm, Name StrCat[pp, "5Sample outer radius [m]"], Highlight Str[colorpp]}
H = {2*cm, Name StrCat[pp, "6Sample height [m]"], Highlight Str[colorpp]}
];
ECORE = 1000;
SKINECORE = 11110;
REFCORE = 11111;
REFCOREBOTTOM = 11115;
REFCOREPOINT = 11112;
REFCOREOUT = 11113;
REFCOREIN = 11114;
// Smp
vE()+=newv; Cylinder(newv) = {0,-H/2,0, 0,H,0, RintSMP};
vE()+=newv; Cylinder(newv) = {0,-H/2,0, 0,H,0, RextSMP};
vCoreE()+=newv; BooleanDifference(newv) = { Volume{vE(1)}; Delete; }{ Volume{vE(0)}; Delete; };
lc2 = RextSMP/4; // smp
Characteristic Length { PointsOf{ Volume{vCoreE()}; } } = lc2;
Recursive Color SteelBlue {Volume{vCoreE()};}
Physical Volume(ECORE) = vCoreE();
bnd_vCoreE() = CombinedBoundary{Volume{vCoreE()};};
Physical Surface(SKINECORE) = bnd_vCoreE();
Physical Surface(REFCORE)= bnd_vCoreE(1);
Physical Surface(REFCOREBOTTOM)= bnd_vCoreE(2);
Physical Surface(REFCOREOUT)= bnd_vCoreE(0);
Physical Surface(REFCOREIN)= bnd_vCoreE(3);
pts() = PointsOf {Surface { bnd_vCoreE(1) }; };
Physical Point(REFCOREPOINT)= pts(0);
Merge "undeformedSMP.msh";
//Include "deformed_mesh_filename.geo";
sl1 = newsl;
Surface Loop(sl1) = Surface{:};
//+
cm = 1e-2; // Unit
pp = "Input/10Geometric dimensions/0";
pp2 = "Input/10Geometric dimensions/1Shell radius/";
colorpp = "Ivory";
//Coil
//+
SetFactory("Built-in");
//+
p1 = newp;
Point(p1) = {0, 0.5, 0, 1.0};
//+
p2 = newp;
Point(p2) = {0.045, 0.5, 0, 1.0};
//+
p3 = newp;
Point(p3) = {-0.045, 0.5, 0, 1.0};
//+
p4 = newp;
Point(p4) = {0, 0.5, 0.045, 1.0};
//+
p5 = newp;
Point(p5) = {0, 0.5, -0.045, 1.0};
//+
c1 = newc;
Circle(c1) = {p5, p1, p3};
//+
c2 = newc;
Circle(c2) = {p3, p1, p4};
//+
c3 = newc;
Circle(c3) = {p4, p1, p2};
//+
c4 = newc;
Circle(c4) = {p2, p1, p5};
//+
p6 = newp;
Point(p6) = {0.075, 0.5, 0, 1.0};
//+
p7 = newp;
Point(p7) = {-0.075, 0.5, 0, 1.0};
//+
p8 = newp;
Point(p8) = {0, 0.5, 0.075, 1.0};
//+
p9 = newp;
Point(p9) = {0, 0.5, -0.075, 1.0};
//+
c5 = newc;
Circle(c5) = {p9, p1, p7};
//+
c6 = newc;
Circle(c6) = {p7, p1, p8};
//+
c7 = newc;
Circle(c7) = {p8, p1, p6};
//+
c8 = newc;
Circle(c8) = {p6, p1, p9};
//+
cl1 = newcl;
Curve Loop(cl1) = {c5,c6,c7,c8};
//+
cl2 = newcl;
Curve Loop(cl2) = {c1,c2,c3,c4};
//+
s1 = news;
Surface(s1) = {cl1, cl2};
//+
//#Transfinite Curve {c1,c2,c3,c4} = 11 Using Progression 1;
//+
//#Transfinite Curve {c5,c6,c7,c8} = 11 Using Progression 1;
//+
e() = Extrude {0, -1, 0} {
Surface{s1}; //#Layers {10};
};
/*
Printf("top curve = %g", e(0));
Printf("vol = %g", e(1));
Printf("surf = %g", e(2));
Printf("surf = %g", e(3));
Printf("surf = %g", e(4));
Printf("surf = %g", e(5));
Printf("surf = %g", e(6));
Printf("surf = %g", e(7));
Printf("surf = %g", e(8));
Printf("surf = %g", e(9));
*/
//#Transfinite Surface {e(2),e(3),e(4),e(5),e(6),e(7),e(8),e(9)};
sl2 = newsl;
Surface Loop(sl2) = Boundary{Volume{e(1)}; };
//+
Physical Surface(2222) = Boundary{Volume{e(1)}; };
Physical Volume(2000) = {e(1)};
DefineConstant[
Ly = {100*cm, Name StrCat[pp, "3Coil length along y-axis [m]"], Highlight Str[colorpp]}
wcoil = {3*cm, Name StrCat[pp, "4Coil width [m]"], Highlight Str[colorpp]}
RintCoil = {4.5*cm, Name StrCat[pp, "5Coil inner radius [m]"], Highlight Str[colorpp]}
RextCoil = {7.5*cm, Name StrCat[pp, "6Coil outer radius [m]"], Highlight Str[colorpp]}
];
//+
p10 = newp;
Point(p10) = {0, 0, 0, 1.0};
//+
//Air
//+
p11 = newp;
Point(p11) = {2.8, 0, 0, 1.0};
//+
p12 = newp;
Point(p12) = {-2.8, 0, 0, 1.0};
//+
p13 = newp;
Point(p13) = {0, 2.8, 0, 1.0};
//+
p14 = newp;
Point(p14) = {0, -2.8, 0, 1.0};
//+
c9 = newc;
Circle(c9) = {p11, p10, p13};
//+
c10 = newc;
Circle(c10) = {p13, p10, p12};
//+
c11 = newc;
Circle(c11) = {p12, p10, p14};
//+
c12 = newc;
Circle(c12) = {p14, p10, p11};
//+
//+
p15 = newp;
Point(p15) = {0, 0, 2.8, 1.0};
//+
p16 = newp;
Point(p16) = {0, 0, -2.8, 1.0};
//+
c13 = newc;
Circle(c13) = {p12, p10, p15};
//+
c14 = newc;
Circle(c14) = {p15, p10, p11};
//+
c15 = newc;
Circle(c15) = {p11, p10, p16};
//+
c16 = newc;
Circle(c16) = {p16, p10, p12};
//+
//+
c17 = newc;
Circle(c17) = {p13, p10, p15};
//+
c18 = newc;
Circle(c18) = {p15, p10, p14};
//+
c19 = newc;
Circle(c19) = {p14, p10, p16};
//+
c20 = newc;
Circle(c20) = {p16, p10, p13};
//+
//+
cl3 = newcl;
Curve Loop(cl3) = {c9,c17,c14};
//+
s2 = news;
Surface(s2) = {cl3};
//+
cl4 = newcl;
Curve Loop(cl4) = {c9,-c20,-c15};
//+
s3 = news;
Surface(s3) = {cl4};
//+
cl5 = newcl;
Curve Loop(cl5) = {c12,-c14,c18};
//+
s4 = news;
Surface(s4) = {cl5};
//+
cl6 = newcl;
Curve Loop(cl6) = {c12,c15,-c19};
//+
s5 = news;
Surface(s5) = {cl6};
//+
cl7 = newcl;
Curve Loop(cl7) = {c17,-c13,-c10};
//+
s6 = news;
Surface(s6) = {cl7};
//+
cl8 = newcl;
Curve Loop(cl8) = {c10,-c16,c20};
//+
s7 = news;
Surface(s7) = {cl8};
//+
cl9 = newcl;
Curve Loop(cl9) = {c13,c18,-c11};
//+
s8 = news;
Surface(s8) = {cl9};
//+
cl10 = newcl;
Curve Loop(cl10) = {-c11,-c16,-c19};
//+
s9 = news;
Surface(s9) = {cl10};
//+
sl3 = newsl;
Surface Loop(sl3) = {s2,s3,s4,s5,s6,s7,s8,s9};
v4 = newv;
Volume(v4) = {sl3,sl2,sl1};
Physical Surface(3333) = {s2,s3,s4,s5,s6,s7,s8,s9};
Physical Volume(3000) = v4;
DefineConstant[
Flag_Infinity = {0, Choices{0,1},
Name "Input/01Use shell transformation to infinity"}
];
// radius for surrounding air with transformation to infinity
If(Flag_Infinity==1)
label_Rext = "1Outer [m]";
Else
label_Rext = "1[m]";
EndIf
Rint = 200*cm;
Rext = 280*cm;
DefineConstant[
Rint = {200*cm, Min 0.15, Max 0.9, Step 0.1, Name StrCat[pp2, "7Inner [m]"],
Visible (Flag_Infinity==1), Highlight Str[colorpp] },
Rext = {280*cm, Min Rint, Max 1, Step 0.1, Name StrCat[pp2, StrCat["8", label_Rext]],
Visible 1, Highlight Str[colorpp] }
];
lc0 = Pi*Rint/6; // air
Characteristic Length { PointsOf{ Volume{v4}; } } = lc0;
lc1 = 2*wcoil/2; // coil
Characteristic Length { PointsOf{ Volume{e(1)}; } } = lc1;
//Mesh 2;
//Mesh 3;
//Mesh.SaveAll = 1;
//Save "mergedMesh.msh";
import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
filename = "/IPVolume1000val_EMFIELDSOURCEMean.csv"
f = "./"+filename
newdataframe = pd.read_csv(f,sep=";",header=None).values
print(newdataframe.shape)
print()
# convert array into dataframe
DF = pd.DataFrame(newdataframe)
# save the dataframe as a csv file
#DF.to_csv("./EMFIELDSOURCEMean.csv",sep=";",header=False,index=False)
plt.figure()
plt.title('EMFieldSource Mean vs Time', size=14)
plt.xlabel('Time [s]', size=14)
plt.ylabel('EMFieldSource Mean [W/m$^3$]', size=14)
plt.plot(newdataframe[:,0],newdataframe[:,1])
plt.savefig('./EMFieldSourceMean.png', format='png')
plt.show()
import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
filename = "/IPVolume1000val_TEMPERATUREMean.csv"
f = "./"+filename
newdataframe = pd.read_csv(f,sep=";",header=None).values
print(newdataframe.shape)
print()
# convert array into dataframe
DF = pd.DataFrame(newdataframe)
# save the dataframe as a csv file
#DF.to_csv("./TEMPERATUREMean.csv",sep=";",header=False,index=False)
plt.figure()
plt.title('Temperature Mean vs Time', size=14)
plt.xlabel('Time [s]', size=14)
plt.ylabel('Temperature Mean [K]', size=14)
plt.plot(newdataframe[:,0],newdataframe[:,1])
plt.savefig('./TemperatureMean.png', format='png')
plt.show()
import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
filename = "/IPVolume1000val_VOLTAGEMean.csv"
f = "./"+filename
newdataframe = pd.read_csv(f,sep=";",header=None).values
print(newdataframe.shape)
print()
# convert array into dataframe
DF = pd.DataFrame(newdataframe)
# save the dataframe as a csv file
#DF.to_csv("./VOLTAGEMean.csv",sep=";",header=False,index=False)
plt.figure()
plt.title('Voltage Mean vs Time', size=14)
plt.xlabel('Time [s]', size=14)
plt.ylabel('Voltage Mean [V]', size=14)
plt.plot(newdataframe[:,0],newdataframe[:,1])
plt.savefig('./VoltageMean.png', format='png')
plt.show()
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment