Merge branch 'develop' into feature/fluidized_bed_chemistry_cases

tutorials/cfdemSolverRhoPimpleChem/PolydisperseFluidizedBed/R2_FB/DEM/in.liggghts_run
This commit is contained in:
danielque
2023-02-03 12:37:51 +01:00
75 changed files with 1044 additions and 21134 deletions

View File

@ -1,33 +1,39 @@
clear all;
clc;
%% Script to plot the angular velocities and the particle positions
% plot the angular velocity and the particle position
A = importdata('../../DEM/post/angular_velocity_no_coupling.txt',' ',1);
B = importdata('../../DEM/post/position_no_coupling.txt',' ',1);
C = importdata('../../DEM/post/angular_velocity.txt',' ',1);
D = importdata('../../DEM/post/position.txt',' ',1);
A = importdata('../../DEM/post/angular_velocity.txt',' ',1);
B = importdata('../../DEM/post/position.txt',' ',1);
pos1 = B.data();
omega1 = A.data();
pos2 = D.data();
omega2 = C.data();
omega = A.data();
pos = B.data();
time = omega1(:,1);
omegax1 = omega1(:,2);
omegay1 = omega1(:,3);
omegaz1 = omega1(:,4);
time1 = omega(:,1);
omegax = omega(:,2);
omegay = omega(:,3);
omegaz = omega(:,4);
time2 = omega2(:,1);
omegax2 = omega2(:,2);
omegay2 = omega2(:,3);
omegaz2 = omega2(:,4);
time2 = pos(:,1);
posx = pos(:,2);
posy = pos(:,3);
posz = pos(:,4);
figure
plot(time,omegaz1,'-.-',time2,omegaz2,'Linewidth',1.5)
figure(1)
plot(time1,omegaz,'-.-','Linewidth',1.5)
xlabel('Time (s)')
ylabel('Angular velocity (1/s)')
legend('One-way coupling','Two-way coupling')
legend('Two-way torque coupling')
axis([0 0.5 0 11])
set(gca,'FontSize',12)
print('angular_velocity_compare.eps')
print('angular_velocity.eps')
figure(2)
plot(time2,posy,'-.-','Linewidth',1.5)
xlabel('Time (s)')
ylabel('Y-position (m)')
legend('Two-way torque coupling')
axis([0 0.5 0 0.1])
set(gca,'FontSize',12)
print('position.eps')

View File

@ -0,0 +1,21 @@
{
"type" : "CFDEMcoupling",
"runs" : [
{
"name" : "cfdemrun",
"solver" : "cfdemSolverIB",
"type" : "CFDEMcoupling/mpi",
"nprocs" : 4,
"data" : {
"series" : [
{"name" : "angular_velocity", "file" : "angular_velocity.txt", "columns" : ["t", "omegax", "omegay", "omegaz"]},
{"name" : "position", "file" : "position.txt", "columns" : ["t", "x", "y", "z"]}
],
"plots" : [
{"name" : "z-omega", "title" : "Particle z-Angular Velocity", "xdata" : "angular_velocity.t", "ydata" : ["angular_velocity.omegaz"], "xlabel" : "time [s]", "ylabel" : "angular velocity [1/s]"},
{"name" : "y-position", "title" : "Particle y-Position", "xdata" : "position.t", "ydata" : ["position.y"], "xlabel" : "time [s]", "ylabel" : "position [m]"}
]
}
}
]
}

View File

@ -45,7 +45,7 @@ plt.xlabel('Time (s)')
plt.ylabel('Min. eigen value of G')
plt.minorticks_on()
plt.grid()
plt.savefig('RBC_Posieuille_Gyration.eps')
plt.savefig('RBC_Poiseuille_Gyration.eps')
plt.show()

View File

