nLight Beam Class

applications/solvers/additiveFoam/additiveFoam.C

@@ -42,6 +42,8 @@ Description
 #include "movingHeatSourceModel.H"
 #include "foamToExaCA/foamToExaCA.H"
 
+#include "interface/interfacePoints.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -67,12 +69,16 @@ int main(int argc, char *argv[])
 
  movingHeatSourceModel sources(mesh);
 
+ DynamicList<List<scalar>> interfaceData;
+
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
  Info<< "\nStarting time loop\n" << endl;
 
  while (runTime.run())
  {
+ volScalarField R("R", mag(fvc::ddt(T)));
+
  #include "updateProperties.H"
 
  #include "readTimeControls.H"
@@ -94,6 +100,8 @@ int main(int argc, char *argv[])
  }
 
  #include "thermo/TEqn.H"
+
+ #include "interface/interface.H"
 
  ExaCA.update();
 
@@ -104,6 +112,8 @@ int main(int argc, char *argv[])
  << nl << endl;
  }
 
+ #include "interface/interfaceWrite.H"
+
  ExaCA.write();
 
  return 0;

applications/solvers/additiveFoam/interface/interface.H

@@ -0,0 +1,28 @@
+// write out interface dimensions
+Info<< "interface span " << "(" << Tliq.value() << "): "
+ << interfacePoints(T, Tliq.value()) << endl;
+
+Info<< "interface span " << "(" << Tsol.value() << "): "
+ << interfacePoints(T, Tsol.value()) << endl;
+
+// write out the maximum temperature
+Info<< "maximum temperature: " << gMax(T.internalField()) << endl;
+
+volScalarField G("G", mag(fvc::grad(T)));
+
+forAll(mesh.cells(), celli)
+{
+ if ( (T[celli] <= Tliq.value()) && (T.oldTime()[celli] > Tliq.value()) )
+ {
+ const point pt = mesh.C()[celli];
+
+ interfaceData.append(
+ {
+ pt[0],
+ pt[1],
+ pt[2],
+ R[celli] / G[celli],
+ G[celli]
+ });
+ }
+}

applications/solvers/additiveFoam/interface/interfacePoints.H

@@ -0,0 +1,81 @@
+vector interfacePoints(const volScalarField& field, const scalar& iso)
+{
+ const fvMesh& mesh = field.mesh();
+
+ // set local reference to mesh data
+ const labelUList& owner = mesh.owner();
+ const labelUList& neighbour = mesh.neighbour();
+
+ const volVectorField& cc = mesh.C();
+
+ DynamicList<vector> positions(mesh.nFaces());
+
+ // check internal faces
+ for(label facei=0; facei < mesh.nInternalFaces(); facei++)
+ { 
+ const label own = owner[facei];
+ const label nei = neighbour[facei];
+
+ scalar minFace = min(field[own], field[nei]);
+ scalar maxFace = max(field[own], field[nei]);
+
+ if ((minFace < iso) && (maxFace >= iso))
+ {
+ vector d = cc[nei] - cc[own];
+ vector p = cc[own] + d*(iso - field[own])/(field[nei] - field[own]);
+ positions.append(p);
+ }
+ }
+
+ // check boundary faces
+ const volScalarField::Boundary& fieldBf = field.boundaryField();
+
+ forAll(fieldBf, patchi)
+ {
+ const fvPatchScalarField& fieldPf = fieldBf[patchi];
+
+ const labelUList& faceCells = fieldPf.patch().faceCells();
+
+ if (fieldPf.coupled())
+ {
+ // processor boundary : interpolate across face
+ const vectorField ccn(cc.boundaryField()[patchi].patchNeighbourField());
+ const scalarField fn(fieldPf.patchNeighbourField());
+
+ forAll(faceCells, facei)
+ {
+ label own = faceCells[facei];
+
+ scalar minFace = min(field[own], fn[facei]);
+ scalar maxFace = max(field[own], fn[facei]);
+
+ if ((minFace < iso) && (maxFace >= iso))
+ {
+ vector d = ccn[facei] - cc[own];
+ vector p = cc[own] + d*(iso - field[own])/(fn[facei] - field[own]);
+ positions.append(p);
+ }
+ }
+ }
+ else
+ {
+ // physical boundary : take face centre in liquid
+ const vectorField fc(mesh.boundaryMesh()[patchi].faceCentres());
+ const scalarField pif(fieldPf.patchInternalField());
+
+ forAll(fc, i)
+ {
+ if (pif[i] >= iso)
+ {
+ positions.append(fc[i]);
+ }
+ }
+ }
+ }
+
+ positions.shrink();
+
+ boundBox isoBb(positions);
+
+ return returnReduce(isoBb.span(), maxOp<vector>());
+}

applications/solvers/additiveFoam/interface/interfaceWrite.H

@@ -0,0 +1,22 @@
+const fileName interfacePath
+(
+ runTime.rootPath()/runTime.globalCaseName()/"solidificationData"
+);
+
+mkDir(interfacePath);
+
+OFstream os
+(
+ interfacePath + "/" + "data_" + Foam::name(Pstream::myProcNo()) + ".csv"
+);
+
+for(int i=0; i < interfaceData.size(); i++)
+{
+ int n = interfaceData[i].size()-1;
+
+ for(int j=0; j < n; j++)
+ {
+ os << interfaceData[i][j] << ",";
+ }
+ os << interfaceData[i][n] << endl;
+}

applications/solvers/additiveFoam/movingHeatSource/Make/files

@@ -11,6 +11,7 @@ heatSourceModels/heatSourceModel/heatSourceModel.C
 heatSourceModels/heatSourceModel/heatSourceModelNew.C
 heatSourceModels/superGaussian/superGaussian.C
 heatSourceModels/modifiedSuperGaussian/modifiedSuperGaussian.C
+heatSourceModels/nLight/nLight.C
 
 movingHeatSourceModel/movingHeatSourceModel.C
 

applications/solvers/additiveFoam/movingHeatSource/heatSourceModels/heatSourceModel/heatSourceModel.C

@@ -239,15 +239,15 @@ Foam::heatSourceModel::qDot()
  );
  volScalarField& qDot_ = tqDot.ref();
 
- // sample gaussian distribution at desired resolution
+ // sample distribution at desired resolution
  const scalar power_ = movingBeam_->power();
 
  if (power_ > small)
  {
  const vector position_ = movingBeam_->position();
 
  // udpate the absorbed power and heat source normalization term
- const scalar aspectRatio = 
+ const scalar aspectRatio =
  dimensions_.z() / min(dimensions_.x(), dimensions_.y());
 
  dimensionedScalar absorbedPower
@@ -275,7 +275,7 @@ Foam::heatSourceModel::qDot()
  );
 
  const pointField& points = mesh_.points();
- 
+
  treeBoundBox beamBb
  (
  position_ - 1.5*dimensions_,