initial commit for version 1.6.x patch release
[OpenFOAM-1.6.x.git] / src / engine / include / StCorr.H
bloba776fef5c3d85c636ef5b5297925ef1876fed343
1     dimensionedScalar StCorr("StCorr", dimless, 1.0);
3     if (ign.igniting())
4     {
5         // Calculate volume of ignition kernel
6         dimensionedScalar Vk("Vk", dimVolume, gSum(c*mesh.V().field()));
7         dimensionedScalar Ak("Ak", dimArea, 0.0);
9         if (Vk.value() > SMALL)
10         {
11             // Calculate kernel area from its volume
12             // and the dimensionality of the case
14             switch(mesh.nGeometricD())
15             {
16                 case 3:
17                 {
18                     // Assume it is part-spherical
19                     scalar sphereFraction
20                     (
21                         readScalar
22                         (
23                             combustionProperties.lookup
24                             (
25                                 "ignitionSphereFraction"
26                             )
27                         )
28                     );
30                     Ak = sphereFraction*4.0*mathematicalConstant::pi
31                        *pow
32                         (
33                             3.0*Vk
34                            /(sphereFraction*4.0*mathematicalConstant::pi),
35                             2.0/3.0
36                         );
37                 }
38                 break;
39             
40                 case 2:
41                 {
42                     // Assume it is part-circular
43                     dimensionedScalar thickness
44                     (
45                         combustionProperties.lookup("ignitionThickness")
46                     );
48                     scalar circleFraction
49                     (
50                         readScalar
51                         (
52                             combustionProperties.lookup
53                             (
54                                 "ignitionCircleFraction"
55                             )
56                         )
57                     );
59                     Ak = circleFraction*mathematicalConstant::pi*thickness
60                        *sqrt
61                         (
62                             4.0*Vk
63                            /(circleFraction*thickness*mathematicalConstant::pi)
64                         );
65                 }
66                 break;
68                 case 1:
69                     // Assume it is plane or two planes
70                     Ak = dimensionedScalar
71                     (
72                         combustionProperties.lookup("ignitionKernelArea")
73                     );
74                 break;
75             }
77             // Calculate kernel area from b field consistent with the
78             // discretisation of the b equation.
79             volScalarField mgb =
80                 fvc::div(nf, b, "div(phiSt,b)") - b*fvc::div(nf) + dMgb;
81             dimensionedScalar AkEst = gSum(mgb*mesh.V().field());
83             StCorr.value() = max(min((Ak/AkEst).value(), 10.0), 1.0);
85             Info<< "StCorr = " << StCorr.value() << endl;
86         }
87     }