@ -0,0 +1,27 @@
{
"type" : "CFDEMcoupling",
"runs" : [
{
"name" : "liggghts-init",
"input_script" : "DEM/in.liggghts_init",
"type" : "liggghts/serial"
},
{
"name" : "cfdemrun",
"depends_on" : "liggghts-init",
"solver" : "cfdemSolverIBContinuousForcing",
"type" : "CFDEMcoupling/mpi",
"nprocs" : 2,
"data" : {
"series" : [
{"name" : "position", "file" : "particle_position.txt", "columns" : ["t", "x1", "y1", "z1", "x2", "y2", "z2", "x3", "y3", "z3", "x4", "y4", "z4", "x5", "y5", "z5", "x6", "y6", "z6", "x7", "y7", "z7", "x8", "y8", "z8", "x9", "y9", "z9", "x10", "y10", "z10"]}
],
"plots" : [
{"name" : "x-position", "title" : "Particle x-Position", "xdata" : "position.t", "ydata" : ["position.x1","position.x2","position.x3","position.x4","position.x5","position.x6","position.x7","position.x8","position.x9","position.x10"], "xlabel" : "time [1e-6 s]", "ylabel" : "position [1e-6 m]"},
{"name" : "y-position", "title" : "Particle y-Position", "xdata" : "position.t", "ydata" : ["position.y1","position.y2","position.y3","position.y4","position.y5","position.y6","position.y7","position.y8","position.y9","position.y10"], "xlabel" : "time [1e-6 s]", "ylabel" : "position [1e-6 m]"},
{"name" : "z-position", "title" : "Particle z-Position", "xdata" : "position.t", "ydata" : ["position.z1","position.z2","position.z3","position.z4","position.z5","position.z6","position.z7","position.z8","position.z9","position.z10"], "xlabel" : "time [1e-6 s]", "ylabel" : "position [1e-6 m]"}
]
}
}
]
}

View File

@ -0,0 +1,27 @@
{
"type" : "CFDEMcoupling",
"runs" : [
{
"name" : "liggghts-init",
"input_script" : "DEM/in.liggghts_init",
"type" : "liggghts/serial"
},
{
"name" : "cfdemrun",
"depends_on" : "liggghts-init",
"solver" : "cfdemSolverIBContinuousForcing",
"type" : "CFDEMcoupling/mpi",
"nprocs" : 2,
"data" : {
"series" : [
{"name" : "position", "file" : "particle_position.txt", "columns" : ["t", "x1", "y1", "z1", "x2", "y2", "z2", "x3", "y3", "z3", "x4", "y4", "z4", "x5", "y5", "z5", "x6", "y6", "z6", "x7", "y7", "z7", "x8", "y8", "z8", "x9", "y9", "z9", "x10", "y10", "z10"]}
],
"plots" : [
{"name" : "x-position", "title" : "Particle x-Position", "xdata" : "position.t", "ydata" : ["position.x1","position.x2","position.x3","position.x4","position.x5","position.x6","position.x7","position.x8","position.x9","position.x10"], "xlabel" : "time [1e-6 s]", "ylabel" : "position [1e-6 m]"},
{"name" : "y-position", "title" : "Particle y-Position", "xdata" : "position.t", "ydata" : ["position.y1","position.y2","position.y3","position.y4","position.y5","position.y6","position.y7","position.y8","position.y9","position.y10"], "xlabel" : "time [1e-6 s]", "ylabel" : "position [1e-6 m]"},
{"name" : "z-position", "title" : "Particle z-Position", "xdata" : "position.t", "ydata" : ["position.z1","position.z2","position.z3","position.z4","position.z5","position.z6","position.z7","position.z8","position.z9","position.z10"], "xlabel" : "time [1e-6 s]", "ylabel" : "position [1e-6 m]"}
]
}
}
]
}

View File

@ -173,6 +173,9 @@ FinesFieldsProps
poresizeWidth 0.2;
verbose true;
diffCoeff 0.005;
tauRelease 0.1;
depositionLength 0.1;
kineticClogging true;
}
gravityProps

View File

@ -23,13 +23,13 @@ startTime 0;
stopAt endTime;
endTime 40;
endTime 10;
deltaT 0.0005;
writeControl timeStep;
writeInterval 8000;
writeInterval 1000;
purgeWrite 0;

View File

@ -2,3 +2,4 @@ cd CFD
blockMesh
topoSet -dict system/topoSetDict
cp -r orig.0/ 0
decomposePar

View File

