How do you like your eggs?


Introduction to PDE modeling with FEM

Ricardo López 
Mathematical Computation group

\[RawEscape]

Why Finite Elements?

Arbitrarily shaped regions

A large class of partial differential equations (PDEs)

What Is Needed for a Finite Element Analysis

A discrete representation of a region, i.e. a mesh

A partial differential equation

Boundary conditions that link the equation with the region

For an in depth overview of the Finite Element Method take a look at: FEMDocumentation/tutorial/SolvingPDEwithFEM

Kenji Lopez: https://youtu.be/hb0Elaa6gxY?si=-hrGmWlONQhAewP1

Version 1 (getting a workflow started):

Things to keep in mind

Size of objects

Units

Load the Finite Element package:

Needs["NDSolve`FEM`"]

Set the $HistoryLength to 0:

$HistoryLength=0;

Create a simple geometry of an egg:

radius=0.025;
Ω=RegionUnion[Circle[{0,0},radius,{-π/2,π/2}],Line[{{0,radius},{0,-radius}}]]

Generate a default finite element mesh from the geometry:

mesh=ToElementMesh[Ω]

Visualize the finite element mesh:

mesh["Wireframe"]

guide/PDEModeling

Documentation for the (RefLink[HeatTransferPDEComponent, paclet:ref/HeatTransferPDEComponent]).

PDEModels/tutorial/HeatTransfer/HeatTransfer

ρ (Subscript[C, p])((∂T(t, X))/(∂t))+∇·(-k ∇T(t,X))=0

T [(QuantityTF[, "K", kelvins, "Kelvins"])], r[m], z[m], t[s].

Set up the dependent and independent model variables:

vars={T[t,r,z],t,{r,z}};

Set up the parameters:

pars=<|"RegionSymmetry"->"Axisymmetric"|>;

Set up the PDE operator:

op=HeatTransferPDEComponent[vars,pars]

Boundary conditions

Convert the initial condition value from (QuantityTF[, "°C", degrees Celsius, "DegreesCelsius"]) to (QuantityTF[, "K", kelvins, "Kelvins"]):

initT=QuantityMagnitude[Quantity[20.,"DegreesCelsius"],"Kelvins"]

Set up an initial condition for the temperature inside the egg:

ics={T[0,r,z]==initT}

Specify the boundary temperature in [(QuantityTF[, "K", kelvins, "Kelvins"])]:

bcT=QuantityMagnitude[Quantity[100.,"DegreesCelsius"],"Kelvins"]

Specify the temporal behavior of the boundary condition:

bcFunction=initT+(bcT-initT)Piecewise[{{1/2(1-Cos[10*Pi*t]),t<=1/10}},1]

Visualize the temporal behavior of the boundary condition:

Show[
Plot[bcFunction,{t,0,0.2},PlotRange->All,AxesLabel->{"t", "Boundary Temperature"}]
,Graphics[{Dashed,Red,Line[{{0,bcT},{0.2,bcT}}]}]
,PlotRange->All]

Define the boundary condition:

DirichletCondition[T[t,r,z]==bcFunction,Sqrt[r^2+z^2]==radius^2]

Visualize the point element marker:

mesh["Wireframe"["MeshElement"->"PointElements","MeshElementMarkerStyle"->Blue]]

Specify the boundary condition in [(QuantityTF[, "K", kelvins, "Kelvins"])]:

bc=HeatTemperatureCondition[ElementMarker!=1,vars,pars,<|"SurfaceTemperature"->bcFunction|>]

Mesh visualization tutorial: ElementMeshVisualization

Specify the symmetry boundary condition:

sym=HeatSymmetryValue[ElementMarker==1,vars,pars]

Neumann values:
-FEMDocumentation/tutorial/SolvingPDEwithFEM#1344839347
-ref/NeumannValue#684705718

Specify an end time for the simulation:

tEnd=1;

Solve the equation:

solution=NDSolveValue[{op==sym,bc,ics},T,{t,0,tEnd},{r,z}∈mesh]

