[MFWilliams]

Hi, another thread from me...


So I think I have a better understanding of what is wrong with my simulation and I would like peoples opinion on it and any suggestions to improve it.

Summary of the problem: I am trying to create an FSI simulation for blood flowing through a brain aneurysm. I have an inlet velocity waveform and specify a pressure of 0 Pa at the outlets. The other BC's are below.


What is happening is that I am getting a wave of fluid velocity increase as it moves down the inlet tube, as in the first image. For all of the initial iterations of the FSI sim it converges well (in about 20 FSI corrector iterations) and doesnt have a problem. Then between 0.00025 and 0.0005 suddenly it takes longer to converge in each FSI iteration and eventually after maxing out at 200 FSI corrector iterations it crashes. As it gets slower to converge, it starts to get a lot slower after printing the output:
Time = 0.000XXXX, iteration: X
After it prints that line there is a long pause, which is longer than it takes to run all of the other calculations in that iteration. Could this indicate what is going wrong?


I have also tried to set the coupling time to start at 0.01 seconds. This did roughly the same thing where it was fine with no coupling, then the coupling started and it happily converged for about 0.0002 seconds, then all of a sudden it starts taking longer to converge and finally crashes.


The FSI parameters that I am using are:
FSI relaxation factor = 0.001
FSI tolerance = 1e-06
reuse = 4
Solid density = 1300
Poissions ratio = 0.45
Youngs Modulus = 1.6e6
deltaT = 2.5e-6
maxDeltaT = 1e-5
maxCo = 0.2
Solid domain relaxation = 0.6
(No relaxation for pressure or velocity in fluid domain)


Also what would be the effect of ading relaxation for pressure or velocity in the fluid domain, might this help or produce more problems?


All of these sims have used a meshes with non-orthogonalities no greater than 60 and no more than 80000 elements. Some have used boundary layers some havent. The checkMesh for the fluid and solid from image attached is below:
Fluid mesh:

Time = 0

Mesh stats
    all points:           16283
    live points:          16283
    all faces:            44045
    live faces:           44045
    internal faces:       39371
    cells:                13939
    boundary patches:     4
    point zones:          0
    face zones:           0
    cell zones:           0

Overall number of cells of each type:
    hexahedra:     13765
    prisms:        44
    wedges:        0
    pyramids:      86
    tet wedges:    0
    tetrahedra:    44
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch               Faces    Points   Area [m^2]  Surface topology
    inlet               106      113      8.04803e-06 ok (non-closed singly connected)
    fluidInterface      4514     4539     0.000310381 ok (non-closed singly connected)
    outlet2             38       42       2.36964e-06 ok (non-closed singly connected)
    outlet1             16       18       8.41541e-07 ok (non-closed singly connected)

Checking geometry...
    This is a 3-D mesh
    Overall domain bounding box (-0.0524018 -0.0690049 -0.550233) (-0.0426186 -0.0581131 -0.533035)
    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.37959e-16 -1.15997e-16 9.4579e-17) Threshold = 1e-06 OK.
    Max cell openness = 2.79196e-16 OK.
    Max aspect ratio = 5.45218 OK.
    Minumum face area = 5.24083e-09. Maximum face area = 1.64157e-07.  Face area magnitudes OK.
    Min volume = 2.1492e-13. Max volume = 4.08655e-11.  Total volume = 3.15458e-07.  Cell volumes OK.
    Mesh non-orthogonality Max: 38.8545 average: 5.18479 Threshold = 70
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 0.544832 OK.

Mesh OK.

Solid mesh:

    From function void polyMesh::initMesh()
    in file meshes/polyMesh/polyMeshInitMesh.C at line 81
    Truncating neighbour list at 22506 for backward compatibility
Time = 0

Mesh stats
    all points:           13617
    live points:          13617
    all faces:            31650
    live faces:           31650
    internal faces:       22506
    cells:                9028
    boundary patches:     3
    point zones:          0
    face zones:           0
    cell zones:           0

Overall number of cells of each type:
    hexahedra:     9016
    prisms:        12
    wedges:        0
    pyramids:      0
    tet wedges:    0
    tetrahedra:    0
    polyhedra:     0

Checking topology...
    Boundary definition OK.
    Point usage OK.
    Upper triangular ordering OK.
    Face vertices OK.
    Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces ...
    Patch               Faces    Points   Area [m^2]  Surface topology
    solidEnds           116      174      6.866e-06   ok (non-closed singly connected)
    fluidInterfaceSolid 4514     4539     0.000310381 ok (non-closed singly connected)
    outerSolidFace      4514     4539     0.000360442 ok (non-closed singly connected)

