[thejavarma9]

Originally Posted by amuzeshi

This surface is not oriented and positioned exactly as the left end of your whole geometry; It must be so in order to do cell operations.

Hi!
Yeah, sorry that was a mistake i haven't noticed. But the problem is even the mesh and stl file positioned exactly the same error (outside point is not inside any cell) continued. I tried various points (outside Points) which lies on outer surface of the component to include the inside surface. Could you help me where the mistake was! I have attached mesh file and stl files.


https://drive.google.com/drive/folde...Wf?usp=sharing

[amuzeshi]

Originally Posted by thejavarma9

Hi!
Yeah, sorry that was a mistake i haven't noticed. But the problem is even the mesh and stl file positioned exactly the same error (outside point is not inside any cell) continued. I tried various points (outside Points) which lies on outer surface of the component to include the inside surface. Could you help me where the mistake was! I have attached mesh file and stl files.


https://drive.google.com/drive/folde...Wf?usp=sharing

To the best of my knowledge, 'surfaceToCell' works on volumetric mesh whereas yours is a surface mesh.

[CRI_CFD]

Dear all,


Silly questions: is elasticSolidFoam able to deal with multi-domain, multi-material problems with contact regions?


I am unable to make it run, once I got the meshes by means of snappy+setSet+splitMeshRegions...


Other solvers (chtMultiRegionFoam) allow to remove the contant/polyMesh and it works with meshses defined by regions. elasticSolidFoam crashes once I remove polyMesh (original mesh before splitMeshRegions) folder.


Any hints?


Thanks in advance

[bigphil]

Originally Posted by CRI_CFD

Dear all,


Silly questions: is elasticSolidFoam able to deal with multi-domain, multi-material problems with contact regions?


I am unable to make it run, once I got the meshes by means of snappy+setSet+splitMeshRegions...


Other solvers (chtMultiRegionFoam) allow to remove the contant/polyMesh and it works with meshses defined by regions. elasticSolidFoam crashes once I remove polyMesh (original mesh before splitMeshRegions) folder.


Any hints?


Thanks in advance

What do you mean by "multi-domain"? elasticSolidFoam does not use "regions"; instead it treats multi-materials via the global mesh and a material indicator field (e.g. 0 means the cell is in material 0, 1 means the cell is in material 1, etc.).

