applications/solvers/solidMechanics/elasticNonLinTLSolidFoam/writeFields.H

   1 if (runTime.outputTime())
   2 {
   3     volScalarField epsilonEq
   4     (
   5         IOobject
   6         (
   7            "epsilonEq",
   8            runTime.timeName(),
   9            mesh,
  10            IOobject::NO_READ,
  11            IOobject::AUTO_WRITE
  12         ),
  13         sqrt((2.0/3.0)*magSqr(dev(epsilon)))
  14     );
  15
  16     Info<< "Max epsilonEq = " << max(epsilonEq).value()
  17         << endl;
  18
  19     volScalarField sigmaEq
  20     (
  21         IOobject
  22         (
  23             "sigmaEq",
  24             runTime.timeName(),
  25             mesh,
  26             IOobject::NO_READ,
  27             IOobject::AUTO_WRITE
  28         ),
  29         sqrt((3.0/2.0)*magSqr(dev(sigma)))
  30     );
  31
  32     Info<< "Max sigmaEq = " << max(sigmaEq).value()
  33         << endl;
  34
  35     //- Calculate Cauchy stress
  36     volTensorField F = I + gradU;
  37     volScalarField J = det(F);
  38
  39     //- update density
  40     rho = rho/J;
  41
  42     volSymmTensorField sigmaCauchy
  43     (
  44         IOobject
  45         (
  46            "sigmaCauchy",
  47             runTime.timeName(),
  48             mesh,
  49             IOobject::NO_READ,
  50             IOobject::AUTO_WRITE
  51         ),
  52         (1/J) * symm(F.T() & sigma & F)
  53     );
  54
  55     //- Cauchy von Mises stress
  56     volScalarField sigmaCauchyEq
  57     (
  58         IOobject
  59         (
  60             "sigmaCauchyEq",
  61             runTime.timeName(),
  62             mesh,
  63             IOobject::NO_READ,
  64             IOobject::AUTO_WRITE
  65         ),
  66         sqrt((3.0/2.0)*magSqr(dev(sigmaCauchy)))
  67     );
  68
  69     Info<< "Max sigmaCauchyEq = " << max(sigmaCauchyEq).value()
  70         << endl;
  71
  72
  73     volTensorField Finv = inv(F);
  74     volSymmTensorField epsilonAlmansi
  75     (
  76         IOobject
  77         (
  78             "epsilonAlmansi",
  79             runTime.timeName(),
  80             mesh,
  81             IOobject::NO_READ,
  82             IOobject::AUTO_WRITE
  83         ),
  84         symm(Finv & epsilon & Finv.T())
  85     );
  86
  87 //     volVectorField traction
  88 //       (
  89 //        IOobject
  90 //        (
  91 //     "traction",
  92 //     runTime.timeName(),
  93 //     mesh,
  94 //     IOobject::NO_READ,
  95 //     IOobject::AUTO_WRITE
  96 //     ),
  97 //        mesh,
  98 //        dimensionedVector("zero", dimForce/dimArea, vector::zero),
  99 //        calculatedFvPatchVectorField::typeName
 100 //        );
 101 //     forAll(traction.boundaryField(), patchi)
 102 //       {
 103 //     const tensorField& Fbinv = Finv.boundaryField()[patchi];
 104 //     vectorField nCurrent = Fbinv & n.boundaryField()[patchi];
 105 //     traction.boundaryField()[patchi] =
 106 //       nCurrent & sigmaCauchy.boundaryField()[patchi];
 107 //       }
 108
 109 //     //- write boundary forces
 110 //     //- integrate (sigma2PK & F) over reference area
 111 //     //- which is equivalent to integrating sigmaCauchy
 112 //     //- over the deformed area
 113 //     Info << nl;
 114 //     forAll(mesh.boundary(), patchi)
 115 //       {
 116 //     Info << "Patch " << mesh.boundary()[patchi].name() << endl;
 117 //     const tensorField& Fb = F.boundaryField()[patchi];
 118 //     vectorField totalForce = mesh.Sf().boundaryField()[patchi] & (sigma.boundaryField()[patchi] & Fb);
 119 //     //vectorField totalForce2 = Sf.boundaryField()[patchi] & (sigmaCauchy.boundaryField()[patchi]);
 120
 121 //     vector force = sum( totalForce );
 122 //     //vector force2 = sum( totalForce2 );
 123 //     Info << "\ttotal force is " << force << " N" << endl;
 124 //     //Info << "\ttotal force2 is " << force2 << " N" << endl;
 125
 126 //     const tensorField& Fbinv = Finv.boundaryField()[patchi];
 127 //     vectorField nCurrent = Fbinv & n.boundaryField()[patchi];
 128 //     nCurrent /= mag(nCurrent);
 129 //     scalar normalForce = sum( nCurrent & totalForce );
 130 //     Info << "\tnormal force is " << normalForce << " N" << endl;
 131 //     scalar shearForce = mag(sum( (I - sqr(nCurrent)) & totalForce ));
 132 //     Info << "\tshear force is " << shearForce << " N" << endl;
 133
 134     //if(mesh.boundary()[patchi].name() == "right")
 135     //{
 136     //const vectorField& nOrig = n.boundaryField()[patchi];
 137     //Info << "\tNormal force on right is " << (nCurrent & totalForce) << nl << endl;
 138     //Info << "\tShear force on right is " << ((I - sqr(nCurrent)) & totalForce) << nl << endl;
 139     //Info << "\tpatch gradient is " << U.boundaryField()[patchi].snGrad() << endl;
 140     //Info << "\tpatch gradient (norm) is " << (nCurrent & U.boundaryField()[patchi].snGrad()) << endl;
 141     //Info << "\tpatch gradient (shear) is " << ((I - sqr(nCurrent)) & U.boundaryField()[patchi].snGrad()) << endl;
 142     //Info << "\tpatch Almansi (normal) is " << (nCurrent & (nCurrent & epsilonAlmansi.boundaryField()[patchi])) << endl;
 143     //Info << "\tpatch Almansi (shear) is " << ( (I - sqr(nCurrent)) & (nCurrent & epsilonAlmansi.boundaryField()[patchi])) << endl;
 144     //Info << "\tpatch Green (normal) is " << (nOrig & (nOrig & epsilon.boundaryField()[patchi])) << endl;
 145     //Info << "\tpatch Green (shear) is " << ( (I - sqr(nOrig)) & (nOrig & epsilon.boundaryField()[patchi])) << endl;
 146     //Info << "\tpatch Cauchy stress (normal) is " << (nCurrent & (nCurrent & sigmaCauchy.boundaryField()[patchi])) << endl;
 147     //}
 148
 149 //     if(mesh.boundary()[patchi].type() != "empty")
 150 //       {
 151 //             vector Sf0 = Sf.boundaryField()[patchi][0];
 152 //             symmTensor sigma0 = sigmaCauchy.boundaryField()[patchi][0];
 153 //             Info << "sigmab[0] is " << sigma0 << nl
 154 //                  << "Sfb is " << Sf0 << nl
 155 //                  << "force is " << (Sf.boundaryField()[patchi][0]&sigma.boundaryField()[patchi][0]) << nl
 156 //          << "Sfx*sigmaxx " << (Sf0[vector::X]*sigma0[symmTensor::XX]) <<nl
 157 //                  << "Sfy*sigmaxy " << (Sf0[vector::Y]*sigma0[symmTensor::XY]) << nl
 158 //                  << "Sfx*sigmayx " << (Sf0[vector::X]*sigma0[symmTensor::XY]) << nl
 159 //                  << "Sfy*sigmayy " << (Sf0[vector::Y]*sigma0[symmTensor::YY]) << nl
 160 //                  << endl;
 161 //       }
 162 //    Info << endl;
 163 //    }
 164
 165     runTime.write();
 166 }