tweaked parameters to be more realistic -> coffee cup cools down early (as expected)

mgrid_2/getmatrix.cpp

@@ -434,7 +434,8 @@ double FEM_Matrix::ThermalConductivity(const int subdomain)
     {
         // ceramic mug
         case 0:
-            lambda = 10.0;           // anything from 1 to 4
+            //lambda = 5.0;
+            lambda = 30.0;
             break;
 
         // water
@@ -673,16 +674,16 @@ double FEM_Matrix::Heat_transfer_coefficient(const int subdomain)
 {
     double alpha = 0.0;
 
-    switch (subdomain)
+    switch (subdomain) // radiation!
     {
         // ceramic to air
         case 0:
-            alpha = 1.0;
+            alpha = 15.0;
             break;
 
         // air to air (Free Convection Gas - Free Convection Gas : U = 1 - 2 W/m2K  (typical window, room to outside air through glass))
         case 2:
-            alpha = 15.0;
+            alpha = 30.0;
             break;
 
         default:

mgrid_2/main.cpp

@@ -32,7 +32,7 @@ int main(int argc, char **argv )
     // ##########################################
 
     double dt = 1.0;           // time step
-    int steps = 100;            // number of time iterations
+    int max_it = 600;            // number of max iterations
 
     double u0_mug = 18.0;
     double u0_fluid = 80.0;
@@ -88,9 +88,11 @@ int main(int argc, char **argv )
     double goal_temp = 67.0;
     double goal_perc = 60.0;
 
+    double prev_temp = u0_mug;
+
     double time_count = 0;
-    //while (average_cup_temperature < goal_temp)
-    while (percentage_temp_reached < goal_perc)
+    while (average_cup_temperature < goal_temp && prev_temp <= average_cup_temperature)
+    //while (percentage_temp_reached < goal_perc)
     //for (int step = 0; step < steps; ++step)
     {
         vector<double> G(Mdt.Nrows(), 0.0);        
@@ -105,6 +107,8 @@ int main(int argc, char **argv )
         JacobiSolve(SK, H, uv);          // solve: (M/dt + K + C) * u_{n+1}  =  F + M/dt * u_{n}
                                          //        ----- SK -----               ------ H -------
 
+        prev_temp = average_cup_temperature;
+
         average_cup_temperature = mesh_c.AverageVectorFunction_perSubdomain(uv, 0);
         percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 0, goal_temp);
         cout << "Average cup temperature: " << average_cup_temperature << " after " << time_count << " seconds. " << endl;
@@ -143,8 +147,9 @@ int main(int argc, char **argv )
     percentage_temp_reached = mesh_c.CheckTemp_mult(uv, 1, goal_temp);
 
     time_count = 0;
-    //while (average_coffee_temperature > goal_temp)
-    while (percentage_temp_reached > goal_perc)
+    int step = 0;
+    while (average_coffee_temperature > goal_temp && step < max_it)
+    //while (percentage_temp_reached > goal_perc)
     {
         vector<double> G(Mdt.Nrows(), 0.0);        
         Mdt.Mult(G, uv);                                           // G = M/dt * u_{n}
@@ -163,6 +168,7 @@ int main(int argc, char **argv )
         cout << "Average coffee temperature: " << average_coffee_temperature << " after " << time_count << " seconds. " << endl;
         cout << "% of elements above temperature " << goal_temp << "ºC: " << percentage_temp_reached << endl;
         time_count += dt;
+        ++step;
     }
     t4 = system_clock::now();  // stop timer
 

solid-cpp/ownSolver.cpp

@@ -42,7 +42,7 @@ int main(int argc, char **argv )
     // ##########################################
 
     double dt = 1.0;           // time step
-    int steps = 100;            // number of time iterations
+    int max_it = 300;            // number of time iterations
     double u_out = 18.0;
 
     // ##########################################
@@ -157,8 +157,10 @@ int main(int argc, char **argv )
     goal_perc = 60.0;
 
     double time_count = 0;
-    //while (average_coffee_temperature > goal_temp)
-    while (percentage_temp_reached > goal_perc)
+    int step = 0;
+    while (average_coffee_temperature > goal_temp && step < max_it)
+    //while (percentage_temp_reached > goal_perc)
+    for(int i = 0; i < 100; ++i)
     {
         preciceDt = participantSolid.getMaxTimeStepSize();
         solverDt = 1.0;
@@ -199,7 +201,7 @@ int main(int argc, char **argv )
         for (int i = 0; i < numberOfVertices; ++i)
         {
             int nodeIndex = wetNodes[i];
-            heatFlux[i] = 15.0*(uv[nodeIndex] - 31.0);  // alpha is assumed 15.0 regardless of material
+            heatFlux[i] = 30.0*(uv[nodeIndex] - 31.0);  // alpha is assumed for air material at the top
         } 
         participantSolid.writeData("Solid-Mesh", "Heat-Flux", vertexIDs, heatFlux);
 
@@ -209,6 +211,7 @@ int main(int argc, char **argv )
 
 
         time_count += dt;
+        ++step;
     }
     auto t4 = system_clock::now();  // stop timer
     auto duration = duration_cast<microseconds>(t4 - t3);        // duration in microseconds
@@ -220,6 +223,7 @@ int main(int argc, char **argv )
     participantSolid.finalize();
 
 
+    mesh_c.Visualize(uv);
 
     return 0;
 }