Skip to content

guojuw

Forum Replies Created

Viewing 15 posts - 1 through 15 (of 15 total)
  • Author
    Posts
  • guojuw
    Participant

    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 Guo

    guojuw
    Participant

    Dear 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 Guo

    guojuw
    Participant

    Thanks Jonathan. I worked it out by outputting the WSS[0] to WSS[2].

    Regards,
    Junwei

    in reply to: GSL – GNU Scientific Library #5249
    guojuw
    Participant

    Dear 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 GUO

    guojuw
    Participant

    Dear 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 GUo

    in reply to: SmoothIndicatorEllipsoid3D #5189
    guojuw
    Participant

    Dear 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;
    }

    in reply to: 2D and 3D HLBM object contact stiffness #4683
    guojuw
    Participant

    Dear 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 Guo

    in reply to: 2D and 3D HLBM object contact stiffness #4680
    guojuw
    Participant

    Dear Mathias,

    I mean the parameter gamma in equation(9). How to modify this particle collision parameter?

    Regards,
    Junwei Guo

    Article (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-13

    in reply to: Bug: void ParticleDynamics3D::addSphere() #4646
    guojuw
    Participant

    Dear 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 Guo

    in reply to: Problem in example code porousPoiseuille2d #4573
    guojuw
    Participant

    Dear 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 Guo

    in reply to: Problem in example code porousPoiseuille2d #4568
    guojuw
    Participant

    Dear 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.

    velocity profiles

    error analysis

    Regards,
    Junwei Guo

    guojuw
    Participant

    Dear 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 Guo

    Chai2011 – Multiple Relaxation Time Lattice Boltzmann Model for Generalized Newtonian Fluid Flows.pdf

    in reply to: Printing the relaxation time of each lattice to vtk files #4523
    guojuw
    Participant

    Dear Aurelio,

    Great! Thanks a lot. It helps me.

    Regard,
    Junwei Guo

    in reply to: Printing the relaxation time of each lattice to vtk files #4516
    guojuw
    Participant

    Hi 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 Guo

    guojuw
    Participant

    Dear 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

Viewing 15 posts - 1 through 15 (of 15 total)