In D3Q27_PF_VELOCITY model , how to calculate the forces , viscous force,pressure force in paraview?

[engparth-cmd]

I'm using the D3Q27_PF_VELOCITY model to simulate a droplet impacting an inclined surface. I want to calculate the surface-tension and gravitational forces as the droplet slides. To see these values in ParaView, should I add a global function in Dynamics.R or i will have to implement it in .xml file using R ?

[llaniewski]

I assume @TravisMitchell will be the best person to answer.

But general notes: Global are global integrals over the domain. So they are used to monitor some single thing over time. For example, mass of the droplet to see if stuff is not dissapearing, etc. If you want to see something in Paraview, you probably mean a Quantity, which are fields exported to VTI files you open in Paraview.

Can you elaborate what you mean by the surface-tension forces here?

[engparth-cmd]

sir i am attaching here the viscous dissipation energy equation where all the term are mentioned here

[engparth-cmd]

Ed is dissipation energy and Vd represnts the volume of droplet and Vg represents the volume of air and W(t) is calculated as

[engparth-cmd]

[TravisMitchell]

Hi @engparth-cmd, if I understand these equations correctly -- I do not think you need to add vtk/hdf5 outputs for the forcing terms. This looks like a global quantity integrated over a portion of the domain Vd U Vg. As such, you could create an AddGlobal(W_t), then to get the velocity gradients -- you could consider if these can be determined locally from the distribution functions, or add stencils to the velocity fields and in Dynamics.c.Rt AddToW_t(calculated_value_in_integrand)

[engparth-cmd]

Thank you sir i will try this for calculating the viscous dissipiation energy.
Apart from this as our earlier discussion for calculating the body force,pressure force,surface tension force and viscous force , i have written in following way and it is working.Sir can you please check it is properly implemented?
In dynamics.R i have implemented
AddQuantity(name="Fs", unit="N", vector=T)
AddQuantity(name="Fp", unit="N", vector=T)
AddQuantity(name="Fb", unit="N", vector=T)
AddQuantity(name="Fm", unit="N", vector=T)

and in Dynamics.c.Rt in the following section i have implemented forces as

Get functions to define outputs indicated in .R

// CALCULATION OF BODY FORCE

CudaDeviceFunction vector_t getFb() {
real_t Fx, Fy, Fz;

// local density at the current lattice node:
real_t rh = getRho();       

// compute body-force components:
calc_Fb(&Fx, &Fy, &Fz,rh);

// return as a vector field:
return (vector_t){ Fx, Fy, Fz};

}

/// CALCULATION OF PRESSURE FORCE

CudaDeviceFunction vector_t getFp() {
real_t Fpx, Fpy, Fpz;
real_t pressure = getP();
vector_t gPhi=calcGradPhi();
// compute body-force components:
calc_Fp(&Fpx, &Fpy, &Fpz, pressure, gPhi);

// return as a vector field:
return (vector_t){Fpx, Fpy, Fpz};

}

// CALCULATION OF SURFACE TENSION FORCE

CudaDeviceFunction vector_t getFs() {
real_t Fsx, Fsy, Fsz;
real_t C = getPhaseField();
real_t mu = calcMu(C);
vector_t gPhi=calcGradPhi();
// compute body-force components:
calc_Fs(&Fsx, &Fsy, &Fsz,mu, gPhi);

// return as a vector field:
return (vector_t){Fsx, Fsy, Fsz};

}

// CALCULATION OF VISCOUS FORCE

CudaDeviceFunction vector_t getFm() {
// --- locals from node state
real_t C = getPhaseField();
real_t tau = calcTau(C);
vector_t gradPhi = calcGradPhi();

// ❗ provide u and p for the R-template (it expects these names)
vector_t u = getU();   // current velocity at node
real_t p   = getP();   // pressure (alternatively, we can set p = m0[0] inside R)

// work arrays
real_t stress[6] = {0,0,0,0,0,0};
real_t m0[27], m[27];

<?R
    EQ   = MRT_eq(U, PV(1), u, mat=t(M))
    selR = EQ$order < 10

    # Use the C variable 'p' we declared above:
    EQ$Req[1] = PV("p")

    # load current moments m0 from distributions g
    m0sym = PV("m0[",1:27-1,"]")
    C(m0sym[selR], (M %*% g)[selR])

    # (optional) If you prefer p = m0[0], uncomment the next line and
    # comment 'real_t p = getP();' above:
    # C(PV("p"), m0sym[1])

    # post-collision moments (non-eq part we need for stresses)
    msym = PV("m[",1:27-1,"]")
    Omega      = PV(rep(1, sum(selR)))
    Omega[5:9] = PV("tau")^(-1)
    Omega[1:4] = PV(1)
    C(msym[selR], (m0sym - EQ$Req)[selR] * Omega)

    # back to "populations" representing the non-eq part we project
    new_g = solve(M) %*% msym

    st = PV("stress[",1:6-1,"]")
    C(st[1], sum(U[,1]*U[,1] * new_g)); # XX
    C(st[2], sum(U[,1]*U[,2] * new_g)); # XY
    C(st[3], sum(U[,1]*U[,3] * new_g)); # XZ
    C(st[4], sum(U[,2]*U[,2] * new_g)); # YY
    C(st[5], sum(U[,2]*U[,3] * new_g)); # YZ
    C(st[6], sum(U[,3]*U[,3] * new_g)); # ZZ
?>

// viscous force
real_t dRho = (Density_h - Density_l);
real_t Fx = (0.5 - tau) * dRho * (stress[0]*gradPhi.x + stress[1]*gradPhi.y + stress[2]*gradPhi.z);
real_t Fy = (0.5 - tau) * dRho * (stress[1]*gradPhi.x + stress[3]*gradPhi.y + stress[4]*gradPhi.z);
real_t Fz = (0.5 - tau) * dRho * (stress[2]*gradPhi.x + stress[4]*gradPhi.y + stress[5]*gradPhi.z);

return (vector_t){Fx, Fy, Fz};

}

