guojuw
Forum Replies Created
-
AuthorPosts
-
February 28, 2021 at 1:50 am in reply to: SuperLatticePhysWallShearStress3D() returns norm of wall shear stress tensor #5498guojuwParticipant
Dear Jonathan,
Thanks for your idea. Here is my solution.
1. create SuperVTMwriter3D<T> rawWriter(“rawvtk”,false,false);
2. dump uncompressed data in text, vtkWriter.addFunctor(wallShearStressLidVec);
rawWriter.write(iT);3. modify rawWriter.write(iT) so that it writes data in the region of interest only.
4. process the vtk files directly using python scripts, since they contain all information I need.
I am trying your approach as well.
Regards,
Junwei GuoFebruary 25, 2021 at 4:19 am in reply to: SuperLatticePhysWallShearStress3D() returns norm of wall shear stress tensor #5493guojuwParticipantDear Jonathan,
I have another related problem.
After I have calculated the WSS using SuperLatticePhysWallShearStress3D(). I want to write WSS to disk as a CSV text file for post-processing purposes. I tried two different ways:
1.
a. calculate the WSS using ‘SuperLatticePhysWallShearStress3D<T,DESCRIPTOR> wallShearStressBot(sLattice,superGeometry,4,converter,base)’
b. write data to vtk files, ‘SuperVTMwriter3D<T> vtkWriter(“ellipsoid”);’ then ‘vtkWriter.addFunctor(wallShearStressBotVec);’ and last ‘vtkWriter.write(iT);’
c. extra the WSS using paraview and write WSS as CSV file.This way is very inconvenient.
2.
a. same as 1.a
b. use the interpolate utility AnalyticalFfromSuperF3D<T> intpolateBotWSS(wallShearStressBot,true); then interpolate the WSS at a given location intpolateBotWSS(WSSBot,loc);
c. write WSSBot and loc to disk.This way seems better. However, WSS is calculated at a given 2D plane (usually curved) in SuperLatticePhysWallShearStress3D, AnalyticalFfromSuperF3D<T> interpolates the results incorrectly in 3D (the value of WSS outside the plane is 0).
Do you recommend any better ways to write WSS as CSV text to disk?
Regards,
Junwei GuoFebruary 25, 2021 at 2:20 am in reply to: SuperLatticePhysWallShearStress3D() returns norm of wall shear stress tensor #5492guojuwParticipantThanks Jonathan. I worked it out by outputting the WSS[0] to WSS[2].
Regards,
JunweiguojuwParticipantDear Stephan,
I worked it out in a different but similar way:
1. Add “#include ..” in the source code
2. In $openLBRoot/config.mk, add this line “LDFLAGS :=gsl-config --cflags --libs
3. gsl-config is compiled from GSL source code.Thanks for your help.
Regards,
Junwei GUOOctober 18, 2020 at 11:27 pm in reply to: HLBM: particle-particle interaction and fluid-particle momentum exchange #5234guojuwParticipantDear Nicolas,
I am looking forward to attending the next OpenLB spring school. That would happen next year. Can you please contact the developers to solve this problem?
Thanks again.
Regards,
Junwei GUoguojuwParticipantDear Nicolas,
Thanks again. Looking forward to seeing your new releases. I worked out the SmoothIndicatorEllipsoid3D these few days (need to adjust the epsilon at HLBM object bounds). I will compare my code with your official release.
Regards,
Junwei Guo//Constructor: SmoothIndicatorEllipsoid3D
template <typename T, typename S, bool HLBM>
SmoothIndicatorEllipsoid3D<T,S,HLBM>::SmoothIndicatorEllipsoid3D(Vector<S,3> center, S xLength, S yLength, S zLength, S epsilon, Vector<S,3> theta, S density, Vector<S,3> vel)
: _xLength(xLength),_yLength(yLength),_zLength(zLength)
{
this->_pos = center;
T xr = xLength*.5;
T yr = yLength*.5;
T zr = zLength*.5;
T maxr = max({xr,yr,zr});
//this->_circumRadius = .5*(sqrt(xLength*xLength+yLength*yLength+zLength*zLength))+0.5*epsilon;
this->_circumRadius = maxr + 0.5*epsilon;
this->_name = “1000”; //Ellipsoid particle code;
this->_myMin = {
center[0] – this->getCircumRadius(),
center[1] – this->getCircumRadius(),
center[2] – this->getCircumRadius()
};
this->_myMax = {
center[0] + this->getCircumRadius(),
center[1] + this->getCircumRadius(),
center[2] + this->getCircumRadius()
};
this->_epsilon = epsilon;
this->_theta = {
theta[0] * (M_PI/180.),
theta[1] * (M_PI/180.),
theta[2] * (M_PI/180.)
};
//T mass = xLength*yLength*zLength*density;
T mass = 4./3.*M_PI*xLength*yLength*zLength*density/8.;T xLength2 = xLength*xLength;
T yLength2 = yLength*yLength;
T zLength2 = zLength*zLength;
//https://en.wikipedia.org/wiki/List_of_moments_of_inertia
Vector<S,3> mofi;
mofi[0] = 1./5.*mass*(yLength2+zLength2)/4.;
mofi[1] = 1./5.*mass*(xLength2+zLength2)/4.;
mofi[2] = 1./5.*mass*(yLength2+xLength2)/4.;
this->init(this->_theta, vel, mass, mofi);
}template <typename T, typename S, bool HLBM>
bool SmoothIndicatorEllipsoid3D<T,S,HLBM>::operator()(T output[], const S input[])
{
//1.Calculate distance between point and center of unrotated indicator
//cout << “SmoothIndicatorEllipsoid3D<T,S,HLBM>::operator(). Input=(” <<input[0]<<“,”<<input[1]<<“,”<< input[2]<<“). Output=(“<<output[0]<<“).”<< endl;
//cout << ” ” <<input[0]<<“, “<<input[1]<<“, “<< input[2]<<“, “;
T xDist = input[0] – this->getPos()[0];
T yDist = input[1] – this->getPos()[1];
T zDist = input[2] – this->getPos()[2];//2.Calculate point projected to rotated indicator
// counter-clockwise rotation by _theta=-theta around center// rotating the object BACK to normal??
// T x = this->getPos()[0] + this->getRotationMatrix()[0]*xDist + this->getRotationMatrix()[3]*yDist + this->getRotationMatrix()[6]*zDist;
// T y = this->getPos()[1] + this->getRotationMatrix()[1]*xDist + this->getRotationMatrix()[4]*yDist + this->getRotationMatrix()[7]*zDist;
// T z = this->getPos()[2] + this->getRotationMatrix()[2]*xDist + this->getRotationMatrix()[5]*yDist + this->getRotationMatrix()[8]*zDist;T x = this->getPos()[0]*.0 + this->getRotationMatrix()[0]*xDist + this->getRotationMatrix()[3]*yDist + this->getRotationMatrix()[6]*zDist;
T y = this->getPos()[1]*.0 + this->getRotationMatrix()[1]*xDist + this->getRotationMatrix()[4]*yDist + this->getRotationMatrix()[7]*zDist;
T z = this->getPos()[2]*.0 + this->getRotationMatrix()[2]*xDist + this->getRotationMatrix()[5]*yDist + this->getRotationMatrix()[8]*zDist;T xp = input[0];
T yp = input[1];
T zp = input[2];
T xc = this->getPos()[0];
T yc = this->getPos()[1];
T zc = this->getPos()[2];
T xr = _xLength/2.;
T yr = _yLength/2.;
T zr = _zLength/2.;
T unrotatedInside = ((xp-xc)/xr)*((xp-xc)/xr) + ((yp-yc)/yr)*((yp-yc)/yr) +((zp-zc)/zr)*((zp-zc)/zr) -1.0;T p1 = x/xr;
T p2 = y/yr;
T p3 = z/zr;
T rotatedInside = p1*p1+p2*p2+p3*p3-1.0;if (rotatedInside<.0)
{
output[0]=1.0;
return true;
}
output[0]=.0;
return false;
}guojuwParticipantDear Dr. Mathias,
I think the current hard sphere collision model will result in spheres overlapping if the suspension is not dilute or multiple contacts occur. Would you mind to share your code used for elastic collision (the implementation of parameter gamma) showing in your paper?
Thanks,
Junwei GuoguojuwParticipantDear Mathias,
I mean the parameter gamma in equation(9). How to modify this particle collision parameter?
Regards,
Junwei GuoArticle (Krause2017)
Krause, M. J.; Klemens, F.; Henn, T.; Trunk, R. & Nirschl, H.
Particle flow simulations with homogenised lattice Boltzmann methods
Particuology, Elsevier, 2017 , 34 , 1-13guojuwParticipantDear Nicolas,
I just know the version 1.3r1 is released. Actually I fixed some other small bugs in version 1.3r0.
Do you know how 1.3r1 is changed from 1.3r0 in detail? So I can migrate my code to the newest version?
Regards,
Junwei GuoguojuwParticipantDear Dr. Mathias,
I think the channel width in simulation is correctly set. In spite of the convergence problem, there are still a few other forcing problems.
I used these commands to output the forces:
SuperLatticeGuoZhaoPhysBodyForce2D<T, DESCRIPTOR> bodyForce(sLattice, converter);
SuperLatticePhysDarcyForce2D<T, DESCRIPTOR> darcyForce(sLattice, superGeometry,1 ,converter);
SuperLatticeField2D<T,DESCRIPTOR,olb::descriptors::FORCE> FORCE_field (sLattice);The body force outputted is 10, however the imposed body force is 10000. The output Darcy force is always zero. The external force “olb::descriptors::FORCE>” is also not matching.
Regards,
Junwei GuoguojuwParticipantDear Dr. Mathias,
I tried the resolutions from 20 to 200. In each case, the ratio of dx^2/dt is consistent, which means the kinematic viscosity conversion factors are unique in all cases.
The velocity profiles did not converge to the analytical one, but still converge to the highest resolution case.
I plotted the velocity profiles for all cases here. The L2 error compared with the analytical solution is alose compared.
Regards,
Junwei GuoAugust 23, 2019 at 7:48 pm in reply to: Reference for strain rate tensor prefactor in powerLawBGKdynamics #4526guojuwParticipantDear Davide,
Great thanks for your detailed explanation.
I think the paper by Chai et al. (2011) is more clear. The pi in openlb stands for the second order moment of nonequilibrium distribution (the summation part of eq.(18) in Chai et al. (2011)). And your prefactor ‘pre2’ is the square of the prefactor before the summation part of eq.(18). Other formulation in your codes are consistent with eq.(1-4).Thanks again.
Regards,
Junwei GuoChai2011 – Multiple Relaxation Time Lattice Boltzmann Model for Generalized Newtonian Fluid Flows.pdf
August 22, 2019 at 11:22 pm in reply to: Printing the relaxation time of each lattice to vtk files #4523guojuwParticipantDear Aurelio,
Great! Thanks a lot. It helps me.
Regard,
Junwei GuoAugust 21, 2019 at 8:48 pm in reply to: Printing the relaxation time of each lattice to vtk files #4516guojuwParticipantHi Nicolas,
I mean printing the non-constant relaxation time from non-Newtonian modeling, for example, the power-law case in $openlbRoot/examples/laminar/powerLaw2d.
In other words, I just want to see the relaxation time from each cell at any time step if I am using the DynOmegaD2Q9Descriptor.
Regards,
Junwei GuoAugust 2, 2019 at 10:44 pm in reply to: Core dumped using SmagorinskyPowerLawBGKdynamics for particle sedimenation #4491guojuwParticipantDear Mathias,
Do I need to modify the header/sources file under “${openlbRoot}/src/”, or just need to change the codes within the “int main(int argc, char* argv[]){}”.
Thanks,
Junwei -
AuthorPosts