@ -77,11 +77,11 @@ fix cfd5 ore chem/shrink/core speciesA CO molMassA 0.02801 speciesC CO2 molM
fix cfd6 ore chem/shrink/core speciesA H2 molMassA 0.00202 speciesC H2O molMassC 0.01801 scale_reduction_rate ${rate_scale}
# material properties for chemical reaction
fix k0_CO ore property/atom k0_cfd5 vector yes no no 17 25 2700
fix Ea_CO ore property/atom Ea_cfd5 vector yes no no 69488 95000 130940
fix k0_CO ore property/global k0_cfd5 vector 17 25 2700
fix Ea_CO ore property/global Ea_cfd5 vector 69488 95000 130940
fix k0_H2 ore property/atom k0_cfd6 vector yes no no 30 23 160
fix Ea_H2 ore property/atom Ea_cfd6 vector yes no no 63627 85000 105908
fix k0_H2 ore property/global k0_cfd6 vector 30 23 160
fix Ea_H2 ore property/global Ea_cfd6 vector 63627 85000 105908
# particle parameters
fix porosity ore property/global porosity_ore vector 0.61 0.34 0.19 0.17

View File

@ -62,11 +62,11 @@ fix cfd5 all chem/shrink/core speciesA CO molMassA 0.02801 speciesC CO2 molM
fix cfd6 all chem/shrink/core speciesA H2 molMassA 0.00202 speciesC H2O molMassC 0.01801 scale_reduction_rate 10.0 screen yes
# Chemical properties for unreacted shrink core (activate only when chem/shrink/core is active)
fix k0_CO all property/atom k0_cfd5 vector yes no no 17 25 2700
fix Ea_CO all property/atom Ea_cfd5 vector yes no no 69488 73674 113859
fix k0_CO all property/global k0_cfd5 vector 17 25 2700
fix Ea_CO all property/global Ea_cfd5 vector 69488 73674 113859
fix k0_H2 all property/atom k0_cfd6 vector yes no no 30 23 160
fix Ea_H2 all property/atom Ea_cfd6 vector yes no no 63627 71162 92092
fix k0_H2 all property/global k0_cfd6 vector 30 23 160
fix Ea_H2 all property/global Ea_cfd6 vector 63627 71162 92092
# particle porosity/tortuosity/pore diameter
fix porosity all property/global porosity_all vector 0.65 0.31 0.16 0.15
@ -95,20 +95,6 @@ variable ke_tot equal c_KinEn
# print total kinetic energy
fix printCompute all print ${WI} "${time} ${ke_tot}" file ../DEM/post/printKE.txt title "#time ke_tot"
compute Ea_CO all reduce sum f_Ea_CO[1] f_Ea_CO[2] f_Ea_CO[3]
fix Ea all ave/time 1 1 1 c_Ea_CO[1] c_Ea_CO[2] c_Ea_CO[3]
variable Ea1 equal f_Ea[1]
variable Ea2 equal f_Ea[2]
variable Ea3 equal f_Ea[3]
compute k0CO all reduce sum f_k0_CO[1] f_k0_CO[2] f_k0_CO[3]
fix k0 all ave/time 1 1 1 c_k0CO[1] c_k0CO[2] c_k0CO[3]
variable k01 equal f_k0[1]
variable k02 equal f_k0[2]
variable k03 equal f_k0[3]
fix printk0Ea all print ${WI} "${time} ${Ea1} ${Ea2} ${Ea3} ${k01} ${k02} ${k03}" file ../DEM/post/k0Ea.dat title "#time Ea1 Ea2 Ea3 k01 k02 k03"
###############
# Print out values affecting chemical reduction into specified folder for given time
# Diffusion Coefficient for CO and H2

View File

