CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Error with laplacian of two volScalarFields (http://www.cfd-online.com/Forums/openfoam-programming-development/94982-error-laplacian-two-volscalarfields.html)

Hisham December 1, 2011 13:22

Error with laplacian of two volScalarFields
 
Dear Foamers,

I have this equation:
Code:

const volScalarField TwoMuPlusLambda = 2 * solidShearMod + solidLambdaMod;

//solve
fvVectorMatrix gEqn
      (
      fvm::d2dt2(solidRho, initial_D)
      ==
        fvm::laplacian( TwoMuPlusLambda, initial_D, "laplacian(DD,D)")
        + divSigmaExp
      );

solidRho and TwoMuPlusLambda are both volScalarFields .... The code compiles but it gives infinity for initial_D .... If I use a scalar instead of TwoMuPlusLambda, it gives rather logical output.
It seems there is something wrong with laplacian for a volScalarField Gamma ... Can anyone give me feedback on that.
I use OF 2.0.x on Ubuntu 11.04

Best regards,
Hisham

Hisham December 2, 2011 10:57

I tried to compare two equations one with a scalar and another with a volScalarField for comparisons and they yield different results for the same numerical input ... Results for the scalar are the logical ones:

Code:


volVectorField D
(
 IOobject
 (
  "D",
  runTime.timeName(),
  mesh,
  IOobject::MUST_READ,
  IOobject::AUTO_WRITE
 ),
 mesh
);
Info<< "\nReading field initial_D\n" << endl;

// Initial displacement
volVectorField initial_D
(
 IOobject
 (
  "initial_D",
  runTime.timeName(),
  mesh,
  IOobject::MUST_READ,
  IOobject::AUTO_WRITE
 ),
 mesh
);

volSymmTensorField sigmaD
(
    IOobject
    (
        "sigmaD",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    solidShearMod * twoSymm(fvc::grad(D)) + solidLambdaMod * (I*tr(fvc::grad(D)))
);

volSymmTensorField initialSigmaD
(
    IOobject
    (
        "initialSigmaD",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    sigmaD
);

volVectorField divSigmaExp
(
    IOobject
    (
        "divSigmaExp",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    fvc::div(sigmaD)
);

Info<< "Initializing field of explicit plastic div(sigma) divSigmaExpPlastic\n" << endl;
volVectorField divSigmaExpPlastic
(
    IOobject
    (
        "divSigmaExpPlastic",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    fvc::div(sigmaD)
);

volTensorField gradInitialD = fvc::grad(initial_D);
volTensorField gradD = fvc::grad(D);

Info << "gradD  " << gradD  << endl;
Info << "gradInitialD  " << gradInitialD  << endl;

initialSigmaD =  solidShearMod * twoSymm(gradInitialD) + solidLambdaMod * (I * tr(gradInitialD));
sigmaD = mu * twoSymm(gradD) + lambda * (I * tr(gradD)) ;
Info << "sigmaD  " << sigmaD << endl;
Info << "initialSigmaD  " << initialSigmaD << endl;
divSigmaExp =  fvc::div
              (
                  initialSigmaD
                  - (2 * solidShearMod + solidLambdaMod ) * gradInitialD,
                  "div(sigmaD)"
                )
                + solidRho * Gravity;
Info << "divSigmaExp    " <<  divSigmaExp << endl;
divSigmaExpPlastic =        fvc::div
              (
                  sigmaD
                  - (2 * mu + lambda ) * gradD,
                  "div(sigmaD)"
                )
                + rho * Gravity;                                                    Info << "divSigmaExpPlastic  " << divSigmaExpPlastic << endl;

                                                    // + solidRho * Gravity ;

const volScalarField TwoMuPlusLambda = 2 * solidShearMod + solidLambdaMod;

//solve
fvVectorMatrix gEqn
      (
      fvm::d2dt2(solidRho, initial_D)
      ==
      fvm::laplacian( 2* solidShearMod  + solidLambdaMod, initial_D, "laplacian(DD,D)")
      + divSigmaExp
      );
fvVectorMatrix eEqn
      (
      fvm::d2dt2(rho, D)
        ==
      fvm::laplacian( 2* mu  + lambda, D, "laplacian(DD,D)")
      + divSigmaExpPlastic
      );
initialResidual = eEqn.solve().initialResidual();
initialResidual = gEqn.solve().initialResidual();

Info << "D  " << D << endl;
Info << "initial_D    " << initial_D  << endl;

I applied the solver on the geometry of the plate with a hole tutorial with a fixed value for all patches of a zero vector value.

The results of "D" are logical but for "initial_D" they are not .... the numbers are the same until the gEq and eEq are first solved.

Hisham December 5, 2011 06:11

To summarize:

When I do a laplacian of a volVectorField multiplied by a constant scalar it runs OK and values are logical

But when I replace the constant scalar with a constant volScalarField with a uniform value of the same value as the scalar in the former case, it compiles but gives very large values for the volVectorField solved for that paraFoam views as "inf" or very high values when I change the values.

I need a volScalarField to introduce nonuniform values in practice.

Can someone please help me? It seems like a bug, but I want to make sure I'm not doing it wrong before sending a bug report.

Hisham December 5, 2011 06:46

I found out what the problem was:

The values of the volScalarField at patches needed to be updated to be equal to the values of cells with faces at patch. I tried that and it gave same results as a the scalar. I do not know if this fully solved the problem or not .... Until other bugs appear, I would assume it debugged :rolleyes:

hrushi.397 February 12, 2013 08:04

Quote:

Originally Posted by Hisham (Post 334661)
I found out what the problem was:

The values of the volScalarField at patches needed to be updated to be equal to the values of cells with faces at patch. I tried that and it gave same results as a the scalar. I do not know if this fully solved the problem or not .... Until other bugs appear, I would assume it debugged :rolleyes:

How exactly did you do the update?

Hisham February 12, 2013 08:46

Quote:

Originally Posted by hrushi.397 (Post 407383)
How exactly did you do the update?

Hello

First you do:
Code:

forAll(variable.boundaryField(), patchi)
{
      variable.boundaryField()[patchi] == variable.boundaryField()[patchi].patchInternalField();
}

Then:
Code:

D.correctBoundaryConditions();
Hisham

hrushi.397 February 12, 2013 10:07

Thanks Hisham :). Appreciate your help.

ChrisA February 20, 2013 17:48

Thank you very much for following up with the solution, I never would have figured out my issue otherwise. (I didn't even think of this as an issue)

hrushi.397 August 12, 2013 03:07

Hi Hisham,

I am facing a somewhat similar issue now. It has to do with laplacian. I am trying to play around with laplacianFoam and I tried to change the equation from implicit to explicit(fvm::laplacian(DT,T) to fvc::laplacian(DT,T)) just to see if I get a bit different answer. To my surprise the laplacianFoam crashed. I checked the laplacian values and it seems laplacian explodes and runs into very high values.

What could have gone wrong? Is this a bug in Openfoam?

Thanks,

Hrushi

bigphil August 12, 2013 06:27

Quote:

Originally Posted by hrushi.397 (Post 445084)
Hi Hisham,

I am facing a somewhat similar issue now. It has to do with laplacian. I am trying to play around with laplacianFoam and I tried to change the equation from implicit to explicit(fvm::laplacian(DT,T) to fvc::laplacian(DT,T)) just to see if I get a bit different answer. To my surprise the laplacianFoam crashed. I checked the laplacian values and it seems laplacian explodes and runs into very high values.

What could have gone wrong? Is this a bug in Openfoam?

Thanks,

Hrushi

Hi Hrushi,

Just to clarify, you changed the equation from implicit:
Code:

                fvm::ddt(T) - fvm::laplacian(DT, T)
to explicit like this:
Code:

                fvm::ddt(T) - fvc::laplacian(DT, T)
?

If you did this then, this method will be conditionally stable (as with all explicit methods).
Therefore the results will 'blow up' unless the stability criterion is satisfied. For 1-D heat conduction, the time-step is limited:
dt < rho*c*dx*dx/(2*k)
where dt is the time-step, rho is the density, c is the specific heat, dx is the mesh spacing and k is the conductivity.

Philip

hrushi.397 August 12, 2013 06:50

Thanks for the reply, Phil. I was thinking on similar lines today. I shall definitely try this too.


All times are GMT -4. The time now is 21:33.