[engparth-cmd]

F_mu[0] = (0.5-tau) * (Density_h-Density_l) * (stress[0]*gradPhi.x + stress[1]*gradPhi.y + stress[2]*gradPhi.z);
F_mu[1] = (0.5-tau) * (Density_h-Density_l) * (stress[1]*gradPhi.x + stress[3]*gradPhi.y + stress[4]*gradPhi.z);
F_mu[2] = (0.5-tau) * (Density_h-Density_l) * (stress[2]*gradPhi.x + stress[4]*gradPhi.y + stress[5]*gradPhi.z);
F_total[0] = F_surf[0] + F_pressure[0] + F_body[0] + F_mu[0];
F_total[1] = F_surf[1] + F_pressure[1] + F_body[1] + F_mu[1];
F_total[2] = F_surf[2] + F_pressure[2] + F_body[2] + F_mu[2];

	        AddToF_pressureX(F_pressure[0]);
		AddToF_pressureY(F_pressure[1]);
		AddToF_pressureZ(F_pressure[2]);
		AddToF_surfX(F_surf[0]);
		AddToF_surfY(F_surf[1]);
		AddToF_surfZ(F_surf[2]);
		AddToF_muX(F_mu[0]);
	 	AddToF_muY(F_mu[1]);
	        AddToF_muZ(F_mu[2]);	  

i have added these things in the force section for calculating the forces in the paraview.So does it work?

[TravisMitchell]

Hi @engparth-cmd do you have your own Fork of the code that I can view this on? I think one thing to consider is in the getFm function, there is a dependence of F on U and U also depends on F -> in the code you can see this that there is a fixed point iterator that runs this calculation a couple of times to try and cater for the implicit relationship here. I think you have perhaps removed this in your function above.

[engparth-cmd]

Yes sir, i have created my Fork with the name engparth-cmd. I have included there Globals in Dynamics.c.Rt and Dynamics.R please suggest me for the implementation.Kindly you can check in my Fork where i have included Globals in Dynamics.c.Rt and Dynamics.R. In the above discussion earlier i had tried by using getFM function but as it is integrated over domain so as per your sugggestion i have created Globals for visualizing the forces in the paraview.

[engparth-cmd]

F_mu[0] = (0.5-tau) * (Density_h-Density_l) * (stress[0]*gradPhi.x + stress[1]*gradPhi.y + stress[2]*gradPhi.z);
F_mu[1] = (0.5-tau) * (Density_h-Density_l) * (stress[1]*gradPhi.x + stress[3]*gradPhi.y + stress[4]*gradPhi.z);
F_mu[2] = (0.5-tau) * (Density_h-Density_l) * (stress[2]*gradPhi.x + stress[4]*gradPhi.y + stress[5]*gradPhi.z);
F_total[0] = F_surf[0] + F_pressure[0] + F_body[0] + F_mu[0];
F_total[1] = F_surf[1] + F_pressure[1] + F_body[1] + F_mu[1];
F_total[2] = F_surf[2] + F_pressure[2] + F_body[2] + F_mu[2];

        AddToF_pressureX(F_pressure[0]);
	AddToF_pressureY(F_pressure[1]);
	AddToF_pressureZ(F_pressure[2]);
	AddToF_surfX(F_surf[0]);
	AddToF_surfY(F_surf[1]);
	AddToF_surfZ(F_surf[2]);
	AddToF_muX(F_mu[0]);
 	AddToF_muY(F_mu[1]);
        AddToF_muZ(F_mu[2]);	  

i have added these things in the force section for calculating the forces in the paraview.So does it work?can you please check ?

[TravisMitchell]

hi @engparth-cmd - it is end of semester marking time for me, so I won't get time to test out your implementation during November. However, you are adding to Globals above, so won't be able to see these in Paraview - their integration over the domain would be added to a .csv file if you are using a Log element in the xml input file.

By the looks you have also tried to add these as quantities and added get functions for them, I would suggest you try to design some simple test cases to verify the implementation is working as you expect. Have you tried doing this already?