[qi.yang@polimi.it]

Dear all,


I am using solids4Foam to simulate wave-induced seabed response. Firstly I imported the wave pressure obtained from interfoam results on the seabed surface by using "timeaveragedvaryingmappefield" and then I run the case. However, I don't think the pore water pressure is the right one. What is your opinion? since I am not pretty sure that the numerical settings, especially the boundary condition that I set are correct. I attached all the files here. Any suggestions will be highly appreciated.

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

dimensions      [0 1 0 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            solidSymmetry;
        patchType       symmetryPlane;
        value           uniform (0 0 0);
    }

    outlet
    {
        type            solidSymmetry;
        patchType       symmetryPlane;
        value           uniform (0 0 0);
    }

    sides
    {
        type            solidSymmetry;//fixedDisplacementZeroShear;
        patchType	symmetryPlane;
	value           uniform (0 0 0);
    }

   cylinder// soilStructureInterface
    {
        type            fixedDisplacement;
        value           uniform (0 0 0);
    }

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

    top
    {
        type            solidTraction;
        traction        uniform (0 0 0);
        pressure        uniform 0;
        value           uniform (0 0 0);
    }

}

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "22";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            symmetryPlane;
    }

    outlet
    {
        type            symmetryPlane;
    }

    sides
    {
        type            symmetryPlane;//zeroGradient;
    }

    cylinder//soilStructureInterface
    {
        type            zeroGradient;
    }

    ground
    {
        type            zeroGradient;
    }

    top
    {
      //  type            timefixedValue;
       // value           uniform 0;
    type timeVaryingMappedFixedValue;
        setAverage  off;

	}

}

Thanks,

[qi.yang@polimi.it]

Originally Posted by qi.yang@polimi.it

Dear all,


I am using solids4Foam to simulate wave-induced seabed response. Firstly I imported the wave pressure obtained from interfoam results on the seabed surface by using "timeaveragedvaryingmappefield" and then I run the case. However, I don't think the pore water pressure is the right one. What is your opinion? since I am not pretty sure that the numerical settings, especially the boundary condition that I set are correct. I attached all the files here. Any suggestions will be highly appreciated.

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

dimensions      [0 1 0 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            solidSymmetry;
        patchType       symmetryPlane;
        value           uniform (0 0 0);
    }

    outlet
    {
        type            solidSymmetry;
        patchType       symmetryPlane;
        value           uniform (0 0 0);
    }

    sides
    {
        type            solidSymmetry;//fixedDisplacementZeroShear;
        patchType	symmetryPlane;
	value           uniform (0 0 0);
    }

   cylinder// soilStructureInterface
    {
        type            fixedDisplacement;
        value           uniform (0 0 0);
    }

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

    top
    {
        type            solidTraction;
        traction        uniform (0 0 0);
        pressure        uniform 0;
        value           uniform (0 0 0);
    }

}

FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "22";
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            symmetryPlane;
    }

    outlet
    {
        type            symmetryPlane;
    }

    sides
    {
        type            symmetryPlane;//zeroGradient;
    }

    cylinder//soilStructureInterface
    {
        type            zeroGradient;
    }

    ground
    {
        type            zeroGradient;
    }

    top
    {
      //  type            timefixedValue;
       // value           uniform 0;
    type timeVaryingMappedFixedValue;
        setAverage  off;

	}

}

Thanks,

In terms of the boundary condition of D, should I also apply the Traction boundary condition (Which is the dynamic wave pressure) on the top of seabed, in addition to set the dynamic pressure for the boundary condition of pore pressure p at surface? I saw the tutorial in "minigeotechfoam", solidtraction on the surface is zero w.r.t. dU and only wave pressure boundary is prescribed on surface regarding P. Anyone can give me some tips? Thanks.

[bigphil]

The solidTraction condition will enforce the total traction (as defined by the total stress sigma), where total traction is the sum of the effective traction (soil skeleton) and pore pressure.

You can see how sigma is calculated in solids4foam-release/src/solids4foamModels/materialModels/mechanicalModel/mechanicalLaws/linearGeometryLaws/poroLinearElastic/poroLinearElastic.C in the correct function:

