[charbel]

Hi,

I make pisoFsiFoam to solve fluid structure interaction in turbulent flow basing on icoFsiFoam. In fact I only modified the fluid solver.
The pisoFsiFoam and a 2D test case can be obtained from the following link:

http://dl.free.fr/azoxbwmku

After compiling this solver I tried it for the simple 2D case. But the solver stops after the first time step without showing error message. In fact the message I have from the log.pisoFsiFoam file is below:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Build : 1.6-ext
Exec : pisoFsiFoam
Case : /home/charbel.habchi/OpenFOAM/charbel.habchi-1.6-ext/run/vorflex/pisoFsiFoam2DTestCase/fluid
nProcs : 1
SigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).

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

Create dynamic mesh for time = 0

Selecting dynamicFvMesh dynamicMotionSolverFvMesh
// using new solver syntax:
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0;
}

// using new solver syntax:
U
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0;
}

Selecting motion solver: laplaceFaceDecomposition
// using new solver syntax:
motionU
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0;
}

Selecting motion diffusivity: quadratic
// using new solver syntax:
U
{
solver PCG;
preconditioner DIC;
tolerance 1e-09;
relTol 0.01;
}

Reading transportProperties

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RASModel
Selecting RAS turbulence model realizableKE
realizableKECoeffs
{
Cmu 0.09;
A0 4;
C2 1.9;
sigmak 1;
sigmaEps 1.2;
}

Reading field U

Reading mechanical properties

Normalising E : E/rho

Calculating Lame's coefficients

Plane Stress

mu = 9.60735e+06 Pa
lambda = 9.46396e+06 Pa
Reading coupling properties

Starting time loop

Courant Number mean: 0 max: 0.000205494 velocity magnitude: 0.286
deltaT = 1.19999e-06
Time = 2.19999e-06

Setting pressure
Total pressure force = (0 0 0)
DICPCG: Solving for Ux, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for motionUx, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for motionUy, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for motionUx, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG: Solving for motionUy, Initial residual = 0, Final residual = 0, No Iterations 0
volume continuity errors : volume = 0.000236435, max error = 0, sum local = 0, global = 0
Motion magnitude: mean = 0 max = 0
DILUPBiCG: Solving for Ux, Initial residual = 1, Final residual = 1.05423e-09, No Iterations 1
DILUPBiCG: Solving for Uy, Initial residual = 1, Final residual = 3.21083e-10, No Iterations 1
DICPCG: Solving for p, Initial residual = 1, Final residual = 8.24415e-07, No Iterations 254
DICPCG: Solving for p, Initial residual = 0.028181, Final residual = 9.13188e-07, No Iterations 230
Moving mesh time step continuity errors : sum local = 2.3574e-10, global = -4.60078e-13, cumulative = -4.60078e-13
DICPCG: Solving for p, Initial residual = 0.00804306, Final residual = 8.90858e-07, No Iterations 223
DICPCG: Solving for p, Initial residual = 0.00473995, Final residual = 9.3184e-07, No Iterations 216
Moving mesh time step continuity errors : sum local = 2.64889e-10, global = -1.6892e-12, cumulative = -2.14928e-12
DILUPBiCG: Solving for epsilon, Initial residual = 0.000179671, Final residual = 2.3196e-10, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 1, Final residual = 1.23968e-08, No Iterations 2
ExecutionTime = 1.01 s


End

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

As I am new user of OpenFoam I do not know where is the problem.
I thank you in advance for your help.

Best Regards,
Charbel.

[RugbyGandalf]

Dear Charbel,

i do not know about this solver. Normally, if calculation crashes it will write out an error message.
Did you already check the endTime entry in controlDict (can be found in the system directory)? With entries in the controlDict, you can specify where to start and to end the calculation and also specify the timeStep.

[charbel]

Dear RugbyGandalf,
Thank you.

In fact, I do not think that the controlDict poses problem as you can see it below:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application pisoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 20;
deltaT 0.000001;
writeControl adjustableRunTime;
writeInterval 0.02;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression compressed;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
adjustTimeStep yes;
maxCo 0.5;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

I think the problem is more in one of the attached files that you can also find in http://dl.free.fr/azoxbwmku
Maybe I made some mistake while trying to built the solver!

Best regards,
Charbel.

[RugbyGandalf]

Dear Charbel,

you're right, that wont cause any problems.
I was just overflying the problem and that was my first idea...

Maybe someone knows about the mistakes in your files :?

I am sorry, but at the time I am hardly in stress, cause of writing exams and a thesis

Good luck!

[charbel]

Hi,

I found some errors in solveFluid.H that are now corrected in the attached file below.

The simulations run now but crashes after some time steps showing large courant number. I will try to see the problem.

I have a question remain, in solveFluid.H what most I use:
U -= fvc::grad(p)/UEqn.A() or U -= rUA*fvc::grad(p)

Can you help please!!

Best regards,
Charbel.

[T.D.]

Dear Charbel, kifak mnee7,

look here
Ai Ui = Ri--------> solve(UEqn == -fvc::grad(p))
aP-------------------> AU=UEqn().A()
U = (1/aP) H[U ]-----> U=UEqn().H()/AU
φ = S · H[U ]----------> phi=vc::interpolate(U)& mesh.Sf()
∇ · ( a1 ∇p) = ∇ · (φ) -----> fvm::laplacian(1.0/AU, p)==fvc::div(phi)
U = U − (1/aP) ∇p ---------> U-=fvc::grad(p)/AU


and we know also volScalarField rUA = 1.0/UEqn.A();
so fvc::grad(p)/UEqn.A() is a vector (grad(p)) times scalar (1/UEqn.A()) where the latter is equal to rUA
so vector x scalar = scalar x vector (commutative)
so the same. even if it was in matrix form.

Best Regards,

T.D.

[charbel]

Hello Talib,
Thank you

[T.D.]

you are welcome