Documentation for (RefLink[NDSolveValue, paclet:ref/NDSolveValue])

Inspect the solution:

solution

Inspect the value for the solution at t=0 at half the radius:

solution[0,(radius/2),0]

Extract the minimal and maximal values from the solution:

MinMax[solution["ValuesOnGrid"]]

Compute the offset between Kelvin and degrees Celsius:

degreeCelsiusOffset=QuantityMagnitude[Quantity[0.,"DegreesCelsius"],"Kelvins"]

minmax=MinMax[solution["ValuesOnGrid"]]
minmaxCelsius=minmax-degreeCelsiusOffset

Visualize the solution at half the simulation time:

ContourPlot[solution[(tEnd/2),r,z],{r,z}∈mesh,]

Build a helper function to create contour plots of the temperature at various times:

TemperatureContourPlot[solution_]:=Function[time,
Show[
Module[{minmax=MinMax[solution["ValuesOnGrid"]],minmaxCelsius},
minmaxCelsius=minmax-degreeCelsiusOffset;
ContourPlot[solution[time,r,z],{r,z}∈solution["ElementMesh"],PlotRange->minmax+{-0.01,0.01}minmax,ColorFunction->ColorData[{"TemperatureMap",minmax}],ColorFunctionScaling->False,PlotLegends->BarLegend[{ColorData[{"TemperatureMap",minmaxCelsius}],minmaxCelsius},LegendLabel->"°[C]"],AspectRatio->Automatic,PlotLabel->("Time t = "<>StringPadRight[ToString[N[(time/60)]],5])<>"[min]"]
]
,solution["ElementMesh"]["Edgeframe"]
]
]

Define a Density version of the Helper function.

TemperatureDensityPlot[solution_]:=Function[time,
Show[
Module[{minmax=MinMax[solution["ValuesOnGrid"]],minmaxCelsius},
minmaxCelsius=minmax-degreeCelsiusOffset;
DensityPlot[solution[time,r,z],{r,z}∈solution["ElementMesh"],]
]
,solution["ElementMesh"]["Edgeframe"]
]
]

Create a number of contour plot frames :

nFrames=10;
frames=TemperatureDensityPlot[solution]/@Range[0,tEnd,tEnd/nFrames];