@ -62,11 +62,11 @@ fix cfd5 all chem/shrink/core speciesA CO molMassA 0.02801 speciesC CO2 molM
fix cfd6 all chem/shrink/core speciesA H2 molMassA 0.00202 speciesC H2O molMassC 0.01801 scale_reduction_rate 10.0 screen yes
# Chemical properties for unreacted shrink core (activate only when chem/shrink/core is active)
fix k0_CO all property/atom k0_cfd5 vector yes no no 17 25 2700
fix Ea_CO all property/atom Ea_cfd5 vector yes no no 69488 73674 113859
fix k0_CO all property/global k0_cfd5 vector 17 25 2700
fix Ea_CO all property/global Ea_cfd5 vector 69488 73674 113859
fix k0_H2 all property/atom k0_cfd6 vector yes no no 30 23 160
fix Ea_H2 all property/atom Ea_cfd6 vector yes no no 63627 71162 92092
fix k0_H2 all property/global k0_cfd6 vector 30 23 160
fix Ea_H2 all property/global Ea_cfd6 vector 63627 71162 92092
# particle porosity/tortuosity/pore diameter
fix porosity all property/global porosity_all vector 0.65 0.31 0.16 0.15
@ -95,20 +95,6 @@ variable ke_tot equal c_KinEn
# print total kinetic energy
fix printCompute all print ${WI} "${time} ${ke_tot}" file ../DEM/post/printKE.txt title "#time ke_tot"
compute Ea_CO all reduce sum f_Ea_CO[1] f_Ea_CO[2] f_Ea_CO[3]
fix Ea all ave/time 1 1 1 c_Ea_CO[1] c_Ea_CO[2] c_Ea_CO[3]
variable Ea1 equal f_Ea[1]
variable Ea2 equal f_Ea[2]
variable Ea3 equal f_Ea[3]
compute k0CO all reduce sum f_k0_CO[1] f_k0_CO[2] f_k0_CO[3]
fix k0 all ave/time 1 1 1 c_k0CO[1] c_k0CO[2] c_k0CO[3]
variable k01 equal f_k0[1]
variable k02 equal f_k0[2]
variable k03 equal f_k0[3]
fix printk0Ea all print ${WI} "${time} ${Ea1} ${Ea2} ${Ea3} ${k01} ${k02} ${k03}" file ../DEM/post/k0Ea.dat title "#time Ea1 Ea2 Ea3 k01 k02 k03"
###############
# Print out values affecting chemical reduction into specified folder for given time
# Diffusion Coefficient for CO and H2

View File

@ -61,8 +61,8 @@ fix cfd3 all couple/cfd/chemistry n_species 3 species_names CO CO2 N2 n_diff
fix cfd5 all chem/shrink/core speciesA CO molMassA 0.02801 speciesC CO2 molMassC 0.04401 scale_reduction_rate 10.0 screen yes
# Chemical properties for unreacted shrink core (activate only when chem/shrink/core is active)
fix k0_CO all property/atom k0_cfd5 vector yes no no 17 25 2700
fix Ea_CO all property/atom Ea_cfd5 vector yes no no 69488 73674 113859
fix k0_CO all property/global k0_cfd5 vector 17 25 2700
fix Ea_CO all property/global Ea_cfd5 vector 69488 73674 113859
# particle porosity/tortuosity/pore diameter
fix porosity all property/global porosity_all vector 0.65 0.31 0.16 0.15
@ -91,20 +91,6 @@ variable ke_tot equal c_KinEn
# print total kinetic energy
fix printCompute all print ${WI} "${time} ${ke_tot}" file ../DEM/post/printKE.txt title "#time ke_tot"
compute Ea_CO all reduce sum f_Ea_CO[1] f_Ea_CO[2] f_Ea_CO[3]
fix Ea all ave/time 1 1 1 c_Ea_CO[1] c_Ea_CO[2] c_Ea_CO[3]
variable Ea1 equal f_Ea[1]
variable Ea2 equal f_Ea[2]
variable Ea3 equal f_Ea[3]
compute k0CO all reduce sum f_k0_CO[1] f_k0_CO[2] f_k0_CO[3]
fix k0 all ave/time 1 1 1 c_k0CO[1] c_k0CO[2] c_k0CO[3]
variable k01 equal f_k0[1]
variable k02 equal f_k0[2]
variable k03 equal f_k0[3]
fix printk0Ea all print ${WI} "${time} ${Ea1} ${Ea2} ${Ea3} ${k01} ${k02} ${k03}" file ../DEM/post/k0Ea.dat title "#time Ea1 Ea2 Ea3 k01 k02 k03"
###############
# Print out values affecting chemical reduction into specified folder for given time
# Diffusion Coefficient for CO and H2