void Foam::poroLinearElastic::correct(volSymmTensorField& sigma)
{
    // Calculate effective stress                                                                                                                                                                            
    linearElastic::correct(sigma);

    // Lookup the pressure field from the solver                                                                                                                                                             
    const volScalarField& p = mesh().lookupObject<volScalarField>("p");

    // Calculate the total stress as the sum of the effective stress and the                                                                                                                                 
    // pore-pressure                                                                                                                                                                                         
    sigma -= (p + p0_)*symmTensor(I);
}

What do you want it to do?

Philip

[qi.yang@polimi.it]

Originally Posted by bigphil

The solidTraction condition will enforce the total traction (as defined by the total stress sigma), where total traction is the sum of the effective traction (soil skeleton) and pore pressure.

You can see how sigma is calculated in solids4foam-release/src/solids4foamModels/materialModels/mechanicalModel/mechanicalLaws/linearGeometryLaws/poroLinearElastic/poroLinearElastic.C in the correct function:

void Foam::poroLinearElastic::correct(volSymmTensorField& sigma)
{
    // Calculate effective stress                                                                                                                                                                            
    linearElastic::correct(sigma);

    // Lookup the pressure field from the solver                                                                                                                                                             
    const volScalarField& p = mesh().lookupObject<volScalarField>("p");

    // Calculate the total stress as the sum of the effective stress and the                                                                                                                                 
    // pore-pressure                                                                                                                                                                                         
    sigma -= (p + p0_)*symmTensor(I);
}

What do you want it to do?

Philip

Hi, Prof. Thanks a lot for your reply. Currently, I obtained the dynamic pressure of the wave and I would like to set it as a boundary condition at the top of the seabed with respect to the pore pressure. Afterward, I will use this solver to get the pore pressure distribution and displacement of soil. Like Tang Tian's tutorial did "minigeotechfoam/tutorials/soil/twoSurfacePlastBiotFoam/waveOverPlateOnSeabed" in which he applied wave boundary condition and also the variable is du instead of u used in solids4Foam. Therefore, I am a little confused about how could set the boundary conditions correctly in my case? I tried one case in which the pressure at top of the seabed for p was prescribed but nothing for D. Thanks for your attention.

[bigphil]

The way you are setting the dynamic pressure condition seems fine, but what condition should be applied to the soil skeleton? Should this effective traction be zero? Or should the total traction be zero? Where total traction is the effective traction plus the pore pressure contribution.
Currently the solids4foam solidTraction enforces a total traction condition but it would be straight-forward to addd a flag which removes the pore pressure contribution so it becomes an effective traction condition.

I have not checked through the minigeotechfoam waveOverPlateOnSeabed case in detail so I am not sure which is appropriate in this case.

[qi.yang@polimi.it]

Originally Posted by bigphil

The way you are setting the dynamic pressure condition seems fine, but what condition should be applied to the soil skeleton? Should this effective traction be zero? Or should the total traction be zero? Where total traction is the effective traction plus the pore pressure contribution.
Currently the solids4foam solidTraction enforces a total traction condition but it would be straight-forward to addd a flag which removes the pore pressure contribution so it becomes an effective traction condition.

I have not checked through the minigeotechfoam waveOverPlateOnSeabed case in detail so I am not sure which is appropriate in this case.

Hi, prof. This is what I am still concerning. According to my experience, only dynamic wave pressure should be given in terms of boundary conditions, which I did in COMSOL a long time ago where Biot's consolidation equation can be solved by specifying zero displacements at top of the seabed. Here, in solids4foam, it seems to me that also solid traction at top should be given, which is equal to the value of dynamic pressure but if so, the sigma is zero. How do you think?

[bigphil]

In COMSOL, I assume they also solve two equations: pore pressure equation, and displacement equation. In that case, they also need two boundary conditions.
However, if as a user you only specify one condition then COMSOL is assuming the other condition for you in a (hopefully) smart manner.

In your case, I guess that the total traction on the seabed is the dynamic pressure you have (ignoring viscous forces from the sea). In that case, the total traction is your dynamic pressure, which you specify as the pore pressure; then the solidTraction condition (which specifies the total traction) needs to also specify this same pressure. In that way, the effective traction on the soil skeleton is zero at the seabed.

However, I suggest you find a simple representative case where you know the correct answer (e.g. you could set up a case in COMSOL) and then verify you get the same answers in solids foam, and you may need to try a variety of boundary conditions combinations until you are confident about the correct usage.

