[DOMaier]

Hey all,

I am working on the poro-elastoplasitc solvers by Tian Tang.
For a multitude if situations concerning geotechnical problems I need to specify a friction boundary condition.

The BC should either specify the displacement or traction (depending on the situation) for the normal component of the DU field.
For the tangential component DU should be 0 as long as the resulting shear stress doesn't exceed the maximum value given by Coulomb's friction law (maxShear proportional to normal traction).
If it does exceed, it should chance to a fixed (tangential) traction BC for that face and iteration.

I was looking at the solidContact Boundary condition, but i think it is a bit overkill for what i need, since i dont need the ability of detachment and re-attachment of the material (yet).

So i was thinking, I should be able to build on the fixedDisplacementZeroShear BC.
However i can not seem to find where the tangential components are specified.
I think the confusion stems from the whole valueFraction-issue that is quite hard to get the head around at the start.

Can someone point me to the right lines in the fixedDisplacementZeroShear to modify?

Greetings
Denis

[DOMaier]

I managed to gain a little bit more insight in how the boundary condition works, in that i now understand that valueFraction is a Tensor(Field) that tells me the mixture of fixedValue and fixedGradient in all directions. (I dont know if thats common practice but thats a beautiful solution)

so then I altered the updateCoeffs() function like this:




void fixedSoilDisplacementCoulombShearFvPatchVectorFiel d::updateCoeffs()
{
if (this->updated())
{
return;
}

vectorField n = patch().nf(); //face normals

const fvMesh& mesh = patch().boundaryMesh().mesh(); //mesh reference

const volSymmTensorField& sigmaV =
mesh.objectRegistry::lookupObject<volSymmTensorFie ld>("sigma"); // (effective) stress field

label patchID = patch().index();

symmTensorField sigmaP = sigmaV.boundaryField()[patchID]; // stress boundary field

vectorField t = sigmaP & n; // traction on boundary
scalarField smallField = scalarField(SMALL, sigmaP.size()) ;
scalarField nt = (t & n) + smallField ; // normal traction magnitude
vectorField st = t - ((t & n) * n) ; // shear traction

vectorField tt ;

forAll(patch(),iface)
{
if (mag(st[iface]) > (muFric_ * nt[iface]))
{
this->valueFraction()[iface] = sqr(n[iface]); // fixedValue only in n-direction
tt[iface] = (st[iface]/mag(st[iface])) * (muFric_ * nt[iface]); // friction traction
}
else
{ this->valueFraction()[iface] = I; // fixedValue in all directions
tt[iface] =0; // won't get used
}
}

bool incremental_ = 1;
scalarField pp = scalarField(sigmaP.size(), 0.0);


// gradient in DU for given traction
refGrad() = soilTractionBoundaryGradient::snGrad
(
tt, //traction
pp, //pressure
fieldName_,
"U",
patch(),
orthotropic_,
nonLinear_,
incremental_
);

directionMixedFvPatchVectorField::updateCoeffs();
}


What do u think about that, would that give me what i want or am I on a wrong track?
the scalar muFrac is read from dictionary as private data


It complies nicely but i get a segmentation error when running a simple case with the BC..
Probably i missed a pointer to data conversion somewhere?

[bigphil]

Hi,

I can see one possible cause of the segmentation fault, coming from the following line:

vectorField tt ;

You should initialise the size of the array otherwise the null constructor will set the size to zero i.e. it should be:

vectorField tt(sigmaP.size(), vector::zero);

Also, the following line:

scalarField smallField = scalarField(SMALL, sigmaP.size()) ;

should be changed to be:

scalarField smallField = scalarField(sigmaP.size(), SMALL) ;

as the field constructor requires SIZE as the first argument and the VALUE as the second.

I would also add a SMALL to the division when calculating the shear traction direction to avoid a floating point, and change "= 0" to "= vector::zero":

forAll(patch(),iface)
{
    if (mag(st[iface]) > (muFric_ * nt[iface])) 
    { 
    this->valueFraction()[iface] = sqr(n[iface]); // fixedValue only in n-direction
    tt[iface] = (st[iface]/(mag(st[iface]) + SMALL)) * (muFric_ * nt[iface]); // friction traction
}
else 
{
    this->valueFraction()[iface] = I; // fixedValue in all directions
    tt[iface] = vector::zero; // won't get used
}
}

Apart from that, it makes sense to me.

Philip

[DOMaier]

thank you a lot for the helpful remarks!


the last few days i was testing out the BC and i was making some minor changes.


The problem i have right now is that i get a snapping behavior of the friction boundary.
I think it comes from the fact, that i use the stress field from the old time step to determine the frictional stresses.
However, in the meantime i solved the DU-Eqn which is leading to a knew stress field.


Calculating an intermediate stress field and using that one to determine the frictional stresses destroys convergence of the solution completely.


Underrelaxation of the "intermediate" stress field doesn't help either.



Do you know how that problem is solved or different in the solidContact boundary?

[bigphil]

Hi Denis,

solidContact calculates the current stress and the friction condition is based on this; however, it uses a so-called penalty formulation (even stuck faces slip a small numerical amount) and does not use the valueFraction field.

But my feeling is that it should be possible to get it to converge using your approach; you may need to relax the valueFraction field as well as the value of the friction stress. Small under-relaxation values may be needed e.g. 0.01.

Philip