View File

@ -61,8 +61,8 @@ fix cfd3 all couple/cfd/chemistry n_species 3 species_names CO CO2 N2 n_diff
fix cfd5 all chem/shrink/core speciesA CO molMassA 0.02801 speciesC CO2 molMassC 0.04401 scale_reduction_rate 10.0 screen yes
# Chemical properties for unreacted shrink core (activate only when chem/shrink/core is active)
fix k0_CO all property/atom k0_cfd5 vector yes no no 17 25 2700
fix Ea_CO all property/atom Ea_cfd5 vector yes no no 69488 73674 113859
fix k0_CO all property/global k0_cfd5 vector 17 25 2700
fix Ea_CO all property/global Ea_cfd5 vector 69488 73674 113859
# particle porosity/tortuosity/pore diameter
fix porosity all property/global porosity_all vector 0.65 0.31 0.16 0.15
@ -91,20 +91,6 @@ variable ke_tot equal c_KinEn
# print total kinetic energy
fix printCompute all print ${WI} "${time} ${ke_tot}" file ../DEM/post/printKE.txt title "#time ke_tot"
compute Ea_CO all reduce sum f_Ea_CO[1] f_Ea_CO[2] f_Ea_CO[3]
fix Ea all ave/time 1 1 1 c_Ea_CO[1] c_Ea_CO[2] c_Ea_CO[3]
variable Ea1 equal f_Ea[1]
variable Ea2 equal f_Ea[2]
variable Ea3 equal f_Ea[3]
compute k0CO all reduce sum f_k0_CO[1] f_k0_CO[2] f_k0_CO[3]
fix k0 all ave/time 1 1 1 c_k0CO[1] c_k0CO[2] c_k0CO[3]
variable k01 equal f_k0[1]
variable k02 equal f_k0[2]
variable k03 equal f_k0[3]
fix printk0Ea all print ${WI} "${time} ${Ea1} ${Ea2} ${Ea3} ${k01} ${k02} ${k03}" file ../DEM/post/k0Ea.dat title "#time Ea1 Ea2 Ea3 k01 k02 k03"
###############
# Print out values affecting chemical reduction into specified folder for given time
# Diffusion Coefficient for CO and H2

View File

@ -31,6 +31,14 @@ dumpIndexDisplacementIncrement 8000;
dumpIndexInputIncrement 4000;
posIndex 0;
posX 1;
posY 2;
posZ 3;
nNeighMin 3;
timePerDisplacementStep 0.1;

View File

@ -52,7 +52,7 @@ surfaceScalarFields
//verbose true;
couplingSubStep 3;
couplingSubStep 7;
initialRecSteps 150;

View File

@ -1,49 +0,0 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.0 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
thermoType
{
type hePsiThermo;
mixture pureMixture;
transport const;
thermo eConst;
equationOfState perfectGas;
specie specie;
energy sensibleInternalEnergy;
}
mixture
{
specie
{
nMoles 1;
molWeight 42.1;
}
thermodynamics
{
Cv 1310;
Hf 0;
}
transport
{
mu 1e-05;
Pr 0.72;
}
}
// ************************************************************************* //

View File

@ -27,4 +27,8 @@ Cv Cv [ 0 2 -2 -1 0 0 0 ] 1310;
molM molM [1 0 0 0 -1 0 0 ] 0.0421;
TMax TMax [0 0 0 1 0 0 0] 400;
TMin TMin [0 0 0 1 0 0 0] 250;
// ************************************************************************* //

File diff suppressed because it is too large Load Diff

View File