Philip

[qi.yang@polimi.it]

Originally Posted by bigphil

In COMSOL, I assume they also solve two equations: pore pressure equation, and displacement equation. In that case, they also need two boundary conditions.
However, if as a user you only specify one condition then COMSOL is assuming the other condition for you in a (hopefully) smart manner.

In your case, I guess that the total traction on the seabed is the dynamic pressure you have (ignoring viscous forces from the sea). In that case, the total traction is your dynamic pressure, which you specify as the pore pressure; then the solidTraction condition (which specifies the total traction) needs to also specify this same pressure. In that way, the effective traction on the soil skeleton is zero at the seabed.

However, I suggest you find a simple representative case where you know the correct answer (e.g. you could set up a case in COMSOL) and then verify you get the same answers in solids foam, and you may need to try a variety of boundary conditions combinations until you are confident about the correct usage.

Philip

Thanks a lot. I understood the procedure however, I found that setting boundary condition of porepressure at seabed surface can be achieved by using " timeVaryingMappedFixedValue" whereas for the solidtraction, how can I specify a pressure field dependent on positions and time?

[qi.yang@polimi.it]

I rechecked that the effective stress at seabed surface should be zero and pore pressure is specified as dynamic wave pressure. I am thinking that if I apply fixedDisplacement for D at seabed surface (but in reality, there is no limitation on displacement) and dynamic wave pressure for p at the top, the total solid traction sigma should be equal to the wave pressure. So that there is no need to use solidtraction BC? I obtained the porepressure is the same as the condition of solidtraction used. Otherwise, it is possible to set the solidtraction equal to p?

[bigphil]

Originally Posted by qi.yang@polimi.it

I rechecked that the effective stress at seabed surface should be zero and pore pressure is specified as dynamic wave pressure. I am thinking that if I apply fixedDisplacement for D at seabed surface (but in reality, there is no limitation on displacement) and dynamic wave pressure for p at the top, the total solid traction sigma should be equal to the wave pressure. So that there is no need to use solidtraction BC? I obtained the porepressure is the same as the condition of solidtraction used. Otherwise, it is possible to set the solidtraction equal to p?

If you set the displacement to fixedDisplacement "(0 0 0)" then the effective traction will probably not be equal to zero, as the boundary will be able to support an external load.

The simplest thing will be to allow solidTraction to equal p; I will add an option to solidTraction called "setEffectiveTraction" which will default to "off". But you can then set it to "on" and it will then look-up the p field and force the total traction to p. Do you have any simple/quick test case that we can check it on?

Philip

[bigphil]

I have added the option to commit aed78b8f. It is currently on the development branch but will be pushed to the master in the coming days.

It seems to work correctly e.g. if I set "setEffectiveTraction" to "on" and set the traction/pressure in the solidTraction condition to zero, then I can see that the normal component of the stress tensor at the boundary is equal to the dynamic pressure. I tried it on the poroelasticity/stripFooting tutorial.

[qi.yang@polimi.it]

Originally Posted by bigphil

I have added the option to commit aed78b8f. It is currently on the development branch but will be pushed to the master in the coming days.

It seems to work correctly e.g. if I set "setEffectiveTraction" to "on" and set the traction/pressure in the solidTraction condition to zero, then I can see that the normal component of the stress tensor at the boundary is equal to the dynamic pressure. I tried it on the poroelasticity/stripFooting tutorial.

Thanks a lot for your contribution, Prof. So, I have to recompile the solver, right? After the new option was updated successfully, I will try my case again.

[bigphil]

Yes, you should change to the development branch, pull the latest changes, and then compile the code; to do this, run the following commands from the main solids4foam-release directory:

$> git checkout development
$> git pull
$> ./Allwmake

[qi.yang@polimi.it]

Hi Prof,

I found I made a mistake before that I took the total pressure into account. In fact, the dynamic pressure should be equal to the total pressure output from interfoam subtracted by the hydrostatic pressure (rho*g*h). After I used correct wave dynamic pressure. Both two options ( setEffectiveTraction on / off) returned similar pore pressure, velocity and, displacement distribution (reasonable), but the quantity of displacement and velocity is different while the pore pressure is the same (which makes me confused). I will validate additional simple cases later. Thank you so much.

