repulsionConstant estimation

[mojtabaAmiraslanpour]

Dear Hemocell community,

I was wondering if the default value set in the config.xml file (pipeflow) for repulsionConstant is in a normal range. It is set to be 2e-22 which if in lattice units, it would correspond to a minimal force value.

dx = 0.5e-6
dt = 1e-7
rho = 1e3

Considering the above values, the default repulsion force would be on the order of 1e-30 N which is way small. I was expecting a force value in the range of 1 pN.

Am I doing this right?

[gzavo]

Hi @mojtabaAmiraslanpour,
Indeed you want something in the pN order. The best answer for now is that you will need to look at the code and follow the value of repulsionConstant to make sure how it was implemented. It has not been used actively for a while, and is mostly there to show how to introduce attractive or repulsive forces technically. We might produce a calibrated example for this in the future, but atm it is not priority in any running projects. Btw these can act between cells, or if you put down artificial repulsive particles in the wall then between cells and wall. If you happen to create an example case for this, please do consider contributing it to the repo.

[mojtabaAmiraslanpour]

Thank you @gzavo,
This is the formula I found that is implemented in the code:

When d=< r_cutOff, this force will be activated. I have set F_Const=0.1 pN (1.6e-5 FU) and r_cutOff = 1 LU. So when particles are 0.1 LU apart, ~1 pN force is applied between them. It gives good results in terms of anti-penetration behavior:

However, it should not be too large which would alter the physical flow. This definitely needs more simulation cases with various repulsion forces considered to select the most efficient and correct value.

I would be more than happy to contribute to the repo if I come up with a good example and a comprehensive understanding of this!

[mojtabaAmiraslanpour]

The results of my simulations show that applying this repulsive force will prevent penetrated red blood cells to separate from each other. The definition of the force would not recognize between "penetrated" vs. "not penetrated" couple of cells and it would apply the force between any close nodes even if they are penetrated.

Can you please elaborate on this part of the code so that I can come up with a solution? is distance the distance between the two nodes?

HemoCell/core/hemoCellParticleField.cpp

Lines 677 to 743 in 3189621

	#define inner_loop \
	const int & l_index = grid_index(x,y,z); \
	const int & n_index = grid_index(xx,yy,zz); \
	for (unsigned int i = 0; i < particle_grid_size[l_index];i++){ \
	for (unsigned int j = 0; j < particle_grid_size[n_index];j++){ \
	HemoCellParticle & lParticle = particles[particle_grid[l_index][i]]; \
	HemoCellParticle & nParticle = particles[particle_grid[n_index][j]]; \
	if (&nParticle == &lParticle) { continue; } \
	if (lParticle.sv.cellId == nParticle.sv.cellId) { continue; } \
	const hemo::Array<T,3> dv = lParticle.sv.position - nParticle.sv.position; \
	const T distance = sqrt(dv[0]*dv[0]+dv[1]*dv[1]+dv[2]*dv[2]); \
	if (distance < r_cutoff) { \
	const hemo::Array<T, 3> rfm = r_const * (1/(distance/r_cutoff)) * (dv/distance); \
	lParticle.sv.force_repulsion = lParticle.sv.force_repulsion + rfm; \
	nParticle.sv.force_repulsion = nParticle.sv.force_repulsion - rfm; \
	} \
	} \
	}
	void HemoCellParticleField::applyRepulsionForce(bool forced) {
	const T r_const = cellFields->repulsionConstant;
	const T r_cutoff = cellFields->repulsionCutoff;
	if(!pg_up_to_date) {
	update_pg();
	}
	for (HemoCellParticle & particle : particles) {
	particle.sv.force_repulsion = {0.,0.,0.};
	}
	for (int x = 0; x < atomicLattice->getNx()-1; x++) {
	for (int y = 0; y < atomicLattice->getNy()-1; y++) {
	for (int z = 0; z < atomicLattice->getNz()-1; z++) {
	//Manual finding, we could make a map, but for now this should be fast enough
	//0, 0, 0 //0, 1, 0 //0, 0, 1 //0, 1, 1
	int xx = x;
	for (int yy = y; yy <= y+1; yy++) {
	for (int zz = z; zz <= z+1; zz++) {
	inner_loop
	}
	}
	//1, 0, 0
	//1, 1, 0
	//1, 0, 1
	//1, 1, 1
	//1, 0,-1
	//1,-1, 0
	//1,-1,-1
	//1, 1,-1
	//1,-1, 1
	xx = x+1;
	for (int yy = y-1; yy <= y+1; yy++) {
	for(int zz = z-1; zz <= z+1; zz++) {
	if (yy < 0) {continue;}
	if (zz < 0) {continue;}
	inner_loop
	}
	}
	xx = x;
	int yy = y+1, zz = z-1;
	if (zz<0) { continue; }
	//0, 1,-1
	inner_loop
	}
	}
	}
	}

[gzavo]

You can see how distance is defined in the code-segment you linked:

HemoCell/core/hemoCellParticleField.cpp

Line 687 in 3189621

const T distance = sqrt(dv[0]*dv[0]+dv[1]*dv[1]+dv[2]*dv[2]); \

This is the distance of surface particles, with no notion of whether penetrating or not.
Checking for every particle if they are inside another cell would have a huge associated computational cost.