@ -51,6 +51,63 @@ libs (
functions
{
Tave
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
log true;
writeControl timeStep;
writeInterval 10;
writeFields false;
regionType all;
operation volAverage;
weightField voidfractionRec;
fields
(
T
);
}
fluxprobe0
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 100;
log true;
writeFields false;
regionType faceZone;
name fzcenter;
operation sum;
fields
(
phi
);
}
fluxprobe1
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 100;
log true;
writeFields false;
regionType faceZone;
name fz0;
operation sum;
fields
(
phi
);
}
fieldAverage1
{
type fieldAverage;
@ -125,13 +182,12 @@ functions
);
}
fieldOutput
{
type writeObjects;
functionObjectLibs ( "libutilityFunctionObjects.so" );
exclusiveWriting true;
objects ("rhoRecMean" "voidfraction" "voidfractionRec" "voidfractionMean" "voidfractionRecMean" "addSourceMean" "phiSMean" "TMean" "T" "particleTemp" "particleTempMean");
objects ("rhoRecMean" "voidfraction" "voidfractionRec" "voidfractionMean" "voidfractionRecMean" "addSourceMean" "TMean" "T" "partTemp" "partTempMean");
writeControl timeStep;
writeInterval 40000;
}
@ -177,10 +233,10 @@ functions
scalarList newWeights(2);
newWeights[0] = w0;
newWeights[1] = 1-w0;
/*
scalarList newWeights(1);
newWeights[0] = 1.0;
*/
// scalarList newWeights(1);
// newWeights[0] = 1.0;
weightDict.set("weights",newWeights);
#};
}

View File

@ -28,13 +28,6 @@ gradSchemes
divSchemes
{
default Gauss linear;
div(phi,U) Gauss limitedLinearV 1;
div(phi,k) Gauss limitedLinear 1;
div(phi,epsilon) Gauss limitedLinear 1;
div(phi,R) Gauss limitedLinear 1;
div(R) Gauss linear;
div(phi,nuTilda) Gauss limitedLinear 1;
div((nuEff*dev(T(grad(U))))) Gauss linear;
div(phi,T) Gauss limitedLinear 1;
// div(phi,T) Gauss upwind;
}
@ -57,7 +50,6 @@ snGradSchemes
fluxRequired
{
default no;
p ;
T ;
}

View File

@ -17,14 +17,6 @@ FoamFile
solvers
{
"(p|rho)"
{
solver PCG;
preconditioner DIC;
tolerance 1e-6;
relTol 0.01;
}
"(correctedField)"
{
solver PCG;
@ -33,26 +25,6 @@ solvers
relTol 0.05;
}
"(p|rho)Final"
{
$p;
relTol 0;
}
"(U|k|e|epsilon|R|nuTilda|c)"
{
solver smoothSolver;
smoother GaussSeidel;
tolerance 1e-05;
relTol 0;
}
"(U|e|k|nuTilda)Final"
{
$U;
relTol 0;
}
T
{
solver PBiCG;

View File

@ -39,7 +39,6 @@ variable rfPW2 equal 0.1
variable Tpart equal 330
variable dt equal 0.0025
#0.001
variable skin equal 0.0005
###############################################
@ -112,13 +111,9 @@ variable Tave equal v_sumT1/v_np
fix printheat all print 10 "${time} ${Tave}" file ../DEM/temp_ave.txt title "#time T_ave"
#fix tdist1 all ave/histo 1000 1 1000 324 380 56 f_Temp[0] mode vector file ../DEM/temp_histo.txt title1 "Particle temperatures [K]"
#fix tdist2 all ave/histo 1000 1 1000 324 380 448 f_Temp[0] mode vector file ../DEM/temp_histo_fine.txt title1 "Particle temperatures [K]"
fix tdist1 all ave/histo 400 1 400 324 380 56 f_Temp[0] mode vector file ../DEM/temp_histo.txt title1 "Particle temperatures [K]"
fix tdist2 all ave/histo 400 1 400 324 380 448 f_Temp[0] mode vector file ../DEM/temp_histo_fine.txt title1 "Particle temperatures [K]"
fix integr all nve/sphere
#screen output
@ -131,3 +126,6 @@ dump dmp all custom/vtk 10000 ../DEM/post/dump*.liggghts_coupled.vtk id type x
run 1
set region total property/atom Temp ${Tpart}
# heat sources are already set in the restart file; if they would not have been - or if restarting in a different way -, set them now
#set region insReg property/atom heatSource 3.3379e-04 # = V_p x \dot{q}