The multi-material procedure is implemented differently in solids4foam (https://bitbucket.org/philip_cardiff...s4foam-release), where it uses cell zones to define different materials; you may find that easier to use.

[CRI_CFD]

Originally Posted by bigphil

What do you mean by "multi-domain"? elasticSolidFoam does not use "regions"; instead it treats multi-materials via the global mesh and a material indicator field (e.g. 0 means the cell is in material 0, 1 means the cell is in material 1, etc.).

The multi-material procedure is implemented differently in solids4foam (https://bitbucket.org/philip_cardiff...s4foam-release), where it uses cell zones to define different materials; you may find that easier to use.

Hi big_phil,


Thanks for your answer.


I meant I have a multi-region domain, after meshing two stl geometries with an interface zone in common with snappy (I don't know the interface region in advance). After that I use cellSet + splitMeshRegions to remove background mesh from blockMesh and keep only my two regions of interest.


After that, I can easily assign materials to regions using materials file and cellToRegion content on constat folder (paying attention to the number of materials, as of course one of them is not present anymore in the meshed domain)...


The thing is that (I read the code yesterday) elasticSolidFoam doesn't work with regions (there is no vector solids[i] or fluids[i] as in chtMultiRegionFoam). So if I want to proceed, I will have to code by myself the solver including the presence of several solid regions, following the style of chtMultiRegionFoam...



I am also open to use solids4foam but...still can't find the way to deal with the interface to be able to assign a contact condition on it. The tutorial slidingFrictionBall is based on an easy geometry created with blockMesh (and after with setSet the contact zones are assigned, but I need snappy because of the complexity of the geometry).


Here is a link to my case, just in case you or someone can give me any hints, I feel really stuck into this now (I am planning to use elasticSolidFoam and elasticPlasticSolidFoam with a contact model at the interface between both regions):


https://1drv.ms/u/s!AjenQBoiljuzgoQs...rRL4Q?e=KqbJla


Many thanks in advance! Any idea will be really appraciated

[bigphil]

Have you seen the "solids4foam-release/tutorials/solids/linearElasticity/punch" tutorial?
It demonstrates a contact case for two objects that come from STLs. In that case, cfMesh is used but it is the same method when using snappyHexMesh.

[CRI_CFD]

Originally Posted by bigphil

Have you seen the "solids4foam-release/tutorials/solids/linearElasticity/punch" tutorial?
It demonstrates a contact case for two objects that come from STLs. In that case, cfMesh is used but it is the same method when using snappyHexMesh.

Hi bigphil,


I was checking it now and applying it to my case (with foam extend, but reading the Allrun seems that we are doing the same, both solids4foam and foam extend 4.1, the tutorial creates two meshes with defined patches in the geometry, reasign patches and merge meshes). Problem still remains:


Create time

Create mesh for time = 0


Reading g
Reading field U

Patch load Traction boundary field: U
--> FOAM Warning : The gradScheme for U should be "extendedLeastSquares 0" for the boundary non-orthogonal correction to be right
Patch oldInternalFaces Traction boundary field: U
--> FOAM Warning : The gradScheme for U should be "extendedLeastSquares 0" for the boundary non-orthogonal correction to be right
Creating solidContact patch
Normal contact model: dirichletNeumann
Friction contact model: dirichletNeumann
Friction law: coulomb
Interpolation of traction from slave to master: ggi
--> FOAM Warning : The gradScheme for U should be "extendedLeastSquares 0" for the boundary non-orthogonal correction to be right
Creating solidContact patch
Selecting rheology model linearElastic
Creating constitutive model
--> FOAM Warning : history patch support_mandible not found. Force-displacement will not be written
Selecting divSigmaExp calculation method surface

Starting time loop

Time = 1


Predicting U, gradU and snGradU based on V,gradV and snGradV

Slave contact patch esfera_acoplada1_to_esfera_completa1 grabbing normalContactModel pointer from master
Slave contact patch esfera_acoplada1_to_esfera_completa1 grabbing frictionContactModel pointer from master


--> FOAM FATAL ERROR:

request for regIOobject solidInterface from objectRegistry region0 failed
available objects of type regIOobject are
0()

From function objectRegistry::lookupObject<Type>(const word&) const
in file /home/isma/foam/foam-extend-4.1/src/foam/lnInclude/objectRegistryTemplates.C at line 198.

FOAM aborting

Aborted (core dumped)





I think it only works when you clearly define your patches within your stl/ftr/geometry files (I suppose a perfect matchs is needed, I mean region1.stl with interface_master patch and region2.stl with interface_slave patch mathing perfectly with interface_master). But I do not know my interface regions in advance, and also number of triangles of master and slave can differ.



The crucial step for me is when I create regions using splitMeshRegions. This is the point where interfaces are automatically created. But then I have 2 regions and elasticSolidFoam fails (works only in 1 region as you pointed). I tried to reverse the procdure using mergeMeshes, and even if the interface is kept, when I run checkMesh I notice I still have two domains.


So my guesses are:


1. To find a way to obtain my interfaces from cellZones (without the need of running splitMeshRegions) after setSet+subsetMesh. Thus, I will have 1 region and elasticSolidFoam should work properly. Also contacts can be defined properly, as the interfaces are available under constant/boundary

2. To find a way to merge regions created with splitMeshRegions without removing interfaces. Hence, 1 region will be available, but interfaces as well, so I can apply contact conditions in 0/U
3. To keep the idea of splitMeshRegions and modify elasticSolidFoam to be able to run in several regions (sounds too complex for my OpenFOAM level)



What do you think? I am trying with a simplified case of 2 spheres in contact (I created them with MeshMixer)

[CRI_CFD]

Originally Posted by bigphil

Have you seen the "solids4foam-release/tutorials/solids/linearElasticity/punch" tutorial?
It demonstrates a contact case for two objects that come from STLs. In that case, cfMesh is used but it is the same method when using snappyHexMesh.

BTW bigphil,


I have run the punch tutorial as you mentioned (with the Allrunscript).


When I check the mesh:


/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
Build : 4.1-70b064d0f326
Exec : checkMesh
Date : May 16 2022
Time : 19:49:55
Host : "isma-Super-Server"
PID : 28091
CtrlDict : "/home/isma/foam/solids4foam/solids4foam-release/tutorials/solids/linearElasticity/punch/system/controlDict"
Case : /home/isma/foam/solids4foam/solids4foam-release/tutorials/solids/linearElasticity/punch
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
all points: 1173
live points: 1173
all faces: 2943
live faces: 2943
internal faces: 2303
cells: 900
boundary patches: 8
point zones: 0
face zones: 0
cell zones: 2

Overall number of cells of each type:
hexahedra: 774
prisms: 28
wedges: 0
pyramids: 70
tet wedges: 0
tetrahedra: 28
polyhedra: 0

Checking topology...
Boundary definition OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
*Number of regions: 2
The mesh has multiple regions which are not connected by any face.
<<Writing region information to "0/cellToRegion"
Nuumber of cells per region:
0 120
1 780


Checking patch topology for multiply connected surfaces ...
Patch Faces Points Area [m^2] Surface topology
punchContact 4 9 0.00114293 ok (non-closed singly connected)
cylinderTractionFree200 229 0.0313897 ok (non-closed singly connected)
symmx 142 173 0.0248622 ok (non-closed singly connected)
cylinderContact 48 59 0.00783092 ok (non-closed singly connected)
punchLoading 16 22 0.00194111 ok (non-closed singly connected)
cylinderFixed 48 59 0.00783092 ok (non-closed singly connected)
symmy 142 173 0.0248605 ok (non-closed singly connected)
punchTractionFree 40 53 0.00741268 ok (non-closed singly connected)

Checking geometry...
This is a 3-D mesh
Overall domain bounding box (0 0 -0.2) (0.1 0.1 0.1)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Mesh (non-empty, non-wedge) dimensions 3
Boundary openness (1.09252e-16 -2.27411e-17 -2.46282e-17) Threshold = 1e-06 OK.
Max cell openness = 2.5034e-16 OK.
Max aspect ratio = 5.11641 OK.
Minumum face area = 8.73766e-06. Maximum face area = 0.000326572. Face area magnitudes OK.
Min volume = 3.22562e-08. Max volume = 3.80673e-06. Total volume = 0.00174183. Cell volumes OK.
Mesh non-orthogonality Max: 33.4544 average: 7.46978 Threshold = 70
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 0.650548 OK.

Mesh OK.

End









So after simulation, mesh still have two regions....so I guess solids4foam overperforms elasticSolidFoam in terms of management of regions (elasticSolidFoam only works with one region, on the contrary, solids4foam is able to run with 2 regions).


Is this correct? If so, I will try to run my mesh with solids4foam and see if it works...

[bigphil]

I think there is some confusion around the word "regions" here.

Solvers like chtMultiRegionSimpleFoam use multiple regions (e.g. solid and fluid) where each region is a separate mesh stored in, e.g. constant/solid/polyMesh and constant/fluid/polyMesh. When these types of solvers run, they write the results to the individual regions, e.g. 100/fluid/U and 100/solid/T.
Similarly, solids4foam uses this approach for FSI cases. In that case, there is no global mesh, i.e. constant/polyMesh does not exist or is empty.
Sometimes these "regions" are called sub-meshes.

When you run checkMesh on a case with regions like this then checkMesh will ignore these regions by default; if you run checkMesh on a solids4foam FSI case then it will complain that no mesh exists. Instead you would need to use the "-region" option of checkMesh to check an individual region, e.g. checkMesh -region solid.

Now, confusingly, when you run checkMesh on any mesh, whether it is the global mesh or a region, checkMesh will examine the connectivity of the cells. If it finds groups of cells which are disconnected then it reports these as "Number of regions:". These "regions" are just groups of cells within the mesh and are not mesh regions as per above. These are not regions as per multi-region solvers. Does that make sense?

So, to answer your question, for solid-only simulations (not FSI) then both elasticSolidFoam and solids4foam use only one region, as per the definition of regions/sub-meshes above; however, both elasticSolidFoam and solids4foam are able to deal with multiple disconnected groups of cells within the mesh. For example, the slidingFrictionBall case works with both solvers and contains only one region/sub-mesh (the default constant/polyMesh region), where there are two disconnected groups of cells as per the "Number of regions:" checkMesh output.

[CRI_CFD]

Originally Posted by bigphil

I think there is some confusion around the word "regions" here.

Solvers like chtMultiRegionSimpleFoam use multiple regions (e.g. solid and fluid) where each region is a separate mesh stored in, e.g. constant/solid/polyMesh and constant/fluid/polyMesh. When these types of solvers run, they write the results to the individual regions, e.g. 100/fluid/U and 100/solid/T.
Similarly, solids4foam uses this approach for FSI cases. In that case, there is no global mesh, i.e. constant/polyMesh does not exist or is empty.
Sometimes these "regions" are called sub-meshes.

When you run checkMesh on a case with regions like this then checkMesh will ignore these regions by default; if you run checkMesh on a solids4foam FSI case then it will complain that no mesh exists. Instead you would need to use the "-region" option of checkMesh to check an individual region, e.g. checkMesh -region solid.

Now, confusingly, when you run checkMesh on any mesh, whether it is the global mesh or a region, checkMesh will examine the connectivity of the cells. If it finds groups of cells which are disconnected then it reports these as "Number of regions:". These "regions" are just groups of cells within the mesh and are not mesh regions as per above. These are not regions as per multi-region solvers. Does that make sense?

So, to answer your question, for solid-only simulations (not FSI) then both elasticSolidFoam and solids4foam use only one region, as per the definition of regions/sub-meshes above; however, both elasticSolidFoam and solids4foam are able to deal with multiple disconnected groups of cells within the mesh. For example, the slidingFrictionBall case works with both solvers and contains only one region/sub-mesh (the default constant/polyMesh region), where there are two disconnected groups of cells as per the "Number of regions:" checkMesh output.

Hi bigphil,

Thanks for your detailed explanations. And as you mention sounds weird same software is calling "regions" to two different things.

After understanding this, I am even more confused about my case...

Basically, after meshing with snappy my 2 regions (block and lower_teeth) I have the whole thing meshed (block, lower_teeth and background coming from blockMesh). I need to remove the background cells (I call this region "air"). Hence, I have 3 cellZones.

block and lower_teeth are in contact, so I also need an interface to impose the contacts.

So the procedure I considered was to isolate my tow regions of interest using cellSet, then subsetMesh and then splitMeshRegions (of course, different output, depending if I use cellZones or cellZOnesOnly options). In any case following this procedure I got 2 regions using cellZonesOnly (two folders inside constant directory, each one with a meshed region as you mentioned). Hence, in order to group everything, I split my case in two cases (case 1 with polyMesh coming from block folder and case 2 with polyMesh coming from lower_teeth folder), removing the original polyMesh. Then, using mergeMeshes, I recover (in a folder callled "1") the "merged" mesh with created interfaces and patches of both meshes, so contacts can be applied in 0/D (or 0/U in elasticSolidFoam). When I ran the case, during execution I can see that contacts are being imposed. But, for example, in solids4foam, following warning arises:



Create time

/*---------------------------------------------------------------------------*\
| For further information on the solids4Foam toolbox implementations, |
| please see the following publications: |
| |
| P. Cardiff, A Karac, P. De Jaeger, H. Jasak, J. Nagy, A. Ivankovic, |
| Z. Tukovic: An open-source finite volume toolbox for solid mechanics and |
| fluid-solid interaction simulations. arXiv:1808.10736v2, 2018, available |
| at https://arxiv.org/abs/1808.10736. |
| |
| Z. Tukovic, A. Karac, P. Cardiff, H. Jasak, A. Ivankovic: OpenFOAM |
| finite volume solver for fluid-solid interaction. Transactions of |
| Famena, 42 (3), pp. 1-31, 2018, 10.21278/TOF.42301. |
\*---------------------------------------------------------------------------*/

Selecting physicsModel solid

Selecting solidModel linearGeometryTotalDisplacement
Selecting dynamicFvMesh staticFvMesh
Creating fixedDisplacement boundary condition
Creating solidTraction boundary condition
limiter coefficient: 1
under-relaxation method: fixed
Creating the mechanicalModel
Number of cells in new mesh : 2931583
Number of faces in new mesh : 9692536
Number of points in new mesh: 3829338
Cannot read "/home/consolids4foam/unidas/system/blockSubset". Copy from base
Cannot read "/home/consolids4foam/unidas/system/blockSubset". Copy from base
Number of cells in new mesh : 545699
Number of faces in new mesh : 1882096
Number of points in new mesh: 791502
Cannot read "/home/consolids4foam/unidas/system/lower_teethSubset". Copy from base
Cannot read "/home/consolids4foam/unidas/system/lower_teethSubset". Copy from base
Selecting mechanical law linearElastic
Selecting mechanical law linearElastic
Time = 1

Evolving solid solver
Solving the momentum equation for D
Reading individual shadowPatch
Normal contact model: standardPenalty
Initialising stored previous normalPenaltyFactors
Friction contact model: standardPenalty
Friction law: coulomb
solidContact: block_to_lower_teeth
useNewPointDistanceMethod: false
projectPointsToPatchBoundary: false
usePrevCandidateMasterNeighbors: false
--> FOAM Warning :
From function newGGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file lnInclude/newGGIInterpolationWeights.C at line 463
polygonIntersection is returning a zero surface area between
Master face: 35 and Neighbour face: 224 intersection area = 0
Please check the two quick-check algorithms for newGGIInterpolation. Something is missing.
--> FOAM Warning :
From function newGGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file lnInclude/newGGIInterpolationWeights.C at line 463
polygonIntersection is returning a zero surface area between
Master face: 222 and Neighbour face: 224 intersection area = 0
Please check the two quick-check algorithms for newGGIInterpolation. Something is missing.
--> FOAM Warning :
From function newGGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file lnInclude/newGGIInterpolationWeights.C at line 463
polygonIntersection is returning a zero surface area between
Master face: 224 and Neighbour face: 35 intersection area = 0
Please check the two quick-check algorithms for newGGIInterpolation. Something is missing.
--> FOAM Warning :
From function newGGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file lnInclude/newGGIInterpolationWeights.C at line 463
polygonIntersection is returning a zero surface area between
Master face: 224 and Neighbour face: 222 intersection area = 0
Please check the two quick-check algorithms for newGGIInterpolation. Something is missing.

Simulation is running (I can see residuals with correctors), but after some time, I get a floating point exception. So at this point I have no idea if the problem comes from the mesh or from my set up....I found this the easiest way to create the interfaces (they are not included within stl files).

When I started, I found that snappy , setSet...from OF 9 release is working better to remove the background mesh and differentiate my two regions...So I do everything with OF9 and then I switch to fe4.1/solids4foam. I suppose this can't be the mistake, but no 100% sure...

With elasticSolidFoam (when I just did splitmeshregions use, I can't see that the solver is considering the contact regions during runtime, so seems it is using original polyMesh, without interfaces, so usually it runs for a few iterations and then max stress/strain becomes huge)...but in this case, after split and mergeMeshes, the output is:


Create time

Create mesh for time = 0


Reading g
Reading field U

--> FOAM Warning : The gradScheme for U should be "extendedLeastSquares 0" for the boundary non-orthogonal correction to be right
Patch oldInternalFaces Traction boundary field: U
--> FOAM Warning : The gradScheme for U should be "extendedLeastSquares 0" for the boundary non-orthogonal correction to be right
Creating solidContact patch
Normal contact model: standardPenalty
Friction contact model: frictionless
Interpolation of traction from slave to master: ggi
Patch load Traction boundary field: U
--> FOAM Warning : The gradScheme for U should be "extendedLeastSquares 0" for the boundary non-orthogonal correction to be right
Creating solidContact patch
Selecting rheology model linearElastic
Creating constitutive model
Force-displacement for patch support will be written to forceDisp.dat
Selecting divSigmaExp calculation method surface

Starting time loop

Time = 1


Predicting U, gradU and snGradU based on V,gradV and snGradV

Slave contact patch lower_teeth_to_block grabbing normalContactModel pointer from master
Slave contact patch lower_teeth_to_block grabbing frictionContactModel pointer from master


--> FOAM FATAL ERROR:

request for regIOobject solidInterface from objectRegistry region0 failed
available objects of type regIOobject are
0()

From function objectRegistry::lookupObject<Type>(const word&) const
in file /home/isma/foam/foam-extend-4.1/src/foam/lnInclude/objectRegistryTemplates.C at line 198.

FOAM aborting

Aborted (core dumped)



I would appreciate any tip or idea. Would be great to know if I am failing in the mesh generation, the solver set-up or the contact regions definition and modelling. I can share simulation files if needed.


Many thanks

[bigphil]

Your mesh workflow seems fine to me.

I suggest sticking with solids4foam.

To check the meshes are "good", I suggest using a fixedDisplacement patch on each of the two regions and then applying constant traction to some other patch on each of the two reasons. This should work if the mesh is fine.

Then, you can switch to using solidContact conditions. When using these solidContact conditions, it is important that rigid body motions are constrained, otherwise one (or both) of the parts may "float away" causing divergence. The best setup for convergence would be to have fixedDisplacement patches on both regions, which displace them into contact.

[CRI_CFD]

Originally Posted by bigphil

Your mesh workflow seems fine to me.

I suggest sticking with solids4foam.

To check the meshes are "good", I suggest using a fixedDisplacement patch on each of the two regions and then applying constant traction to some other patch on each of the two reasons. This should work if the mesh is fine.

Then, you can switch to using solidContact conditions. When using these solidContact conditions, it is important that rigid body motions are constrained, otherwise one (or both) of the parts may "float away" causing divergence. The best setup for convergence would be to have fixedDisplacement patches on both regions, which displace them into contact.

Hi big phil,


I will try your strategy and let you know. I really want to use solids4foam, looks amazing!


Just to make sure about your suggestion:


1. Change solidContact BC on interface (replace them by fixedDisplacement)
2. I have only one patch for imposing loads, so I will use here solidTraction
3. I am fixing the geometry on the bottom patch using fixedDisplacement to (0,0,0) trying to avoid the geometry to float away
4. All other patches are solidTractionFree.


Before doing the suggested changes, this is my D file:


support
{
type fixedDisplacement;
value uniform (0 0 0);
}
load
{
type solidTraction;
traction uniform ( 0 0 10000 );
pressure uniform 0;
value uniform (0 0 10000);
}
oldInternalFaces
{
type solidTractionFree;
value uniform (0 0 0);

}
block_to_lower_teeth
{
type solidContact;
master yes;
contactActive yes;
rigidMaster no;
shadowPatch lower_teeth_to_block;
interpolationMethod ggi;
//interpolationMethod patchToPatch;
projectionAlgo visible;
projectionDir contactSphere;
//projectionDir vector;
correctionFrequency 1;
normalContactModel standardPenalty;
standardPenaltyNormalModelDict
{
relaxationFactor 0.1;
penaltyScale 1;
returnScale 1.5;
contactGapTol 1e-7;
limitPenetration no;
penetrationLimit -0.03;
limitPressure no;
pressureLimit 1e12;
correctMissedVertices yes;
distanceMethod point;
//distanceMethod face;
oscillationCorrection yes;
smoothingSteps 1;
infoFrequency 10;
}
//frictionContactModel frictionless;
frictionContactModel standardPenalty;
standardPenaltyFrictionModelDict
{
relaxationFactor 0.1;
penaltyScale 1;
oscillationCorrection yes;
smoothingSteps 10;
infoFrequency 10;

frictionLaw coulomb;
frictionLawDict
{
frictionCoeff 0.1;
}
}

value uniform (0 0 0);
}

lower_teeth_to_block
{
type solidContact;
master no;
contactActive yes;
shadowPatch block_to_lower_teeth;
value uniform (0 0 0);
}
}








so after your suggestions:


support
{
type fixedDisplacement;
value uniform (0 0 0);
}
load
{
type solidTraction;
traction uniform ( 0 0 10000 );
pressure uniform 0;
value uniform (0 0 10000);
}
oldInternalFaces
{
type solidTractionFree;
value uniform (0 0 0);

}
block_to_lower_teeth
{
type fixedDisplacement;
value uniform (0 0 0);
}

lower_teeth_to_block
{
type fixedDisplacement;
value uniform (0 0 0);
}
}


Is this what you are suggesting?



Fingers cross...hope the mesh is fine!


Thanks for the advise!

[CRI_CFD]

Well, I have tried several things during this time, based on bigphil suggestions.


Here are my findings (if any).


1. I tested a simplified geometry with 3 zones and 3 contact regions. I set all interfaces to a fixedValue BC (value (0 0 0 )). If I go for a serial ran, I get the log attached. Everything seems to converge properly, so...does this mean the mesh is fine? I would say yes, but...


2. If I go for a parallel run, with exactly the same settings I got the following log (different and no convergence is reached)


3. To be honest, the step 2 is not exactly the same case, as I needed to specify a value for my solidTraction BC on oldInternal faces (actually the outer boundary of the geometry, that is supposed to be free of traction). If not, I got this error when I try to run in parallel:
.......

Processor 43
Number of cells = 36102
Number of faces shared with processor 42 = 8398
Number of faces shared with processor 44 = 6951
Number of processor patches = 2
Number of processor faces = 15349
Number of boundary faces = 22640

Processor 44
Number of cells = 36144
Number of faces shared with processor 43 = 6951
Number of processor patches = 1
Number of processor faces = 6951
Number of boundary faces = 30615

Number of processor faces = 544514
Max number of processor patches = 2
Max number of faces between processors = 32095

Wrote decomposition to "/home/constant/cellDecomposition" for use in manual decomposition.

Wrote decomposition as volScalarField to cellDist for use in post-processing.


--> FOAM FATAL IO ERROR:

Cannot find 'value' entry on patch oldInternalFaces of field D in file "/0/D"
which is required to set the values of the generic patch field.
(Actual type solidTractionFree)

Please add the 'value' entry to the write function of the user-defined boundary-condition
or link the boundary-condition into libfoamUtil.so

file: /0/D::boundaryField:ldInternalFaces from line 214 to line 214.

From function genericFvPatchField<Type>::genericFvPatchField(con st fvPatch&, const Field<Type>&, const dictionary&)
in file fields/fvPatchFields/basic/generic/genericFvPatchField.C at line 72.

FOAM exiting


Why do I have to set a value, considering that a run withg 1 core is not requesting it? Could this be the reason of getting different results?






4. On top of that, when I try to change (in the serial run) the fixedValue Bcs on interfaces to contact BCs, I got the the following warning for a long time period and then solver crashes: first iteration is done, and on Time=2 solver starts to complain about ggi (and the log weight is several GB)


--> FOAM Warning :
From function newGGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file lnInclude/newGGIInterpolationWeights.C at line 463
polygonIntersection is returning a zero surface area between
Master face: 2332 and Neighbour face: 2334 intersection area = 0
Please check the two quick-check algorithms for newGGIInterpolation. Something is missing.
--> FOAM Warning :
From function newGGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file lnInclude/newGGIInterpolationWeights.C at line 463
polygonIntersection is returning a zero surface area between
Master face: 2334 and Neighbour face: 2332 intersection area = 0


I googled this and no info was obtained. Any ideas on how to deal with this warning? It is making the simulation explode...


My "D" file is:


FoamFile
{
format ascii;
class volVectorField;
location "0";
object D;
version 2.0;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 1 0 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{

support
{
type fixedDisplacement;
value uniform (0 0 0);
}
load
{
type solidTraction;
traction uniform ( 0 0 0 );
pressure uniform 10000;
value uniform (0 0 0);
}

oldInternalFaces
{



type solidTractionFree;
value uniform (0 0 0);

/*traction uniform ( 0 0 0 );
pressure uniform 0;
value uniform (0 0 0);*/
}



//////////////////////////////////////////////////////////////////////////////////////////////////
screw_to_hex
{
type solidContact;
master yes;
rigidMaster no;
shadowPatch hex_to_screw;
normalContactModel standardPenalty;
standardPenaltyNormalModelDict
{
relaxationFactor 0.01;
penaltyScale 1;
}
//frictionContactModel frictionless;
frictionContactModel standardPenalty;
standardPenaltyFrictionModelDict
{
relaxationFactor 0.01;
penaltyScale 1;
frictionLaw coulomb;
frictionLawDict
{
frictionCoeff 0.1;
}
}
value uniform (0 0 0);
}

hex_to_screw
{
type solidContact;
shadowPatch screw_to_hex;
master no;
value uniform (0 0 0);

}

///////////////////////////////////////////////////////////////////////////////////////////////



screw_to_abutment
{
type solidContact;
master yes;
rigidMaster no;
shadowPatch abutment_to_screw;
normalContactModel standardPenalty;
standardPenaltyNormalModelDict
{
relaxationFactor 0.01;
penaltyScale 1;
}
//frictionContactModel frictionless;
frictionContactModel standardPenalty;
standardPenaltyFrictionModelDict
{
relaxationFactor 0.01;
penaltyScale 1;
frictionLaw coulomb;
frictionLawDict
{
frictionCoeff 0.1;
}
}
value uniform (0 0 0);


}

abutment_to_screw
{
type solidContact;
shadowPatch screw_to_abutment;
master no;
value uniform (0 0 0);

}




hex_to_abutment
{
type solidContact;
master yes;
rigidMaster no;
shadowPatch abutment_to_hex;
normalContactModel standardPenalty;
standardPenaltyNormalModelDict
{
relaxationFactor 0.01;
penaltyScale 1;
}
//frictionContactModel frictionless;
frictionContactModel standardPenalty;
standardPenaltyFrictionModelDict
{
relaxationFactor 0.01;
penaltyScale 1;
frictionLaw coulomb;
frictionLawDict
{
frictionCoeff 0.1;
}
}
value uniform (0 0 0);


}

abutment_to_hex
{
type solidContact;
shadowPatch hex_to_abutment;
master no;
value uniform (0 0 0);

}

}



If files or case is needed for testing, just let me know. Any hints are really appreciated.


Best

[CRI_CFD]

BTW, trying to mimic the procedure with elasticSolidFoam (foam extend 4.1), no matter if I define solid contacts or fixedValue Bcs, I got the following error:


--> FOAM FATAL ERROR:

request for regIOobject solidInterface from objectRegistry region0 failed
available objects of type regIOobject are
0()

From function objectRegistry::lookupObject<Type>(const word&) const
in file /home/foam-extend-4.1/src/foam/lnInclude/objectRegistryTemplates.C at line 198.

FOAM aborting



Any suggestions?

[CRI_CFD]

Also I am attaching the result of checkMesh -allTopology -allGeometry (never used before):



Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
all points: 2476007
live points: 2476007
all faces: 5724647
live faces: 5724647
internal faces: 4773130
cells: 1624632
boundary patches: 9
point zones: 0
face zones: 9
cell zones: 3

Overall number of cells of each type:
hexahedra: 1413014
prisms: 4
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 211614

Checking topology...
Boundary definition OK.
Point usage OK.
Upper triangular ordering OK.
Topological cell zip-up check OK.
Face vertices OK.
Face-face connectivity OK.
<<Writing 15 cells with with single non-boundary face to set oneInternalFaceCells
<<Writing 2693 cells with with single non-boundary face to set twoInternalFacesCells
*Number of regions: 3
The mesh has multiple regions which are not connected by any face.
<<Writing region information to "0/cellToRegion"
Nuumber of cells per region:
0 481487
1 238127
2 905018


Checking patch topology for multiply connected surfaces ...
Patch Faces Points Area [m^2] Surface topology Bounding box
load 13868 14641 3.12562 ok (non-closed singly connected) (-1.5911 -1.59413 8.99977) (1.60861 1.60489 9.02198)
oldInternalFaces 750272 753190 238.778 multiply connected (shared edge) (-2.24111 -2.24595 -4.55002) (2.26071 2.25467 9.02198)
abutment_to_screw 11002 11837 2.29938 ok (non-closed singly connected) (-1.64789 -1.67049 3.94957) (1.68761 1.68823 3.96079)
abutment_to_hex 6638 7355 1.53075 ok (non-closed singly connected) (-1.41761 -1.43288 3.95345) (1.44458 1.45025 3.96615)
hex_to_screw 73102 73676 26.1481 multiply connected (shared edge) (-1.41761 -1.43288 1.14249) (1.44458 1.45025 3.95769)
hex_to_abutment 6638 7355 1.53075 ok (non-closed singly connected) (-1.41761 -1.43288 3.95345) (1.44458 1.45025 3.96615)
support 5893 6101 1.33194 ok (non-closed singly connected) (-0.689553 -0.634327 -4.55512) (0.686041 0.664033 -4.49743)
screw_to_hex 73102 73676 26.1481 multiply connected (shared edge) (-1.41761 -1.43288 1.14249) (1.44458 1.45025 3.95769)
screw_to_abutment 11002 11837 2.29938 ok (non-closed singly connected) (-1.64789 -1.67049 3.94957) (1.68761 1.68823 3.96079)
<<Writing 52 conflicting points to set nonManifoldPoints

Checking geometry...
This is a 3-D mesh
Overall domain bounding box (-2.24111 -2.24595 -4.55512) (2.26071 2.25467 9.02198)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Mesh (non-empty, non-wedge) dimensions 3
Boundary openness (1.143e-15 1.41986e-17 1.34325e-15) Threshold = 1e-06 OK.
Max cell openness = 4.21797e-16 OK.
Max aspect ratio = 11.3371 OK.
Minumum face area = 7.01425e-06. Maximum face area = 0.0106761. Face area magnitudes OK.
Min volume = 4.12618e-07. Max volume = 0.000672995. Total volume = 74.6727. Cell volumes OK.
Mesh non-orthogonality Max: 50.2305 average: 16.8299 Threshold = 70
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 2.04439 OK.
Min/max edge length = 0.000144859 0.159832 OK.
Writing 183438 near (closer than 1.49953e-05 apart) points to set nearPoints
*There are 11336 faces with concave angles between consecutive edges. Max concave angle = 74.5451 degrees.
Writing 11336 faces with concave angles to set concaveFaces
Face flatness (1 = flat, 0 = butterfly) : average = 0.998331 min = 0.599533
*There are 40 faces with ratio between projected and actual area < 0.8
Minimum ratio (minimum flatness, maximum warpage) = 0.599533
Writing 40 warped faces to set warpedFaces
Cell determinant (wellposedness) : minimum: 0 average: 11.6112
***Cells with small determinant found, number of cells: 22034
Writing 22034 under-determined cells to set underdeterminedCells

Failed 1 mesh checks.

End

[bigphil]

Originally Posted by CRI_CFD

1. I tested a simplified geometry with 3 zones and 3 contact regions. I set all interfaces to a fixedValue BC (value (0 0 0 )). If I go for a serial ran, I get the log attached. Everything seems to converge properly, so...does this mean the mesh is fine? I would say yes, but...

Yes, the mesh is probably good enough.

Originally Posted by CRI_CFD

2. If I go for a parallel run, with exactly the same settings I got the following log (different and no convergence is reached)

Hmnnn, that's strange. That means something is not correct for parallel. I would suggest trying it with solids4foam (or apologies, maybe this is solids4foam?).

Originally Posted by CRI_CFD

3. To be honest, the step 2 is not exactly the same case, as I needed to specify a value for my solidTraction BC on oldInternal faces (actually the outer boundary of the geometry, that is supposed to be free of traction). If not, I got this error when I try to run in parallel:

Why do I have to set a value, considering that a run withg 1 core is not requesting it? Could this be the reason of getting different results?

For boundary conditions that decomposePar is unaware of, it expects them all to have values specified. Using a value of "uniform (0 0 0)" should be fine.

Alternatively, the libraries defining these boundary conditions must be loaded in controlDict (not needed here).

Originally Posted by CRI_CFD

4. On top of that, when I try to change (in the serial run) the fixedValue Bcs on interfaces to contact BCs, I got the the following warning for a long time period and then solver crashes: first iteration is done, and on Time=2 solver starts to complain about ggi (and the log weight is several GB)


--> FOAM Warning :
From function newGGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file lnInclude/newGGIInterpolationWeights.C at line 463
polygonIntersection is returning a zero surface area between
Master face: 2332 and Neighbour face: 2334 intersection area = 0
Please check the two quick-check algorithms for newGGIInterpolation. Something is missing.
--> FOAM Warning :
From function newGGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file lnInclude/newGGIInterpolationWeights.C at line 463
polygonIntersection is returning a zero surface area between
Master face: 2334 and Neighbour face: 2332 intersection area = 0


I googled this and no info was obtained. Any ideas on how to deal with this warning? It is making the simulation explode...

I suggest following the approach of "make things simpler until it works"; from that, it should be possible to figure out the problem and possibly its cause.
For example, use fricitonless instead of friction. In addition, make sure each separate region has a fixedDisplacement condition to stop it floating away. Then try a low value for penaltyScale (this scales the stiffness of the contact springs), e.g. try 1e-6. You could even try 0.0, which should give the same behaviour as solidTractionFree.

[CRI_CFD]

Hi again,


Thanks for the answer.


Based on your suggestions, I have created a very simple case (two cubes in contact). A load is applied on top of the smaller cube ("tool") and zero displacement is fixed on the bottom face of the bigger cube ("base").


I am trying to run solids4foam. I followed the procedure mentioned, got my mesh, remove background mesh, splitMeshRegions...


When I ran solids4foam, I got the same warning and then simulation stops:


--> FOAM Warning :
From function newGGIInterpolation<MasterPatch, SlavePatch>::calcAddressing()
in file lnInclude/newGGIInterpolationWeights.C at line 463
polygonIntersection is returning a zero surface area between
Master face: 399 and Neighbour face: 288 intersection area = 0
Please check the two quick-check algorithms for newGGIInterpolation. Something is missing.
Floating point exception (core dumped)



I am totally stuck on this, I do not understand where is the problems with this simplified geometry. I tried the fixed Bcs and solver runs, but when I introduce contacts (no matter if I use frictionless or others, or if I reduce penaltyScale) newGGIInterpolation is there again.


The example is uploaded here:


https://drive.google.com/file/d/1Syb...ew?usp=sharing


Any hints will be really appreciated.



PD: BTW, just in case it help others when meshing several volumes, I realised that meshing different volumes in contact was not working in foam extend. I bet the reason for this is that I was not using closed stl files. Now, once they have been closed, the zoneInside=true does the job)

[CRI_CFD]

Still stuck on this...


I wonder if the reason of this warning/error can be:


1. Concave cells: is GGI able to deal with concave cells created by snappy?
2. mergeMeshes does not merge the faces between the two set of cells (is stitchMesh also needed? Bit I guess if I use stitchMesh I will loose my interface, so the faces will become a unique faceZone, so no contacts are allowed, becoming an internal surface)
3. The use of splitMeshRegions with cellZones option (i am using now cellZonesOnly)


Also...if I want to establish contacts, I am not able to have fixedValue BC on both regions (only the one with the support will be fixed to (0 0 0) displacement). If I fix the contact regions to avoid regions "floating away", then no tractions will be transferred from the slave surface to the master surface.