[qi.yang@polimi.it]

Dear Prof.,

I compared the results with literature, in terms of a simple case, namely seabed response under a wave. As you found also, the effective normal stress now is equal to the dynamic wave pressure. However, the pore pressure and vertical effective stress don't agree with each other (the maximum and minimum values of pore pressure are basically consistent due to the right value of dynamic wave pressure I think). Additionally, I found the displacement seems not converged. I attached my "testcase" here. Hope you can help me have a look at it if possible. Thanks.

[bigphil]

There is no mesh in the attached case. Can you please check this?

Also, the "sigma" field is the total stress, not the effective stress.

You need to subtract the pressure term from the sigma field to get the effective stress.
You can create the ZZ effective stress in ParaView using the "Calculator" field (Filters -> Alphabetical -> Calculator). You can then create a point or cell field as "sigma + p".

Philip

[qi.yang@polimi.it]

Originally Posted by bigphil

There is no mesh in the attached case. Can you please check this?

Also, the "sigma" field is the total stress, not the effective stress.

You need to subtract the pressure term from the sigma field to get the effective stress.
You can create the ZZ effective stress in ParaView using the "Calculator" field (Filters -> Alphabetical -> Calculator). You can then create a point or cell field as "sigma + p".

Philip

Hi Prof,

I misunderstood that option you added (if turned on), which indicates only effective stress applied on the boundary, instead of sigma being the effective stress in the field, isn't it? As you can see, the effective normal stress (In fact, I applied sigma_ZZ-p) gets similar to the result found in the literature.
However, the quantities deviate to some extent. Due to the limitation of file size, I transferred the mesh file to you via Wetransfer. Sorry for that. Of course, I can create a blockMeshDict file for this simple case if you need it.
Thanks for your attention.

[bigphil]

Originally Posted by qi.yang@polimi.it

I misunderstood that option you added (if turned on), which indicates only effective stress applied on the boundary, instead of sigma being the effective stress in the field, isn't it?

Yes you are correct.

The sigma tensor always represents the total stress, regardless of this option being on or off.

This option for the solidTraction boundary condition just allows the user to enforce a total traction condition or an effective traction condition.

For convenience, I have just added the field "sigmaEff" (effective stress) to the base mechanical law so the pro-elastic mechanical laws will now write this and you can view it in ParaView. See commit 21081a3e on the development branch; this will be merged to the master branch in the coming week.

As to why your case is not the exact same as your reference from literature, I suggest you double check the geometry, materials, loading conditions, mesh sensitivity, and time-step sensitivity.

Philip

[qi.yang@polimi.it]

Originally Posted by bigphil

Yes you are correct.

The sigma tensor always represents the total stress, regardless of this option being on or off.

This option for the solidTraction boundary condition just allows the user to enforce a total traction condition or an effective traction condition.

For convenience, I have just added the field "sigmaEff" (effective stress) to the base mechanical law so the pro-elastic mechanical laws will now write this and you can view it in ParaView. See commit 21081a3e on the development branch; this will be merged to the master branch in the coming week.

As to why your case is not the exact same as your reference from literature, I suggest you double check the geometry, materials, loading conditions, mesh sensitivity, and time-step sensitivity.

Philip

Thanks a lot. I will check it as soon as possible.

[qi.yang@polimi.it]

I tried to use a smaller time step, which is 0.1 (previously 0.5), with the same mesh (it is fine enough, I think). However, it seems to me D is not converged at each time step. After I adopted a coarse mesh, D had a smaller residual and the solution is similar. I also attached the displacement of soil calculated by COMSOL, with general partial differential equations, which I expect the equations used in the two codes are consistent. You can simply see the difference in the horizontal displacement. What is surprising to me is that the displacement on the surface is especially strange. I am considering the boundary conditions whether I specified are correct or not. What is your idea?
After I checked the source code, I found some differences. e.g.

    rKprime_
    (
        (saturation_/KWater_)
      + (1.0 - saturation_)
        /dimensionedScalar("atmosphericPressure", dimPressure, 1e+05)
    )

the formulation in literature is

    rKprime_
    (
        (1/KWater_)
      + (1.0 - saturation_)
        /(dimensionedScalar("atmosphericPressure", dimPressure, 1e+05)+gamma_w*g*h)
    )