[Hisham]

Dear Foamers,

I have this equation:

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]

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:

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]

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]

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

[hrushi.397]

Originally Posted by Hisham

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

How exactly did you do the update?

[Hisham]

Originally Posted by hrushi.397

How exactly did you do the update?

Hello

First you do:

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

Then:

D.correctBoundaryConditions();

Hisham

[hrushi.397]

Thanks Hisham . Appreciate your help.

[ChrisA]

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]

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]

Originally Posted by hrushi.397

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:

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

to explicit like this:

                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]

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