Checking geometry...
    This is a 3-D mesh
    Overall domain bounding box (-0.0527121 -0.0693202 -0.550553) (-0.0422981 -0.0577922 -0.532975)
    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.14638e-16 3.75844e-16 -2.4529e-16) Threshold = 1e-06 OK.
    Max cell openness = 2.0773e-16 OK.
    Max aspect ratio = 4.89239 OK.
    Minumum face area = 1.33733e-08. Maximum face area = 2.47861e-07.  Face area magnitudes OK.
    Min volume = 2.31903e-12. Max volume = 4.31077e-11.  Total volume = 1.07163e-07.  Cell volumes OK.
    Mesh non-orthogonality Max: 55.353 average: 8.25151 Threshold = 70
    Non-orthogonality check OK.
    Face pyramids OK.
    Max skewness = 1.24647 OK.

Mesh OK.

Boundary conditions:
U:

internalField   uniform (0 0 0);

boundaryField
{
    outlet2
    {
        type            fluxCorrectedVelocity;
    }
    inlet
    {
        type            groovyBC;
        refValue        uniform (0 0 0);
        refGradient     uniform (0 0 0);
        valueFraction   uniform 1;
        value           uniform (0 0 0);
        valueExpression "(a1*sin(b1*x+c1) + a2*sin(b2*x+c2) + a3*sin(b3*x+c3) + a4*sin(b4*x+c4) + a5*sin(b5*x+c5) + a6*>        gradientExpression "vector(0,0,0)";
        fractionExpression "1";
        evaluateDuringConstruction 0;
        cyclicSlave     0;
        variables
25
(
"x=time();"
"a1=1.476;"
"b1=1.011;"
"c1=0.4674;"
"a2=1.206;"
"b2=1.132;"
"c2=3.477;"
"a3=0.05853;"
"b3=8.393;"
"c3=-0.2345;"
"a4=0.04648;"
"b4=16.76;"
"c4=-0.6503;"
"a5=0.00276;"
"b5=5.108;"
"c5=2.317;"
"a6=0.02639;"
"b6=33.52;"
"c6=-1.617;"
"a7=0.0228;"
"b7=25.14;"
"c7=-1.522;"
"a8=0.01827;"
"b8=41.87;"
"c8=-2.746;"
)
;
        timelines       (
);
        lookuptables    (
);
    }
    outlet1
    {
        type            fluxCorrectedVelocity;
    }
    fluidInterface
    {
        type            newMovingWallVelocity;
        value           uniform (0 0 0);
    }
}

P:

boundaryField
{
    outlet1
    {
        type        fixedValue;
        value       uniform 0;
    }
    outlet2
    {
        type        fixedValue;
        value       uniform 0;
    }
    inlet
    {
        type       zeroGradient;
    }
    fluidInterface
    {
        type       extrapolatedPressure;
        value      uniform 0;
    }

D:

internalField   uniform (0 0 0);

boundaryField
{
    solidEnds
    {
        type            fixedDisplacement;
        value           uniform (0 0 0);
    }
    fluidInterfaceSolid
    {
        type            solidTraction;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        value           uniform (0 0 0);
    }
    outerSolidFace
    {
        type            solidTraction;
        traction        uniform ( 0 0 0 );
        pressure        uniform 0;
        value           uniform (0 0 0);
    }
}

Also I have seen the BC zeroGradientCorrected in the microBeamInCrossflow tutorial, what is this and might it help?


In my very novice opinion it seems that as all of the simulations are failing after roughly the same amount of time after being coupled, there must be some sort of numerical error that starts insignificant but then takes 0.00025-0.0005 seconds to increase to an extent where it affects the solution, and the time it takes to crash the simulation is proportional to the refinement of the mesh, ie more refined takes less time to break. If this is correct how might I find out what is happening? I also think it must be to do with the solid settings or mesh as I can run a CFD sim with the same fluid parameters and it will converge perfectly.


Here is a log file for a simulation which used a slightly more refined mesh, but which crashed in a very similar fashion.

https://drive.google.com/file/d/10J8...ew?usp=sharing




EDIT: I have also just realised that there is a missing '}' at the end of the pressure BC for the fluid. Could this tiny error be a problem, becuase it is also missing on the 3dTube tutorial case which I built my fsi case from.


EDIT2: here is a link to the most recent case file: https://drive.google.com/drive/folde...mu?usp=sharing

[Neb]

Hi,
have you solved your problem? i have the exact same block in the simulation.