frames=Rasterize[#,ImageResolution->144]&/@frames;

The frames are rasterized to make the notebook more lightweight. But this optional.

Animate the frames:

ListAnimate[frames,SaveDefinitions->True]

Version 2 (an example of what can go wrong):

ρ (Subscript[C, p])((∂T(t, X))/(∂t))+∇·(-k ∇T(t,X))=0

Set up the material parameters:

pars["MassDensity"]=1000;
pars["SpecificHeatCapacity"]=3000;
pars["ThermalConductivity"]=1/2;

Inspect the parameters:

pars

Set up the simulation end time for 10 minutes:

tEnd=QuantityMagnitude[Quantity[10,"Minutes"],"Seconds"]

Regenerate the equations with the new parameters:

op=HeatTransferPDEComponent[vars,pars]

Solve the PDE:

solution=NDSolveValue[{op==sym,bc,ics},T,{t,0,tEnd},{r,z}∈mesh]

Look at the minimum and maximum values in the solution in degrees Celsius:

MinMax[solution["ValuesOnGrid"]]-degreeCelsiusOffset

Find the position of the minimal value in the solution data:

badPos=Position[solution["ValuesOnGrid"],Min[solution["ValuesOnGrid"]]]

Extract the time steps taken during the solution:

timeSteps=solution["Coordinates"][[1]]

Find the time step at which the minimal value is stored:

badPos[[1,1]]

badTime=timeSteps[[18]]

Visualize the solution at the time step with the extreme values:

Plot3D[solution[badTime,r,z],{r,z}∈mesh,]

Zoom into the solution plot and visualize the underlying finite element mesh:

Show[
ElementMeshProjection[mesh,{#[[1]],#[[2]],initT}&]["Wireframe"],
Plot3D[solution[badTime,r,z],{r,z}∈mesh,],]

Documentation for ElementMeshProjection

We need a refinement function:

ToElementMesh[Ω,MeshRefinementFunction->refinementFunction]

Create a signed distance function of the whole egg region:

sdf=SignedRegionDistance[Disk[{0,0},radius]]

Visualize the refinement function over the meshed geometry:

Plot3D[-sdf[{x,y}],{x,y}∈mesh]

Create the mesh refinement function:

refinementFunction=Function[{vertices,area},Evaluate[area>-0.0001*sdf[Mean[vertices]] ]]

Redefine the mesh:

mesh =ToElementMesh[Ω,MeshRefinementFunction->refinementFunction]
mesh["Wireframe"]

Solve the PDE and monitor the progress:

Monitor[solution=NDSolveValue[{op==sym,bc,ics},T,{t,0,tEnd},{r,z}∈mesh,EvaluationMonitor:>(monitor=Row[{"t = ",CForm[t]}])],monitor]

Documentation for Monitor

Inspect the temperature range of the newly found solution:

MinMax[solution["ValuesOnGrid"]]-degreeCelsiusOffset

Create the frames for the solution:

nFrames=10;
frames=TemperatureDensityPlot[solution]/@Range[0,tEnd,tEnd/nFrames];

Rasterize the frames:

frames=Rasterize[#,ImageResolution->144]&/@frames;

Animate the frames:

ListAnimate[frames,SaveDefinitions->True]

Inspect the egg's temperature at r=0 and z=0 after 6 minutes:

solution[6 60,0,0]-degreeCelsiusOffset

Version 3 (represent multiple material regions):

Create a multi-material geometry:

innerRadius=radius/2;
Ω=RegionUnion[Circle[{0,0},innerRadius,{-π/2,π/2}],Circle[{0,0},radius,{-π/2,π/2}],Line[{{0,radius},{0,0},{0,-radius}}]]

Specify material region markers:

materialRegions=<|"EggYolk"->1,"EggWhite"->2|>;

Create a mesh with material markers:

mesh=ToElementMesh[Ω,"RegionMarker"->{{{innerRadius/2,0},materialRegions["EggYolk"]},{{(radius+innerRadius)/2,0},materialRegions["EggWhite"]}}]

Visualize the multi-material mesh:

mesh["Wireframe"["MeshElementStyle"->FaceForm/@{Yellow,White}]]

Visualize the point element marker:

mesh["Wireframe"["MeshElement"->"PointElements","MeshElementMarkerStyle"->Blue]]

Regenerate the boundary condition:

bc=HeatTemperatureCondition[ElementMarker==2||ElementMarker==5,vars,pars,<|"SurfaceTemperature"->bcFunction|>]

Create a mesh with material markers and a mesh refinement:

mesh=ToElementMesh[Ω,"RegionMarker"->{{{innerRadius/2,0},materialRegions["EggYolk"]},{{(radius+innerRadius)/2,0},materialRegions["EggWhite"]}},"MeshRefinementFunction"->refinementFunction]

Visualize the multi-material mesh:

mesh["Wireframe"["MeshElementStyle"->FaceForm/@{Yellow,White}]]

Coimbra et al., “Density, Heat Capacity and Thermal Conductivity of Liquid Egg Products” https://doi.org/10.1016/j.jfoodeng.2005.01.043.

Specify material parameters:

pars["MassDensity"]=Piecewise[{
{1023.,ElementMarker==materialRegions["EggYolk"]},
{1143.,ElementMarker==materialRegions["EggWhite"]}
}];
pars["ThermalConductivity"]=Piecewise[{
{0.4,ElementMarker==materialRegions["EggYolk"]},
{0.6,ElementMarker==materialRegions["EggWhite"]}
}];
pars["SpecificHeatCapacity"]=Piecewise[{
{2600.,ElementMarker==materialRegions["EggYolk"]},
{3700.,ElementMarker==materialRegions["EggWhite"]}
}];

Regenerate the partial differential equation:

op=HeatTransferPDEComponent[vars,pars]

Solve the PDE:

Monitor[solution=NDSolveValue[{op==sym,bc,ics},T,{t,0,tEnd},{r,z}∈mesh,EvaluationMonitor:>(monitor=Row[{"t = ",CForm[t]}])],monitor]

Inspect the temperature range of the newly found solution:

MinMax[solution["ValuesOnGrid"]]-degreeCelsiusOffset

Create the frames for the solution:

nFrames=10;
frames=TemperatureDensityPlot[solution]/@Range[0,tEnd,tEnd/nFrames];

Rasterize the frames:

frames=Rasterize[#,ImageResolution->144]&/@frames;

Animate the frames:

ListAnimate[frames,SaveDefinitions->True]

Inspect the egg’s temperature at r=0 and z=0 after 6 minutes:

solution[6 60,0,0]-degreeCelsiusOffset

Version 4 (Refine the model further)

“Density, Heat Capacity and Thermal Conductivity of Liquid Egg Products,” https://doi.org/10.1016/j.jfoodeng.2005.01.043.

“Thermophysical and Rheological Properties of Liquid Egg White and Yolk during Thermal Pasteurization of Intact Eggs”. (WebLink[doi:10.1007/s11694-014-9183-6, https://link.springer.com/article/10.1007/s11694-014-9183-6]).

Set up measurement data for the mass density of egg yolk and white:

ρYolkData={{273.15,1037.11},{278.15,1036.39},{283.15,1035.52},{288.15,1034.57},{293.15,1033.48},{298.15,1032.27},{303.15,1030.94},{308.15,1029.47},{313.15,1027.91},{318.15,1026.21},{323.15,1024.47},{328.15,1022.54},{333.15,1020.57}};
ρWhiteData={{273.15,1043.03},{278.15,1042.87},{283.15,1042.5},{288.15,1042.},{293.15,1041.2},{298.15,1040.18},{303.15,1039.},{308.15,1037.56},{313.15,1035.97},{318.15,1034.18},{323.15,1032.13},{328.15,1029.97},{333.15,1027.54}};

Create (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) for the mass density data:

{ρYolkInterpolation,ρWhiteInterpolation}=(Interpolation[#,InterpolationOrder->1,"ExtrapolationHandler"->{Automatic,"WarningMessage"->False}]&/@{ρYolkData,ρWhiteData})

Visualize the measured data and the interpolated functions:

Show[
ListPlot[{ρYolkData,ρWhiteData}],
Plot[{ρYolkInterpolation[T],ρWhiteInterpolation[T]},{T,273.15,373.15},PlotLegends->LineLegend[{"Yolk","White"},LegendLabel->"Density"]],PlotRange->All
]

Set up measurement data for the thermal conductivity of egg yolk and white:

κYolkData={{283.15`,0.408`},{288.15`,0.41`},{293.15`,0.414`},{298.15`,0.416`},{303.15`,0.419`},{308.15`,0.421`},{313.15`,0.422`},{318.15`,0.427`},{323.15`,0.435`},{328.15`,0.44`},{333.15`,0.443`}};
κWhiteData=;

Create an (RefLink[InterpolatingFunction, paclet:ref/InterpolatingFunction]) for the thermal conductivity data:

{κYolkInterpolation,κWhiteInterpolation}=(Interpolation[#,InterpolationOrder->1,]&/@{κYolkData,κWhiteData})

Visualize the measured data and the interpolated functions:

Show[
ListPlot[{κYolkData,κWhiteData}],
Plot[{κYolkInterpolation[T],κWhiteInterpolation[T]},{T,273.15,373.15},PlotLegends->LineLegend[{"Yolk","White"},LegendLabel->"Conductivity"]],PlotRange->All
]

Set up measurement data for the heat capacity of the egg yolk and white:

CpYolkData={{273,2575.`},{275,2612.`},{278,2718.`},{281,2614.`},{285,2621.`},{291,2664.`},{295,2714.`},{298,2679.`},{301,2739.`},{307,2732.`},{311,2772.`}};

CpWhiteData={{273,3551},{275,3539},{278,3665},{281,3539},{285,3530},{291,3552},{295,3588},{298,3539},{301,3590},{307,3557},{311,3584}};

Create an LinearModelFit for the heat capacity data:

CpYolkLinearFit=LinearModelFit[CpYolkData,T,T]//Normal

CpWhiteLinearFit=LinearModelFit[CpWhiteData,T,T]//Normal

Documentation for LinearModelFit.

Visualize the measured data and the linear function that fits it:

Show[
ListPlot[{CpYolkData,CpWhiteData}],
Plot[{CpYolkLinearFit,CpWhiteLinearFit},{T,273.15,373.15},PlotLegends->LineLegend[{"Yolk","White"},LegendLabel->"(Subscript[C, p])"]],PlotRange->All
]

Set up the "MassDensity":

pars["MassDensity"]=({{\[Piecewise], ({{ρYolkInterpolation[T[t,r,z]], ElementMarker==materialRegions["EggYolk"]}, {ρWhiteInterpolation[T[t,r,z]], ElementMarker==materialRegions["EggWhite"]}})}});

Set up the "ThermalConductivity":

pars["ThermalConductivity"]=({{\[Piecewise], ({{κYolkInterpolation[T[t,r,z]], ElementMarker==materialRegions["EggYolk"]}, {κWhiteInterpolation[T[t,r,z]], ElementMarker==materialRegions["EggWhite"]}})}});

Set up the "SpecificHeatCapacity":

pars["SpecificHeatCapacity"]=({{\[Piecewise], ({{CpYolkLinearFit/.{T->T[t,r,z]}, ElementMarker==materialRegions["EggYolk"]}, {CpWhiteLinearFit/.{T->T[t,r,z]}, ElementMarker==materialRegions["EggWhite"]}})}});

Regenerate the partial differential equation:

op=HeatTransferPDEComponent[vars,pars]

Solve the PDE, monitor the progress and measure the computational time and memory it takes:

measure=AbsoluteTiming[MaxMemoryUsed[Monitor[solution=NDSolveValue[{op==sym,bc,ics},T,{t,0,tEnd},{r,z}∈mesh,EvaluationMonitor:>(monitor=Row[{"t = ",CForm[t]}])
],monitor]]/(1024.^2)];
Print["Time [s] -> ",measure[[1]],"\nMemory [MB] -> ", measure[[2]]]

Inspect the temperature range of the newly found solution:

MinMax[solution["ValuesOnGrid"]]-degreeCelsiusOffset

“Thermophysical and Rheological Properties of Liquid Egg White and Yolk during Thermal Pasteurization of Intact Eggs”. (WebLink[doi:10.1007/s11694-014-9183-6, https://link.springer.com/article/10.1007/s11694-014-9183-6]).

Set the denature temperature for egg white:

denatureTempWhite= QuantityMagnitude[(Quantity[60., "°C", degrees Celsius, "DegreesCelsius"]),"Kelvins"]

Set the denature temperature for the egg yolk:

denatureTempYolk= QuantityMagnitude[(Quantity[76., "°C", degrees Celsius, "DegreesCelsius"]), "Kelvins"]

Build a helper function to density plot the temperature distribution and visualize the denature temperatures:

TemperatureDenatureDensityPlot[solution_]:=Function[time,Show[TemperatureDensityPlot[solution][time],ContourPlot[{solution[time,r,z]==denatureTempWhite,solution[time,r,z]==denatureTempYolk},{r,z}∈solution["ElementMesh"],ContourStyle->{Directive[Green,Dashed],Directive[Magenta,Dashed]}]]]

Create the frames for the solution:

nFrames=10;
frames=TemperatureDenatureDensityPlot[solution]/@Range[0,tEnd,tEnd/nFrames];

Rasterize the frames:

frames=Rasterize[#,ImageResolution->144]&/@frames;

Animate the frames:

ListAnimate[frames,SaveDefinitions->True]

Kenji Lopez: https://youtu.be/hb0Elaa6gxY?si=-hrGmWlONQhAewP1

Inspect the egg’s temperature at r=0 and z=0 after 6 minutes:

solution[6 60,0,0]-degreeCelsiusOffset

Inspect the denaturation temperature of the yolk:

denatureTempYolk-degreeCelsiusOffset

Version 5 (Build a realistic egg geometry)

Create a helper function to compute coordinates of the shape of an egg geometry:

ClearAll[eggShapePoints]
 eggShapePoints[L_, B_,w_] := Table[Evaluate[{Re[B (Sqrt[((z (L-z))/((L^2)+4 w (L-2 z)+4 (w^2)))])/. z-> z + L/2],  z}], 
      {z, -(L/2), (L/2),(L/200)}]

“Delineating an Ovoidal Egg Shape by Length and Breadth: A Novel Two-Parametric Mathematical Model,”  https://doi.org/10.1016/j.biosystemseng.2022.11.003.

Inspect the value of the radius:

radius

Create coordinates for the egg geometry by specifying parameters:

shellCoords=eggShapePoints[2*radius,0.04,-0.007];

Create a spline:

shellSpline=BSplineCurve[shellCoords];
Graphics[{shellSpline},Axes->True]

Create an egg geometry with a subregion:

innerRadius=0.01;
Ω=RegionUnion[Circle[{0,-0.005},innerRadius,{-π/2,π/2}],shellSpline,Line[{shellCoords[[1]], shellCoords[[-1]]}]]

Create a mesh with material markers:

mesh=ToElementMesh[Ω,"RegionMarker"->{{{innerRadius/2,-0.01},materialRegions["EggYolk"]},{{(radius+innerRadius)/2,0},materialRegions["EggWhite"]}}]

Visualize the point element marker:

mesh["Wireframe"["MeshElement"->"PointElements","MeshElementMarkerStyle"->Blue]]

Inspect the boundary condition:

bc

inspect the refinementFunction

refinementFunction

Do an interpolation of the Shell’s coordinates.

rShell=Interpolation[shellCoords/. {r0_,z0_}:>{z0,r0}]

Plot the function that gives the rShell function:

Plot[rShell[z],{z,-0.025,0.025},AxesLabel->{"z","r"}]

Define a function for the radial distance to the shell and plot it:

radialDistanceToShell[{r_,z_}]:=rShell[z]-Abs[r];
Plot3D[radialDistanceToShell[{r,z}],{r,z}∈mesh]

Plot the distance to the shell cubed:

Plot3D[((radialDistanceToShell[{x,y}])^3),{x,y}∈mesh]

Define the refinement function:

refinementFunction=Function[{vertices,area},Module[{p=Mean[vertices],d},d=Evaluate@radialDistanceToShell[p];
area>(d^3) +(10^(-8))]];

Define the mesh:

mesh=ToElementMesh[Ω,"RegionMarker"->{{{innerRadius/2,-0.01},materialRegions["EggYolk"]},{{(radius+innerRadius)/2,0},materialRegions["EggWhite"]}},MeshRefinementFunction->refinementFunction]

Visualize the multi-material mesh:

mesh["Wireframe"["MeshElementStyle"->FaceForm/@{Yellow,White}]]

Solve the PDE, monitor the progress and measure the computational time and memory it takes:

measure=
AbsoluteTiming[
MaxMemoryUsed[
Monitor[solution=NDSolveValue[{op==sym,bc,ics},T,{t,0,tEnd},{r,z}∈mesh,EvaluationMonitor:>(monitor=Row[{"t = ",CForm[t]}])
],monitor]
]/(1024.^2)
];
Print["Time [s] -> ",measure[[1]],"\nMemory [MB] -> ", measure[[2]]]

Inspect the temperature range of the newly found solution:

MinMax[solution["ValuesOnGrid"]]-degreeCelsiusOffset

Create the frames for the solution:

nFrames=10;
frames=TemperatureDenatureDensityPlot[solution]/@Range[0,tEnd,tEnd/nFrames];

Rasterize the frames:

frames=Rasterize[#,ImageResolution->144]&/@frames;

Animate the frames:

ListAnimate[frames,SaveDefinitions->True]

Inspect the egg’s temperature at r=0 and z=0 after 6 minutes:

solution[6 60,0,0]-degreeCelsiusOffset

Version 6 (Use the model)

Set up the initial value of the egg in kelvins:

initT=QuantityMagnitude[Quantity[8.,"DegreesCelsius"],"Kelvins"]

Set up an initial condition for the temperature inside the egg:

ics={T[0,r,z]==initT}

Specify the temporal behavior of the boundary condition:

bcFunction=initT+(bcT-initT)Piecewise[{{1/2(1-Cos[10*Pi*t]),t<=1/10}},1]

Visualize the temporal behavior of the boundary condition:

Show[
Plot[bcFunction,{t,0,0.2},PlotRange->All]
,Graphics[{Dashed,Red,Line[{{0,bcT},{0.2,bcT}}]}]
,PlotRange->All]

Regenerate the boundary condition:

bc=HeatTemperatureCondition[ElementMarker==2||ElementMarker==5,vars,pars,<|"SurfaceTemperature"->bcFunction|>]

Solve the PDE, monitor the progress and measure the computational time and memory it takes:

measure=
AbsoluteTiming[
MaxMemoryUsed[
	Monitor[solutionFridge=NDSolveValue[{op==sym,bc,ics},T,{t,0,tEnd},{r,z}∈mesh,EvaluationMonitor:>(monitor=Row[{"t = ",CForm[t]}])],monitor]
]/(1024.^2)
];
Print["Time [s] -> ",measure[[1]],"\nMemory [MB] -> ", measure[[2]]]

Inspect the temperature range of the newly found solution:

MinMax[solutionFridge["ValuesOnGrid"]]-degreeCelsiusOffset

Create the frames for the solution:

nFrames=10;
frames=TemperatureDenatureDensityPlot[solution]/@Range[0,tEnd,tEnd/nFrames];

Rasterize the frames:

frames=Rasterize[#,ImageResolution->144]&/@frames;

Animate the frames:

ListAnimate[frames,SaveDefinitions->True]

Write a helper function to create a line of a denature temperature of a solution:

denatureLine[sol_,temp_]:=Module[{val},
val=Quiet[r/.FindRoot[sol==temp,{r,innerRadius,0,rBound-0.01rBound}]];
Line[{{val,281-degreeCelsiusOffset},{val,373.15-degreeCelsiusOffset}}]
]

Create a helper function to plot the denature positions for the two solutions at various times:

denaturePlot[t_]:=Show[
Plot[{solution[t,r,0]-degreeCelsiusOffset,solutionFridge[t,r,0]-degreeCelsiusOffset},{r,0,rBound},],
Graphics[{
{Darker@Green,denatureLine[solution[t,r,0.005],denatureTempWhite],Dashed,denatureLine[solutionFridge[t,r,0.005],denatureTempWhite]},
{Red,denatureLine[solution[t,r,0.005],denatureTempYolk],Dashed,denatureLine[solutionFridge[t,r,0.005],denatureTempYolk]},
{Black,Dashed,Line[{{innerRadius,281-degreeCelsiusOffset},{innerRadius,373.15-degreeCelsiusOffset}}]}}
]]

Visualize the temperatures at z=0.005 for the eggs starting from room and refrigerated initial temperatures and the corresponding denature positions:

denaturePlot[6 60]

Visualize the temperatures at z=0.005 for the eggs at 8 minutes:

denaturePlot[8 60]

Conclusions

Creating a simple version at first is very useful.

Increasing step by step the complexity of the Model can make things easier.

Sharp transitions or discontinuities in solution variables can cause numerical instability or accuracy loss—so use local mesh refinement or smoothing in those regions.

Having a way to compare our results with real-world data is very useful.

PDEModels/tutorial/HeatTransfer/ModelCollection/EggBoiling

ricardol@wolfram.com