TFG_VICTOR_SILIANG_RUAN_ZHAO

Vista previa del material en texto

GASTRONOMÍA MOLECULAR: 
ESTUDIO DE LAS CONDICIONES DE 
CONTORNO DE UN MODELO 
TERMOMECÁNICO BIDIMENSIONAL 
DEL SUFLÉ 
TRABAJO FIN DE GRADO PARA 
LA OBTENCIÓN DEL TÍTULO DE 
GRADUADO EN INGENIERÍA EN 
TECNOLOGÍAS INDUSTRIALES 
Febrero 2018 
Victor Siliang Ruan Zhao 
DIRECTORES DEL TRABAJO FIN DE GRADO: 
Tutor: Katerina Foteinopoulou 
Cotutor: Manuel Laso 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Daría todo lo que sé, por la mitad de lo que ignoro. 
(René Descartes) 
 
 
Agradecimientos
La realización de este trabajo no hubiera sido posible sin la siempre rápida y clara
ayuda de mi tutora Katerina Foteinopoulou a quien quiero dar las gracias profundamente
por esta oportunidad y por su paciencia con todas las explicaciones que me ha repetido
muchas veces. También a mi cotutor Manuel Laso, que ha sabido poner paz y orden a
mi mal estructurada planificación. Esta oportunidad tampoco hubiera sido posible sin
la colaboración de mi amigo Álvaro, quien en un momento de oscuridad dio luz a este
Trabajo de Fin de Grado.
Quiero agradecer a mis padres por todo lo que han hecho y sacrificado por mí. Sin sus
esfuerzos no sería la persona que soy actualmente.
A mi hermana Elisa, con quien comparto confidencias desde que tengo uso de razón,
gracias por todo.
A mis compañeros y amigos de la escuela, sin vuestra compañía y ayuda estos años
de duro esfuerzo no hubieran sido posibles. A Juan y a Nacho, por ser los amigos que
todos querrían tener en sus vidas. A Fer, del que todos los días aprendo algo. A Belén y
a Dani, por haberme enseñado tanto aunque ellos no lo sepan. A Andrés y Alberto, que
han sufrido más de una vez mi mal humor. A Hector, Iván, Alex, Patricia, Paula y Ana.
Y a Naomi con quien he compartido un largo camino en estos cortos meses, sin tu ayuda
y compañía no hubiera sido posible.
Se me olvidarán personas importantes seguro, pero a todos vosotros gracias.
1
2 Escuela Técnica Superior de Ingenieros Industriales (UPM)
Resumen
El presente Trabajo de Fin de Grado (TFG) tiene por objeto el estudio de condiciones
de contorno realistas en un modelo bidimensional termomecánico del suflé. Continua
el trabajo de dos trabajos anteriores [1] y [2] que realizaron un modelo unidimensional
y bidimensional, respectivamente. Su realización se ha llevado a cabo en el laboratorio
de Simulación de Materiales no Metálicos de la Escuela Técnica Superior de Ingenieros
Industriales, centro de la Universidad Politécnica de Madrid (ETSII-UPM).
El suflé es una preparado alimenticio que se puede asimilar a una mezcla de varias
fases en coexistencia, líquido (principalmente agua), gas (vapor de agua y aire) y sólido
(elementos como proteínas, polisacáridos, etc.). La evolución de las fases en el cocinado
real provoca la expansión del suflé a causa de la expansión del aire inicial, introducido al
batir las claras de huevo, y la generación de vapor de agua. Los experimentos llevados a
cabo en [1] con distintos grados de firmeza de la clara montada (mayor firmeza, mayor
aire introducido) confirman que el gas (aire y vapor) es el responsable del crecimiento,
siendo mayor la contribución del vapor generado.
Del trabajo de [2], se consiguió realizar un modelo de simulación del suflé bidimensio-
nal con muy buenos resultados en la predicción del cambio de fases que tiene lugar en el
crecimiento del suflé. Sin embargo, sus condiciones de contorno eran ideales, no conside-
raba flujo de calor a través de sus fronteras ni deslizamiento, por lo que era susceptible
de mejora y análisis. El modelo realizado no utiliza parámetros ajustables, es decir, se
calculan todas las incógnitas del problema en cada punto y paso temporal.
Figura 1: Varios suflés reales en un experimento de cocinado
3
El modelo está constituido por varias ecuaciones que se pueden clasificar en dos gru-
pos distintos: ecuaciones del problema del flujo, con las incógnitas ur, uz (velocidades),
T (temperatura), ρagua, ρvapor, ρaire, ρresto (densidades), Pg (presión del gas), Cp (calor
específico), k (conductividad térmica), η (viscosidad), fuente (termino de generación de
vapor) y xVg (fracción volumétrica del gas en la mezcla), y ecuaciones del problema del
mallado, con las incógnitas R y Z, figura 3.8. Las primeras se han constituido mediante
las ecuaciones de conservación de masa energía y momento (3.19) y sus correspondientes
condiciones de contorno e iniciales. Su resolución se ha realizado mediante el método de
Elementos Finitos (Finite Element Method, FEM), concretamente con la aproximación
de Galerkin de los residuos ponderados, en un código propio escrito en lenguaje Fortran
90.
Para tratar el dominio variable en el tiempo, se adoptó un esquema de generación de
malla, que se basa en métodos de mapeo de correspondencia de límites, y el mapeo se
realiza resolviendo un conjunto de ecuaciones diferenciales. Como primer paso, el dominio
físico variable en el tiempo se mapea en un dominio computacional bastante simple, fácil
de definir y constante. A continuación, los valores de las coordenadas físicas en el interior
del dominio se obtienen como las soluciones de un sistema de ecuaciones diferenciales
elípticas parciales. Debido a esta constante reconstrucción / movimiento de la malla, este
método pertenece al grupo de métodos ALE (Arbitrary-Langrangian-Eulerian), [3] [4]
[5]. La formulación de estas ecuaciones es especialmente fiable y rápida para problemas
axisimétricos con fronteras libres, grandes deformaciones y con adaptación del dominio
con el tiempo.
El dominio del flujo variable con el tiempo se mapea a un dominio constante (figura
2), el dominio de simulación, sobre el cual se calcula la evolución del mallado mediante las
ecuaciones de la sección 3.4.3.1. El dominio computacional coincide con el dominio físico
real para t = 0s.
Figura 2: Mapeo bidimensional del suflé
Se pretende realizar un estudio de las condiciones más realistas posibles por lo que
se han estudiado, el flujo de calor a través de las fronteras y deslizamiento en la pared
lateral, por efecto de la viscosidad de la mantequilla que se utiliza normalmente antes de
la introducción del suflé. La implementación numérica de estas condiciones se ha realizado
como condiciones de contorno naturales o tipo Neumann en los términos de superficie que
surgen de las ecuaciones de los residuos.
Se realizaron diferentes simulaciones con combinaciones de parámetros de operación
4 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
del horno (i.e. temperatura del parte radiativo) y de parámetros reológicos de la película
de grasa entre el soufflé y la pared. Los resultados obtenidos fueron bastante satisfac-
torios dado que se identificaron dos fenómenos competitivos en las simulaciones: (1) el
crecimiento por efecto el flujo de calor y (2) el crecimiento por efecto de la viscosidad.
Figura 3: Efectos del flujo de calor (1) y del deslizamiento (2) sobre el crecimiento del
suflé.
De forma general, el aumento de la viscosidad de la mantequilla inhibe la expansión
de los puntos de contacto de la pared siendo el valor máximo decreciente con la misma.
La temperatura de las resistencias realiza el efecto contrario, su aumento provoca el des-
plazamiento del máximo punto superior hacia el borde del molde (hacia la derecha). Su
efecto no solo desplaza de esta manera la frontera libre sino que también eleva la altura
de todo el suflé.
En conjunto, se puede afirmar que el aumento de la viscosidad reduce el efecto de la
expansión lateral favoreciendo la expansión de la parte central. La temperatura en cambio
tiene un efecto global más relevante ya que, si bien la presencia de la viscosidad es notable
en el crecimiento, la expansión es creciente con la temperatura en todos los casos, siendo
más prominente en la región central.
Los resultados de las demás incógnitas del modelo también fueron consistentes con
los fenómenos de transporte de calor y masa que ocurrenen la realidad. Se distinguieron
3 regiones diferentes dentro del suflé entre las isotermas de 70 oC y la de 100 oC, 4. La
primera marca la transición de la fase sólida, por desnaturalización de las proteínas y su
coagulación y la segunda la transición de la fase líquida de la mezcla a fase gaseosa.
La realización de este modelo de simulación en el laboratorio de Simulación de Ma-
teriales no Metálicos ha permitido el conocimiento y un análisis en profundidad de los
fenómenos físicos que tienen lugar en el cocinado real del suflé, los fenómenos que im-
pulsan su crecimiento, y su dependencia de los valores iniciales. Este modelo también es
aplicable a otros productos similares, como bizcochos o galletas, simplemente adaptando
la geometría y la composición de la receta.
Palabras clave: Material compuesto, mecánica de fluidos, métodos numéricos, trans-
ferencia de calor, ingeniería alimentaria, fenómenos de transporte, física de procesos culi-
narios, gastronomía molecular.
Víctor Siliang Ruan Zhao 5
Figura 4: Densidad del vapor para una simulación donde la temperatura radiativa era de
500oC y la viscosidad de 20 Pa s. Tiempo de simulacion t = 360s.
Códigos UNESCO:
331208 Propiedades de los Materiales
330900 Tecnología de los Alimentos
220404 Mecánica de Fluidos
220932 Termodinámica
221311 Fenómenos de Transporte
250121 Simulación numérica
120600 Análisis numérico
6 Escuela Técnica Superior de Ingenieros Industriales (UPM)
Índice general
1. Introducción 11
1.1. Motivaciones del trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2. Descripción del caso del suflé . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.3. Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
1.4. Estructura del trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2. Antecedentes 15
2.1. Antecedentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3. Modelado del suflé. Ecuaciones, condiciones iniciales y de contorno 21
3.1. Ecuaciones fundamentales del modelo . . . . . . . . . . . . . . . . . . . . . 21
3.1.1. Ecuación de conservación de la masa . . . . . . . . . . . . . . . . . 21
3.1.2. Conservación de la energía . . . . . . . . . . . . . . . . . . . . . . . 22
3.1.3. Ecuación de conservación del momento cinético . . . . . . . . . . . 23
3.2. Variación de las propiedades físicas del suflé . . . . . . . . . . . . . . . . . 25
3.2.1. Fracción volumétrica de los gases . . . . . . . . . . . . . . . . . . . 25
3.2.2. Viscosidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2.3. Calor específico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.2.4. Conductividad térmica . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.2.5. Presión de la fase gaseosa . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.6. Término fuente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3. Valores iniciales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.4. Resolución numérica del modelo . . . . . . . . . . . . . . . . . . . . . . . . 32
3.4.1. Método de los Elementos Finitos . . . . . . . . . . . . . . . . . . . 32
7
3.4.2. Aplicación del método de residuos ponderados-Galerkin . . . . . . . 33
3.4.3. Mallado. Despcripción y desarrollo . . . . . . . . . . . . . . . . . . 36
3.4.3.1. Ecuaciones del mallado . . . . . . . . . . . . . . . . . . . . 37
3.4.3.2. Frontera libre . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.5. Condiciones iniciales y de contorno . . . . . . . . . . . . . . . . . . . . . . 39
3.5.1. Condición de simetría axial. Frontera B4 . . . . . . . . . . . . . . . 40
3.5.2. Transferencia de calor. Conducción, convección y radiación . . . . . 41
3.5.3. Condición de deslizamiento . . . . . . . . . . . . . . . . . . . . . . . 44
3.5.4. Equilibrio de presiones . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.5.5. Ecuación cinemática . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.5.6. Resumen de las condiciones de contorno . . . . . . . . . . . . . . . 47
3.5.7. Implementación de las condiciones de contorno . . . . . . . . . . . . 48
3.6. Diagrama de la resolución matemática . . . . . . . . . . . . . . . . . . . . 52
4. Resultados 53
4.1. Validación del modelo con las condiciones de contorno del flujo de calor . . 53
4.2. Comparación de los efectos de la condición de deslizamiento . . . . . . . . 56
4.3. Análisis del modelo validado con condiciones de flujo de calor y deslizamiento 58
4.4. Modelo final. Empleo de flujo radiativo con factores de forma . . . . . . . . 69
4.4.1. Evolución de la velocidad . . . . . . . . . . . . . . . . . . . . . . . 70
4.4.2. Evolución de la temperatura . . . . . . . . . . . . . . . . . . . . . . 71
4.4.3. Evolución de las densidades . . . . . . . . . . . . . . . . . . . . . . 72
4.4.4. Evolución del término fuente . . . . . . . . . . . . . . . . . . . . . . 75
4.4.5. Evolución de la fracción volumétrica de gás . . . . . . . . . . . . . . 76
4.4.6. Evolución de la viscosidad . . . . . . . . . . . . . . . . . . . . . . . 77
4.4.7. Evolución de la presión . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.4.8. Evolución del calor específico . . . . . . . . . . . . . . . . . . . . . 79
4.4.9. Evolución de la conductividad térmica . . . . . . . . . . . . . . . . 80
5. Conclusiones 81
5.1. Líneas futuras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
8 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
6. Planificación temporal y espacial. Valoración de impactos 85
6.1. Estructura de Descomposición del Proyecto (EDP) . . . . . . . . . . . . . 85
6.2. Diagrama de Gantt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.3. Presupuesto económico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7. Valoración de los impactos. Aspectos relacionados con la responsabilidad
social, ética y profesional 89
Apéndices 90
Índice de Figuras 95
Índice de Tablas 97
Bibliografía 102
Víctor Siliang Ruan Zhao 9
10 Escuela Técnica Superior de Ingenieros Industriales (UPM)
Capítulo 1
Introducción
Suflé:
del Fr. soufflé, Alimento preparado con claras de huevo a punto de nieve y cocido en
el horno para que adquiera una consistencia esponjosa. [6].
1.1. Motivaciones del trabajo
Hasta hace un par de décadas los mecanismos que responden a preguntas como ¿por-
qué se vuelve blanca la clara del huevo? o ¿porqué se endurece la corteza del pan?, estaban
sin resolver o incluso no se habían llegado a formular. Surge de esta manera el término
de Gastronomía Molecular acuñado por el físico Hervé This, en colaboración con Ni-
cholas Kurti, en 1988. Esta nueva disciplina científica se definió como la "búsqueda de
los mecanismos de los fenómenos de preparación y consumo", [7], [8]. La publicación de
varios textos científicos desde entonces ha supuesto un gran avance en la ciencia de los
alimentos, una mejor compresión de los fenómenos del transporte y el desarrollo de nuevas
técnicas como el uso de aginatos, sifones o la cocina a baja temperatura.
This propuso una serie de retos para poder dar respuesta a parte de esas preguntas
[9]. Una de ellas fue la del suflé. Algunas de las incógnitas que se plantearon fueron, si
se podía hacer un suflé siempre perfecto dadas unas condiciones, si existe un límite a la
expansión o si es existe un valor de aguas mínimo para hacer el volumen máximo. Se
trata por tanto de cuestiones de gran interés académico y científico que darán explicación
a interesantes fenómenos físicos y que el presente trabajo intentará resolver.
Hasta la fecha, no existen trabajos teóricos que hayan tratado el caso del suflé que
permitan responder a estas preguntas. Sí se ha avanzado en la creación de modelos para
el pan y otros productos en donde se han identificado fenómenos físicosacoplados, es-
pecialmente fenómenos del transporte (calor y masa) cuya descripción no es sencilla y
requiere de simulaciones mediante modelos matemáticos que permitirán resolver los retos
crecientes de la ingeniería de los alimentos y que podrán aportar soluciones en la industria
alimenticia reduciendo el coste y la pérdida de energía.
El desarrollo de este trabajo también persigue una un objetivo de mayor alcance que
son la creación de herramientas que permitan realizar simulaciones futuras en el campo de
11
la ingeniería de los alimentos, desarrollando un código de elementos finitos y llegar a un
modelo validado que de respuesta a fenómenos físicos que se suceden durante los procesos
de horneado.
1.2. Descripción del caso del suflé
Un suflé está compuesto principalmente por claras de huevo (batidas a punto de nieve),
leche, harina y mantequilla, aunque otras recetas puedan incluir otros ingredientes, como
queso o chocolate. Existen, por tanto, 3 fases, líquida, sólida y gases, diferenciadas a
la hora del cocinado. Estas, tienen una fracción volumétrica inicial que variará con el
tiempo debido a los flujos de calor en el calentamiento que sucede dentro del horno. Son
de especial atención del isoterma de 100oC, por el cambio de fase del agua, contenida en
la leche y en los huevos, a vapor, y la de 70oC por el cambio de las claras de huevo a
estado sólido. La aparición de nuevas fases resulta en una clara distinción de zonas dentro
del suflé, especialmente las burbujas de vapor de agua que son responsables parciales del
crecimiento del suflé.
No es el único motivo del crecimiento del suflé (contribuye entorno al 25% [10]), la
consistencia de la claras también juega un papel muy importante a la hora del crecimiento
ya que una consistencia mayor se alcanza cuando más batidas están y mayor aire ha sido
retenido (también se introduce cierto grado de humedad, del propio aire que retrasará
la recondensación del vapor formado). Experimentalmente se demuestra que se puede
alcanzar un volumen dos o tres veces mayor que el inicial [9], [1].
Las etapas del cocinado, tras la introducción en el horno, se pueden resumir de manera
breve en los siguientes puntos:
1. Calentamiento de la paredes del molde mediante convección-conducción y de la parte
superior por radiación. Se establece un gradiente de temperatura en el interior del
suflé.
2. Distribución de la temperatura.
3. Los gradientes de presión producen un movimiento de la parte superior, frontera
libre, de acuerdo con la condición cinemática impuesta, y de las capas interiores del
suflé. Ambas ocurren por formación de burbujas de vapor.
4. La capa superior pierde todo su contenido en agua y se forma la corteza. Por debajo,
la fase vapor mantiene una presión ligeramente superior.
5. La llegada a los 70oC desde el exterior al interior del molde produce el cambio de
estado de las claras de huevo de líquido a sólido formándose la estructura interior,
terminando la cocción cuando se ha transformado por completo.
La descripción física del suflé no es trivial y su modelo de simulación, complejo. Los
distintos fenómenos físicos que se suceden en cada punto y en cada instante de tiempo
son descritos por ecuaciones diferenciales de conservación, masa, energía y momento.
Se tratan de ecuaciones diferenciales no lineales que requieren de técnicas de análisis
12 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
numérico para su resolución. Por esta complejidad se han realizado dos simplificaciones
que se han realizado en el modelo son: se considera unmodelo axisimétrico que no tiene
en cuenta la formación de nuevas burbujas y que dejerían formas toroidales alrededor del
eje de simetría, y se considera un sistema cerrado sin intercambio de masa con el
exterior y que no deja escapar los gases una vez terminada la cocción; como consecuencia
la masa, y la densidad total como suma de densidades aparentes, se conserva y se puede
considerar el suflé como un continuo. La consideración de las burbujas individuales llevaría
a la realización de un modelo tridimensional. Otra suposición, de menor impacto, es la
consideración despreciable del efecto de la gravedad sobre nuestro sistema.
En la realización de este Trabajo de Fin de Grado se ha partido de modelo previo
realizado en el departamento [2] que proponía un modelo ideal de cocción y al que se
le han incorporado nuevas condiciones de contorno (radiación, conducción, convección
y deslizamiento) más realistas con el objetivo de realizar una simulación más fiel de las
etapas de cocción del suflé. El modelo propuesto y validado en este presente trabajo no ha
sido una simple simulación de las ecuaciones sino que ha supuesto un análisis exhaustivos
de las propiedades químicas y físicas de los elementos que componen la mezcla.
Para la simulación se ha utilizado el Método de los Elementos Finitos, más concre-
tamente el método de Galerkin mixto junto con el método de Euler implícito para la
integración temporal, para la discretización espacio-temporal y un mallado estructurado,
implementado con un código propio, en lenguaje Fortran-90, escrito completamente en
el departamento de Simulación de Materiales no Metálicos y no se ha utilizado ningún
sofware comercial para el cálculo y simulación del mismo. El modelo no contiene ningún
parámetro ajustable y todas las propiedades físicas son calculadas teniendo en cuenta si
dependencia con la temperatura el tiempo y la posición.
Víctor Siliang Ruan Zhao 13
1.3. Objetivos
Se persiguen los siguientes objetivos:
Comprender los fenómenos físicos que suceden en el proceso de cocinado real del
sufle.
Desarrollar un código de implementación de las condiciones de contorno realista
mediante elementos finitos.
Validar el nuevo modelo con las nuevas condiciones de contorno incorporadas (ra-
diación , conducción, convección y deslizamiento)
Simular el modelo y realizar un análisis paramétrico para investigar el efecto de las
paredes y de las condiciones de operación.
Verificar la validez de las simulaciones.
1.4. Estructura del trabajo
En este Trabajo de Fin de Grado sigue la siguiente estructura:
En el capítulo 2, se realiza una descripción de los artículos más destacados y rele-
vantes en este campo.
En el capítulo 3, se detallan las ecuaciones (masa, energía, momento), las propieda-
des físicas y las condiciones iniciales y de contorno.
En el capítulo 4, se presentan las principales simulaciones realizadas y se analizan
los efectos de la aplicación de las nuevas condiciones de contorno introducidas.
En el capítulo 5, se discutirán las resultados obtenidos para validar el modelo, así
como las posibles líneas futuras de investigación
En el capítulo 6, se concluirá con las planificación temporal y espacial incluyendo
un breve análisis de los impactos del presente trabajo
14 Escuela Técnica Superior de Ingenieros Industriales (UPM)
Capítulo 2
Antecedentes
2.1. Antecedentes
El uso de modelos de elementos finitos en la industria alimenticia, sobretodo de mo-
delos CFD, está bastante popularizado para ciertos productos alimenticios pero no se ha
realizado, en cuanto al conocimiento del autor, ningún trabajo de igual magnitud para
el caso del suflé. Existen sobretodo artículos acerca del proceso de horneado del pan, un
producto que es similar a nuestro suflé y del que se ha realizado una extensa documen-
tación, [11], [12], [13], [14], [15], [16]. Su similitud se encuentra en la existencia de varias
fases y su evolución temporal así como el movimiento de una frontera libre (corteza del
pan).
Este trabajo es la continuación de dos anteriores TFGs realizados en el departamento
y de los cuales se parte para seguir extendiendo el modelo. El primer modelo propuesto fué
realizado en 2015 por Teresa Díez [1], que realizó dos modelos, en dimensión 0 y 1, y que
fueron validados experimentalmente. El modelo de dimensión cero fué caracterizado por el
volumeny su evolución realizando una primera aproximación, bajos fuertes suposiciones
como transferencia de calor y masa en tiempo finito, que permitieron comprender mejor el
problema aunque se presentaban ciertas limitaciones debido a los fenómenos instantáneos
o nulos que se sucedían. En el modelo unidimensional, se introdujeron el cálculo de las
velocidades, temperaturas, masas y presiones, dependientes de las posición y del tiempo,
que complicaron las ecuaciones del modelo pero que se aproximaba más a la realidad del
proceso de cocción. Se analizó la variación de la altura del suflé que mostraban, cuantitati-
vamente, resultados muy similares con los observados en los experimentos de cocción. Sin
embargo, al no contemplarse la salida de gases tras el proceso de cocción, no se obtuvieron
resultados de a esta parte. Otros parámetros como la variación de la temperatura o de
la masa mostraron resultados congruentes con la realidad, identificándose los fenómenos
de evaporación del agua y generación de vapor que contribuye a la expansión de la que
hablará en los próximos párrafos.
En 2017, se extendió el modelos unidimensional a un modelo bidimensional, en reali-
dad tridimensional con simetría axial, realizado por Sebastian Kramarz [2]. Este modelo
bidimensional añadió dificultad y complejidad a las ecuaciones así como al modelo de si-
mulación de elementos finitos; se impusieron, para un tratamiento inicial, unas condiciones
15
de contorno de tipo Dirichlet 1 en las fronteras del modelo de calentamiento instantáneo
que suponen el punto de partida del presente trabajo. Se consiguió realizar de forma sa-
tisfactoria un marco bidimensional de las ecuaciones de conservación del modelo y de las
ecuaciones del mallado. A pesar de las condiciones ideales, se logró observar que en el
"suflé computacional"se sucedían los fenómenos esperados de expansión en las zonas más
calientes por generación de vapor y zonas de altos gradientes de velocidad que ocurren en
las regiones con más temperatura y libertad de movimiento, la frontera libre. El trabajo
realizado asentó las bases para una mejor compresión de los fenómenos acoplados dentro
del suflé a falta de hacer un estudio de unas condiciones de contorno realistas y de los
efectos de las mismas.
Figura 2.1: Simulación de la velocidad de los puntos del suflé. [2]
Al igual que en los anteriores trabajos, en este no se utilizan parámetros ajustables,
sino que todas las propiedades son calculadas en cada punto del mallado y en cada paso
de tiempo, y son consecuencia de los experimentos sobre las propiedades físicas de los
elementos. Se ha partido también de la receta básica (harina, leche, huevos, sal y man-
tequilla), donde todos los elementos se han asimilado a 4 sustancias: agua, vapor, aire y
resto (conjunto de elementos distintos de los anteriores por fase u otra razón, p.ej. sal,
proteínas, polisacáridos, etc.). Parte de los huevos, las claras fundamentalmente, y la leche
se asimilaron a agua, para que el modelo fuera más sencillo. El vapor es una razón del
crecimiento del suflé y por su generación en el proceso de calentamiento era imprescindi-
ble su consideración. Otra de las razones del crecimiento es la existencia de aire húmedo
1 Nota: Una condición de contorno tipo Dirichlet para una ecuación en derivadas parciales como
∇2y + y = 0 en un dominio Ω ⊂ Rn, es de la forma:
y(x) = f(x) ∀x ∈ ∂Ω (2.1)
donde ∂Ω denota la frontera de Ω. Otro tipo de condición de contorno es el de tipo Neumann, en el que
dada la ecuación en derivadas parciales anterior se impone una condición a la frontera ∂Ω del tipo:
∂ y
∂n
(x) = f(x) ∀x ∈ ∂Ω (2.2)
donde n denota la normal exterior a la frontera ∂Ω.
16 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
dentro de las claras de huevo batidas por lo que también fue considerada. Las proteínas
de las yemas de huevo, que juegan un papel esencial ya que solidifican a los 70oC, así
como otros componentes fueron asimilados a un resto.
En relación al estado del arte de las simulaciones sobre alimentos a través de modelos
matemáticos se ha venido centrando sobre todo en el estudio del horneado del pan y
los fenómenos que lleva asociados. Mondal y Datta realizaron una recopilación del estado
de las investigaciones sobre el pan en 2007 [17] donde recogieron los principales trabajos
de las últimas dos décadas. El horneado de pan es un proceso que envuelve cambios de
temperatura, de volumen y de humedad que se encuentran fuertemente acoplados [13]. Son
de vital importancia los fenómenos del transporte de masa y calor [12], ya que durante el
horneado se suceden distintos estados, la masa húmeda se transforma en miga para después
convertirse en corteza en la parte exterior. El transporte provoca una pérdida de masa,
producción de CO2, expansión de la fase vapor e incremento del volumen, solidificación
de las proteínas, entre otros efectos que se encuentran muy correlacionados en la cocina.
Además, los sistemas físicos que representan no son sólidos continuos y no se tratan de tal
manera sino que debido a las discontinuidades deben ser estudiado como un medio poroso.
Los estudios Datta [18], [19] y de Dhall y Datta , [20], sobre el transporte en materiales
alimentarios, como el pan, propusieron modelos para el tratamiento de estos medios.
Figura 2.2: Mallado axisimétrico del pan [11]
Es muy relevante el trabajo de Ousegui et al [21] que realizó un modelo multifase en
medio poroso con una formulación que permitía observar el efecto de un frente de evapo-
ración y su influencia sobre la expansión así como los fenómenos difusivos de ambas fases.
En la misma línea, los trabajos de Purlis y Salvadori [14], [15] y [16] se habían realizado
previamente un estudio del mecanismo de condensación-evaporación que es el responsable
del rápido calentamiento de la matriz porosa. Concluyeron que su existencia aparece en es-
tructuras porosas tanto cerradas (interior del pan, masa húmeda) como abiertas (corteza)
hasta que se alcanza la temperatura de 100 oC momento en el se produce la evapora-
ción del agua contenida en el interior. El avance de este frente de evaporación es esencial
para el proceso de horneado ya que produce un cambio en la estructura tanto interna,
producción de miga, como externa formación de corteza que previene la deshidratación.
Esta isoterma formará un frontera entre la fase líquida y la fase gaseosa al igual que otras
Víctor Siliang Ruan Zhao 17
como la de 70 oC que separa en alimentos que contienen huevos las fases liquidas de las
sólidas por la desnaturalización (ruptura de la estructura proteica de las claras del huevo
y su posterior coagulación) de las proteínas.
Los modelos propuestos en los trabajos se centraban principalmente en el estudio de los
fenómenos de transporte de masa y calor en medio poroso, de evaporación-condensación
y la sucesión de la frontera de evaporación. Trabajos más recientes han estudiado la
influencia de la viscosidad inicial y la humedad en la deformación de la frontera libre
que constituye la corteza del pan [22]. También se han realizado estudios para conocer la
evolución de las propiedades físicas dentro del pan como contenido y pérdida de agua o
porosidad de la miga [23], [24].
La relación de los fenómenos del transporte con la deformación y consecuente movi-
miento de la frontera libre se ilustra en la figura 2.3.
Figura 2.3: Esquema de relación entre los fenómenos de transporte de calor y la deforma-
ción de [11]
Las ecuaciones que se presentan en los modelos de estudio son, en la mayoría de los
autores, ecuaciones de conservación, masa, energía y momento o ecuaciones como la ley de
Fick o la de Darcy; se tratan pues de ecuaciones diferenciales en derivadas parciales que
requieren de métodos numéricos para su resolución. En la literatura consultada, se emplea
con mucha extensión el método de los elementos finitos (FEM, por si siglas en inglés) y
algunos estudiosmediante el método de diferencias finitas [25]. El empleo de códigos FEM
se encuentra muy extendido en este la simulación en el campo de la ingeniería alimentaria
[26], ya sea mediante el uso de programas comerciales como COMSOL Physics o ANSYS
Fluent, como de códigos propios en el caso de este presente trabajo. Su amplio empleo
es debido a la robustez de las soluciones que se alcanzan, sobretodo en con esquemas de
Euler implícitos en el tiempo, y por la rápida convergencia de los resultados. La forma
de abordar las fronteras libres es una novedad de este trabajo, al igual que los anteriores
[2] y [1], ya que emplea un enfoque ALE (Arbitrary-Lagragian-Eulerian) para capturar el
movimiento de la misma durante el proceso de simulación. La principal ventaja en este
18 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
tipo de método, frente a otros como el Euleriano o Lagrangiano tradicionales, es que no
necesitan generar un nuevo mallado en cada paso de tiempo porque su evolución está
dada por una ecuación diferencial propia.
El uso de un código propio en lugar de uno comercial, de pago de libre uso, radica
en las condiciones de nuestro modelo que necesitaba de condiciones específicas de trabajo
para problemas axisimétricos y con gran deformación de los cuales se disponía de una
literatura amplia como [27], [28], [29] y [5]. Además se trata de códigos adaptados a este
tipo de problema, que no ad hoc para este modelo, y que por tanto son más rápidos y
eficientes en su resolución.
Víctor Siliang Ruan Zhao 19
20 Escuela Técnica Superior de Ingenieros Industriales (UPM)
Capítulo 3
Modelado del suflé. Ecuaciones,
condiciones iniciales y de contorno
Este capítulo está dedicado al desarrollo teórico del modelo. Se empezará por desarro-
llar las ideas del anterior capítulo presentando las ecuaciones fundamentales, conservación
de masa, energía y momento cinético. Las propiedades físicas del modelo definidas se en-
contrarán en el segundo apartado. Se presentarán los valores iniciales recogidos por [1]
en el tercer apartado. El desarrollo matemático, según [2] constituye el cuarto apartado
desarrollando y ampliando las condiciones de contorno realista (transferencia de calor,
deslizamiento, equilibrio de presiones y la ecuación cinemática de la frontera libre). Por
último, se recogen algunos diagramas de flujo para dar explicación a como es la imple-
mentación del código propio realizado.
3.1. Ecuaciones fundamentales del modelo
3.1.1. Ecuación de conservación de la masa
La primera de las ecuaciones del modelo es la conservación de masa [30].
∂ρ
∂t
+∇ · (ρv) = ±fuente (3.1)
El término fuente de esta ecuaciones representa la transformación de la fase líquida
(agua) a la fase vapor que será positivo para la generación de vapor o negativo en el caso
de que el líquido se reduzca.
Considerando las 4 fases de las hipótesis iniciales, agua, vapor, aire y resto la
ecuación para cada una de ellas es:
∂ρagua
∂t
+∇ · (ρaguav) = −fuente (3.2)
∂ρvapor
∂t
+∇ · (ρvaporv) = +fuente (3.3)
21
∂ρaire
∂t
+∇ · (ρairev) = 0 (3.4)
∂ρresto
∂t
+∇ · (ρrestov) = 0 (3.5)
De la consideración de sistema cerrado tendremos que las ρtotal será la suma de las 4
distintas densidades aparentes de las 4 fases, considerando así el suflé como un continuo:
ρtotal = ρagua + ρvapor + ρaire + ρresto (3.6)
El cálculo de estas ecuaciones nos proporcionará información acerca de la como es la
formación del vapor, causa de la elevación del suflé, y también de la evolución de la fase
sólida, por transformación de las proteinas de los huevos, que formará la estructura al
final del cocinado.
El cálculo y la forma del término fuente se recogen en el apartado 3.2.6.
3.1.2. Conservación de la energía
La temperatura superior del interior del horno provoca una transferencia de calor
hacia el suflé y la energía intercambiada en el proceso aumentará la temperatura de los
todos los puntos. El intercambio y aumento de la temperatura obedecerá la ecuación de
conservación:
∂
∂t
(ρtotal Cp T ) +∇ · (ρtotalCp Tv) = −∇ · q +Hv (3.7)
El flujo de calor, q, viene dado por la ley de Fourier:
q = −k∇T (3.8)
donde k es la conductividad térmica y ∇T el gradiente térmico dentro de los puntos del
suflé. Cp es el calor específico del suflé.
El término Hv, es de la forma:
Hv = −fuente Lv (3.9)
donde Lv es el calor latente de vaporización del agua tomado como 561,05 cal/kg y por
tanto la expresión anterior representa la pérdidas de calor por cambio de fase del agua.
Estos parámetros dependerán de la composición y de la temperatura del punto del
suflé (dependencia espacial y temporal) por lo que serán calculados en cada instante de
tiempo y en cada punto aumentando la dificultad del problema.
22 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
3.1.3. Ecuación de conservación del momento cinético
La primera de las simplificaciones fue la consideración de un problema bidimensional de
simetría axial en un sistema de coordenadas cilíndricas, como se puede ver en la siguiente
figura:
Figura 3.1: Varios suflés en un experimento de cocción
Nuestras variables espaciales son por tanto r y z. Suponiendo las fuerzas gravitatorias
despreciables la ecuación de conservación del momento cinético es de la forma:
∂(ρtotal v)
∂t
+∇ · (ρtotal v v) = −∇P +∇ · τ(t, r) (3.10)
Debido a la consideración del suflé como medio continuo, podemos aplicar esta ecuación
al total de la masa e introducir ρtotal de la ecuación (3.6). La velocidad v será una función
de la forma v = v(r, z).
El fluido considerado es un fluido newtoniano y el tensor de tensiones en esta ecuación
estará dado por su ecuación constitutiva [31]:
τ = −η(t, r)[∇ v + (∇ v)t] (3.11)
La viscosidad, η, cambiará como se verá en el apartado de propiedades físicas 3.2
de la concentración volumétrica de gas en la mezcla que cambiará en cada punto por el
calentamiento dentro del horno(dependencia espacial), también dependerá del tiempo por
el aumento de temperatura. La parte isotrópica del tensor de tensiones, proporcionará la
presión interior debida a la fase gaseosa:
σ = −PI + τ (3.12)
Considerando la fase gaseosa como una mezcla de gases ideales, la presión se reduce a
la expresión siguiente:
P = Pgas(t, r) = Paire + Pvapor = RT
(
ρiaire(t, r)
Maire
+
ρivapor(t, r)
Mvapor
)
(3.13)
Víctor Siliang Ruan Zhao 23
Hay que señalar que la densidad que aparece en esta expresión es la densidad intrínseca,
ρi, la densidad de cada uno de los elementos como sustancia que está relacionada con la
densidad aparente por la siguiente expresión:

ρik =
ρk
xVk
k = agua y resto
ρik =
ρk
xVg
k = aire y vapor
(3.14)
Las densidades intrísecas del agua y del resto se tomaron como constantes durante la
simulación pero las del aire y el vapor se calcularon en cada punto y en cada instante de
tiempo como se verá en el apartado de propiedades físicas 3.2.Maire yMvapor, representan
el peso molecular del aire y el vapor.
La ecuación (3.10) se desarrolla según cada componente espacial r , z de forma que:
Componente r:
er ·
∂(ρtotal v)
∂t
+ er · ∇ · (ρtotal v v) = − er · ∇P + er · ∇ · τ (3.15)
Componente z:
ez ·
∂(ρtotal v)
∂t
+ ez · ∇ · (ρtotal v v) = − ez · ∇P + ez · ∇ · τ (3.16)
Finalmente se convierten en:
Componente r:
vr
∂ρtotal
∂t
+ ρtotal
∂vr
∂t
+ 2ρtotal vr
∂vr
∂r
+ ρtotal uz
∂vr
∂z
+ v2
r ρtotal
∂ρtotal
∂r
+
ρtotal v
2
r
r
+ vr uz
∂ρtotal
∂z
+ vr ρtotal
∂ρtotal
∂z
= −∇ ·
(
Per
)
+
P
r
+∇ ·
(
τ · er
)
− τθθ
r
(3.17)
Componente z:
uz
∂ρtotal
∂t
+ ρtotal
∂vz
∂t
+ 2ρtotal uz
∂vz
∂z
+ ρtotal vr
∂vr
∂r
+ v2
z
∂ρtotal
∂z
+
ρtotal vr vz
r
+ vr uz
∂ρtotal
∂r
+ uz ρtotal
∂vz
∂r
= −∇ ·
(
Pez
)
+∇ ·
(
τ · ez
)
(3.18)
Las ecuaciones de conservación, por ser ecuaciones diferenciales, tendrán que estar
definidas en un dominio Ω confrontera Γ en que se definirán las condiciones de contorno
al igual que unas condiciones iniciales. Las condiciones iniciales se resumen en las tablas
3.6 y 6.1 y las condiciones de contorno en tabla 3.12.
24 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
Las siguientes ecuaciones formarán el problema del flujo, representado en la figura 3.8:
∂ρagua
∂t
+∇ · (ρaguav) = −fuente
∂ρvapor
∂t
+∇ · (ρvaporv) = +fuente,
∂ρaire
∂t
+∇ · (ρairev) = 0
∂ρresto
∂t
+∇ · (ρrestov) = 0
∂
∂t
(ρtotal Cp T ) +∇ · (ρtotalCp Tv) = −∇ · q +Hv
vr
∂ρtotal
∂t
+ ρtotal
∂vr
∂t
+ 2ρtotal vr
∂vr
∂r
+ ρtotal uz
∂vr
∂z
+ v2
r ρtotal
∂ρtotal
∂r
+ ρtotal v
2
r
r
+
+vr uz
∂ρtotal
∂z
+ vr ρtotal
∂ρtotal
∂z
= −∇ ·
(
Per
)
+ P
r
+∇ ·
(
τ · er
)
− τθθ
r
uz
∂ρtotal
∂t
+ ρtotal
∂vz
∂t
+ 2ρtotal uz
∂vz
∂z
+ ρtotal vr
∂vr
∂r
+ v2
z
∂ρtotal
∂z
+ ρtotal vr vz
r
+
+vr uz
∂ρtotal
∂r
+ uz ρtotal
∂vz
∂r
= −∇ ·
(
Pez
)
+∇ ·
(
τ · er
)
(3.19)
3.2. Variación de las propiedades físicas del suflé
En el modelado del suflé, ecuaciones (3.19), se incluyen adicionalmente más paráme-
tros que deben ser cálculados en cada punto e instante temporal. Constituyen una parte
esencial que da realismo a las simulaciones y que aporta información fidedigna de como es
la evolución de interior de las misma. Gran parte de estas propiedades han sido calcula-
das mediante correlaciones o fórmulas empíricas que son válidas para el rango de trabajo
de la simulación (sobre todo para la temperatura). Los parámetros inciales utlizados en
este trabajo han sido tomados de [1] y [2] y se amplia con nuevas propiedades y mejores
condiciones de trabajo.
3.2.1. Fracción volumétrica de los gases
La densidad intrínseca definida en la ecuación (3.14) para el aire y el vapor, también se
puede utilizar para calcular las fracciones de vapor del agua y el resto, teniendo en cuenta
que las densidades intrínsecas de estos componentes serán constante en la simulación
aunque no sus densidades aparentes. Mediante estas fracciones volumétricas se puede
calcular la fracción del gas, mezcla de aire y vapor:
xVg = 1− xVagua − xVresto (3.20)
Víctor Siliang Ruan Zhao 25
3.2.2. Viscosidad
La expresión utilizada para el cálculo de la viscosidad se ha obtenido de [32] y [33].
Se trata de una expresión que considera la viscosidad del suflé como dependiente de la
viscosidad del líquido, ηl y como de la fracción volumétrica del gas, xVg , que tendrá un
comportamiento creciente con la misma.
η = ηl
3π
16
[
1−
(
xVg
0,74
)1/3
] (3.21)
Esta expresión depende principalmente de xVg ya que la viscosidad del líquido toma
valores entre [10−104]Pa s. La representación de la expresión se encuentra en la figura 3.2.
Se puede apreciar que para valores del entorno de xVg = 0,74 la expresión diverge; para
estos casos los valores han sido truncados al valor máximo 13017,8Pa s, que es dependerá
del valor escogido de ηl, en el modelo de valor 1000Pa s.
0 0,1 0,2 0,3 0,4 0,5 0,6 0,74 1
0
0,5
1
1,5
2
2,5
3
3,5
4
4,5
5
·104
xVg
η
(P
a
s)
Viscosidad del suflé para ηl = 1000Pa s
ηreal
ηtruncada
Figura 3.2: Viscosidad del suflé en función de la fracción volumétrica de gas
26 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
3.2.3. Calor específico
La correlación para el cálculo del calor específico se tomó de [34] y es de la forma:
Cp,k = Ak +Bk T + Ck T
2 +Dk T
3 k = agua, vapor, aire y resto (3.22)
Ak, Bk, Ck, Dk representa constantes según la tabla 3.1. El rango de validez es amplio
[273,15− 1500]K.
Sustancia A B C D
Agua 18.144 - - -
Vapor de agua 6.97 0.3464E-2 0.04833E-5 -
Aire 6.557 0.1477E-2 0.02148E-5 -
Resto 10.8863 - - -
Tabla 3.1: Constantes de la correlación del calor específico en
kcal
kmolK
, para temperaturas
dadas en K
El calor específico del resto se aproximó al 60% del calor específico del agua: Cp,resto =
0,6Cp,agua. El calor específico del suflé se calculó como la suma ponderada a las fracciones
másicas de las sustancias consideradas (regla de mezcla) [31]
Cp = Xm,aguaCp,agua +Xm,vapor Cp,vapor +Xm,aireCp,aire +Xm,restoCp,resto (3.23)
La fracción másica de cada sustancia se calcula con la siguiente expresión:
Xm,k =
ρk
ρtotal
k = agua, vapor, aire y resto (3.24)
3.2.4. Conductividad térmica
Al igual que con el calor específico, la conductividad térmica se calculará con una regla
de mezla (Reuss), [31]):
ks =
k∑
i=agua
xViki i = agua, vapor, aire y resto (3.25)
Las fracciones volumétricas xVi se calculan con la ecuación 3.14. La conductividad de
cada sustancia se obtuvo de [35] y se resumen en la tabla 3.2.
Las ecuaciones (3.21), (3.23) y (3.25) dependen de las densidades y fracciones volumé-
tricas de los componentes y son parámetros que entran en las ecuaciones de conservación.
Víctor Siliang Ruan Zhao 27
Por tanto, tienen que ser calculados en cada punto y en cada paso de tiempo como se
indica en la figura 3.8.
Sustancia Conductividad térmica
W
mK
Agua 0.6
Vapor 0.0179
Aire 0.023
Resto 0.6
Tabla 3.2: Conductividad térmica de las sustancias según [35]
El cálculo de los términos anteriores permite calcular la presión de la fase gaseosa y el
término fuente que aunque no se traten del propiedades de físicas se han incluido en esta
sección por tener una gran dependencia de las mismas y ser variables de la ecuación de
conservación de la energía.
3.2.5. Presión de la fase gaseosa
Se calcula mediante la ecuación (3.13) y la ecuación (3.14), para el aire y para el vapor.
3.2.6. Término fuente
El término fuente introducido en la ecuación (3.1) representa el cambio de fase del
agua de la mezcla en vapor. Depende fundamentalmente de la temperatura pero también
del gradiente de concentración. Está dado por la siguiente expresión:
fuente = a Jg =
6xVg
dp
kg(csat − c) (3.26)
La tabla 3.4 resume los principales parámetros que interviene en el cálculo del término
fuente. El coeficiente de transferencia de materia, kg, se calcula a través del número de
Sherwood; en una burbuja de dimensión característica dp, se tiene una aproximación de
∼= 10 según [36]. Para la difusividad del vapor en el aire Dvapor,aire encontramos en [35]
un valor aproximado de ∼= 2,92 ∗ 10−5m2/s. Si se toma como valor del diámetro de las
burbujas dp = 10−3m se tiene por tanto que el valor de kg es 0,292m/s.
La concentración de saturación del vapor de agua, csat, se ha calculado mediante la
Ley de Raoult que establece la relación con la presión de vapor de la sustancia en
una mezcla ideal. Por este motivo tampoco se ha considerado el uso de un modelo de
coeficiente de actividad que quedará para futuras línea de investigación y que por lo tanto
en este modelo es de valor igual a la unidad. Se considera que la fracción molar de agua
en líquido, xmolagua , es un valor constante, ∼= 0,85, y que no varía con el tiempo dado que
el resto de componentes tienen un peso molecular muy inferior al del agua.
28 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
Jg Flujo de vapor de agua kg(csat − c)
kg Coeficiente de transferencia de materia
ShDvapor,aire
dp
a Área interfacial específica total entre gas-líquido π d2
p burbujas
burbujas Número de burbujas por m^3
xVg
(π d3
p /6)
csat Concentración de saturación
Pv xmolaguaMv
RT
dp Diámetro media de burbujas ∼= 10−3m
Sh Número adimensional de Sherwood ∼= 10
Dvapor,aire Difusividad del vapor en el aire ∼= 2,92 ∗ 10−5m2/s
xmolagua Fracción molar del agua en el líquido ∼= 0,85
Tabla 3.3: Términos de la ecuación (4.30)
La presión de vapor del agua estará dada por la ecuación de Antoine:
log10Pv(mmHg) = A− B
T + C
(3.27)
Los parámetros utilizados en la ecuación se recogenen la tabla (3.4) de [37] :
Sustancia A B C Tmin (oC) Tmax (oC)
Agua 8.07131 1730.63 233.426 1 100
Agua 8.14019 1810.94 244.485 99 374
Tabla 3.4: Valores de la ecuación de Antoine para el agua
Por último la concentración del vapor de agua es la densidad intríseca que se calcula
mediante la ecuación (3.14).
Víctor Siliang Ruan Zhao 29
Se recogen en la tabla 3.5 las expresiones que constituyen todas las propiedades físicas
del modelo del suflé.
Símbolo Propiedad Expresión
η Viscosidad ηl
3π
16
[
1−
(
xVg/0, 74
)1/3
]
Cp,k Calor específico Ak +Bk T + Ck T
2 +Dk T
3
ks Conductividad térmica
resto∑
i=agua
xVi ki
xVg Fracción volumétrica del gas 1− xVagua − xVresto
Pgas Presión de la fase gaseosa RT
(
ρiaire
Maire
+
ρivapor
Mvapor
)
fuente Término fuente de la ecuación de masa, tabla 3.3 a Jg =
6xVg
dp
kg(csat − c)
Tabla 3.5: Propiedades físicas del modelo y expresiones
30 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
3.3. Valores iniciales
Se recogen los valores iniciales de los parámetros de los trabajos llevados a cabo por
[1] y [2] y que son el resultados del estudio de los elementos básicos de los ingredien-
tes. En las tablas 3.6 6.1 se detallan los valores iniciales a los que se llegaron mediante
experimentación y cálculo de propiedades físicas y que no son parámetros ajustables.
Ingrediente Cantidad(kg) Densidad (kg/m3) Agua (%)
Mantequilla 0.225 911 0.25
Harina 0.02 550 14.1
Leche 0.26 1036 91
Sal 0.045 2200 -
Claras 0.035 1070 88
Claras al punto de nieve 0.16 207.5 -
Total 0.7045 - -
Tabla 3.6: Masas, densidades y% de agua iniciales de los ingredientes de [1]
Las densidades aparentes se calcularon mediante la ecuación 3.14 al igual que se men-
cionó en el apartado 3.1.3. La suma de las densidades aparentes según 3.6 da la ρtotal.
Conocidas las fracciones másicas iniciales, el contenido en agua y que las claras a punto
de nieve el mismo contenido que las claras sin montar, podemos calcular la masa de agua,
magua:
magua =
clara∑
i=mantequilla
miwi (3.28)
Sustancia densidad(kg/m3) m(kg) Xm V (m3) xV
Agua 247.2124 0.27584 0.55446 0.000281966 0.25269
Vapor 0.071454 0.000079731 1.60261 0.000621553 0.55703
Aire 0.5346 0.0005966 0.001991 0.000621553 0.55703
Resto 198.5758 0.22157 0.44537 0.000212303 0.19026
Total 445.85978 - - 0.0011 -
Tabla 3.7: Masas, densidades aparentes, fracciones molares y volumétricas y volumen de
las sustancias del modelo
El suflé inicialmente tiene la forma cilíndrica del molde con una altura de 3 cm y
un radio de 5 cm. Las condiciones iniciales del ambiente son de 25 oC y 1 atm.
Víctor Siliang Ruan Zhao 31
3.4. Resolución numérica del modelo
Se presenta un breve desarrollo teórico de los métodos utilizados para la resolución
numérica de las ecuaciones del modelo así como la forma final a resolver tras haber
implementado el método de elementos finitos.
3.4.1. Método de los Elementos Finitos
La idea fundamental del Método de Elementos Finitos es la división de un dominio
continuo en partes, aproximando las funciones a resolver en cada una por polinomios.
Se trata de un método muy extendido que empezó en el campo del análisis estructural
que transformó el sólido continuo en elementos discretos más fáciles de manejar. En la
actualidad, se utilizan en mucho campos de la ciencia como la ingeniería de materiales o
mecánica de fluidos. Es especialmente adecuado para simulación de procesos en los que
son esperadas grandes deformaciones convirtiéndolo en un método muy apropiado para
problemas del campo de la ingeniería de alimentación.
Generalizando el problema [38], la formulación de elementos finitos busca resolver el
problema:
Sea un conjunto de ecuaciones diferenciales A(u) definido en un dominio Ω tal que
A(u) =

A1(u)
A2(u)
.
.
 = 0 (3.29)
con unas condiciones de contorno B(u), definidas en la frontera de Ω, Γ tal que
B(u) =

B1(u)
B2(u)
.
.
 = 0 (3.30)
se busca encontrar la solución u que satisfaga ambas ecuaciones y las condiciones de con-
torno ∀x ∈ Ω. Como esta formulación es bastante fuerte buscamos formular el problema
de una manera menos restrictiva, la formulación débil.
La formulación débil del problema expresa las ecuaciones diferenciales como integrales
que podrán ser tratadas en un espacio vectorial y en donde la solución u es de la forma:
u ' û =
n∑
i=1
Ni ui = Nû
donde u y û son la solución exacta y aproximada, N las funciones bases (o de forma)
y ui las incógnitas de nuestro problema.
Sean v y v̄ dos conjuntos de funciones arbitrarias de dimensión igual al número de
32 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
ecuaciones diferenciales del problema y las condiciones de contorno, multiplicando e inte-
grando tendremos que:
∫
Ω
vTA(u)dΩ ≡
∫
Ω
[v1A1(u) + v2A2(u) + v3A3(u) + ...]dΩ ≡ 0 (3.31)
∫
Γ
v̄TB(u)dΓ ≡
∫
Γ
[v̄1B1(u) + v̄2B2(u) + v̄BA3(u) + ...]dΓ ≡ 0 (3.32)
De manera que la formulación de elementos finitos queda definida como:∫
Ω
vTA(u)dΩ +
∫
Γ
v̄TB(u)dΓ = 0 (3.33)
El dominio se subdividirá en elementos, Ωe, conectados mediante nodos en los que se
aproximará la solución de nuestras ecuaciones, u.
3.4.2. Aplicación del método de residuos ponderados-Galerkin
Las funciones v y v̄ son arbitrarias por lo que podemos elegir funciones de la forma:
v ' v̂ =
n∑
i=1
wi vi
La ecuación (3.33) quedaría de la forma:
∫
Ω
wT
i A(Nû)dΩ +
∫
Γ
w̄T
i B(Nû)dΓ = 0 (3.34)
donde A(Nû) representa el residuo o error obtenido por sustitución de la aproximación
en la ecuación diferencial y B(Nû) el residuo de las condiciones de contorno. La ecuación
(3.34) es la integral ponderada de los residuos y los términos wi los pesos de integración.
El método de Galerkin es un método de residuos ponderados que utiliza como pesos de
integración las mismas funciones base de la aproximación, Ni. La aproximación realizada
por este método es la más intuitiva y fácil de aplicar siendo la más popularizada tanto en
los códigos propios como en programas comerciales. El mallado utilizado es estructurado
con elementos triangulares cuadráticos, sobre los cuales se definen las funciones base.
Serán cuadráticas para la velocidades y la temperatura y lineales para las densidades.
Las funciones base lineales, dentro del elemento triangular lineal estarán definidas por
3 nodos; las funciones cuadráticas por 6. Toman valor unidad cuando la función está
definida en su nodo y nulo en caso contrario, es decir, φi =
{
1 en el nodo i
0 en el resto
. Las
Víctor Siliang Ruan Zhao 33
derivadas de estas funciones son indeseables ya que, las primeras son continuas a trozos
(discontinuas), y las segundas similares a deltas de Dirac (nulas o infinitas) e introducirían
discontinuidades en el dominio del problema.
Definimos el elemento de referencia T̂ con coordenadas locales como T̂l : [ξ, η] :=
[(0, 0), (0, 1), (1, 0)] en caso de elementos lineales. Para elementos cuadráticos será el T̂q :
[ξ, η] := [(0, 0), (0, 1/2), (0, 1), (1/2, 0), (1/2, 1/2), (1, 0)], como puede verse en la figura 3.3.
Figura 3.3: Elemento de referencia T̂q
Una función base definida en este dominio tendrá la forma:
Funciones lineales
ψi =

1− (x+ y) i = 1
x i = 2
y i = 3
(3.35)
Funciones cuadráticas
φi =

1− (x+ y))(1− 1(x+ y)) i = 1
4x(1− (x+ y)) i = 2
x(2x− 1) i = 3
4y(1− (x+ y)) i = 4
4xy i = 5
y((2x− 1) i = 6
(3.36)
Aplicando la formulación débil de la ecuación (3.33) a las ecuación del modelo, (3.19),
con las funciones base (3.35) y (3.36) obtemos los siguiente residuos cuyo desarrollo:
Residuos de las ecuaciones de conservación de masa:
(3.37)
Ri
ρagua =
∫
Ω
ψi
[
∂ρagua
∂t
+ ur
∂ρagua
∂r
+ uz
∂ρagua
∂z
+ ρagua
ur
r
+ ρagua
∂ur
∂r
+ ρagua
∂uz
∂z
+ fuente
]
dV
34 Escuela Técnica Superior de Ingenieros Industriales (UPM)GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
(3.38)
Ri
ρvapor =
∫
Ω
ψi
[
∂ρvapor
∂t
+ ur
∂ρvapor
∂r
+ uz
∂ρvapor
∂z
+ ρvapor
ur
r
+ ρvapor
∂ur
∂r
+ ρvapor
∂uz
∂z
− fuente
]
dV
Ri
ρaire
=
∫
Ω
ψi
[
∂ρaire
∂t
+ ur
∂ρaire
∂r
+ uz
∂ρaire
∂z
+ ρaire
ur
r
+ ρaire
∂ur
∂r
+ ρaire
∂uz
∂z
]
dV
(3.39)
Ri
ρresto =
∫
Ω
ψi
[
∂ρresto
∂t
+ ur
∂ρresto
∂r
+ uz
∂ρresto
∂z
+ ρresto
ur
r
+ ρresto
∂ur
∂r
+ ρresto
∂uz
∂z
]
dV
(3.40)
Residuo de la ecuación de conservación de la energía:
Ri
T =
∫
Ω
(φi
[
∂ρtotal
∂t
Cp T + ρtotal
∂Cp
∂t
T + ρtotal Cp
∂T
∂t
+ ur Cp T
∂ρtotal
∂r
+ urρtotal T
∂Cp
∂r
+ ur ρtotal Cp
∂T
∂r
+ uz Cp T
∂ρtotal
∂z
+ uzρtotal T
∂Cp
∂z
+ uz ρtotal Cp
∂T
∂z
+ ρtotal Cp T (
ur
r
+
∂ur
∂r
+
∂uz
∂z
) + fuenteLv)− qr
∂φi
∂r
− qz
∂φi
∂z
]
dV −
∫
Γ
φi q · n̂ dS
︸ ︷︷ ︸
TermA
(3.41)
Residuos de la ecuación de conservación del momento:
Ri
ur =
∫
Ω
(φi[ur
∂ρ
∂t
+ ρtotal
∂ur
∂t
+ 2ρtotalur
∂ur
∂r
+ ρtotal uz
ur
∂z
+ u2
r
∂ρtotal
∂r
+
ρtotal
r
u2
r
+ ur uz
ρtotal
∂z
+ uz ρtotal
∂uz
∂z
+
τθθ
r
− P
r
)− ∂φi
∂r
P + τr r
∂φi
∂r
+ τr z
∂φi
∂z
]dV
−
∫
Γ
n̂ φi((τ · er)− P er)dS︸ ︷︷ ︸
TermB
(3.42)
Víctor Siliang Ruan Zhao 35
Ri
uz =
∫
Ω
(φi[uz
∂ρ
∂t
+ ρtotal
∂uz
∂t
+ ρtotalur
∂uz
∂r
+ 2ρtotal uz
uz
∂z
+ uz ur
∂ρtotal
∂r
+ ur uz
ρtotal
r
+ uz ρtotal
∂ur
∂r
+ u2
z
∂ρtotal
∂z
)− P ∂φ
i
∂r
+ τr r
∂φi
∂z
+ τr z
∂φi
∂r
]dV
−
∫
Γ
n̂ φi((τ · ez)− P ez)dS︸ ︷︷ ︸
TermB
(3.43)
Los términos TermA y TermB de las ecuaciones de la energía y del momento han
sido obtenidos mediante integración por partes de los términos en derivadas parciales
segundas, serán los que introducirán las condiciones de contorno según el planteamiento
de la ecuación (3.34).
Adicionalmente, debido a la imposibilidad de computar las derivadas de las variables
en las fronteras por la discontinuidad de la primera derivada de las funciones base se
introdujo una ecuación adicional que resolvería el tensor gradiente de velocidad:
G = ∇ v (3.44)
Debido a la simetría axial del problema, el tensor gradiente de velocidad tiene 5 com-
ponentes distintas: Grr, Grz, Gzr, Gzz, Gθθ. El desarrollo de su formulación débil se realiza
con funciones lineales como con la ecuación de conservación de la masa. Su resolución es
posterior al problema del flujo representado en la figura 3.8.
3.4.3. Mallado. Despcripción y desarrollo
En las simulaciones numéricas de problemas de mecánica de sólidos o de fluidos es
frecuente la distorsión del dominio que se considera en la interacción con sus fronteras
como paredes o fronteras libres, con fases de distinta densidad o con partículas en el
interior del fluido. Es por tanto muy importante el uso de un método capaz de adaptarse
a las distorsiones y que, idealmente, produzca un mallado de alta calidad durante la
evolución temporal del dominio, adaptando los cambios a el del fluido.
Los dos métodos clásicos utilizados para la descripción del movimiento en problemas
de mecánica de fluidos son: (a) método de Lagrange y (b) método de Euler. El método
utilizado en este trabajo es el llamado Arbitrary Lagrangian-Eulerian (ALE) Method que
combina las ventajas de los dos anteriores a la hora de la descripción del movimiento, [39].
En el método de Lagrange, cada punto del mallado es considerado como un nodo
individual asociado a un punto del material o fluido real durante la simulación, es decir,
el movimiento de un nodo seguirá al movimiento del dominio computacional. Es muy
apropiado para la descripción de fronteras libres sin gran distorsión, dado que si el dominio
se distorsiona este tendría que se reconfigurado (remeshing) en cada paso de tiempo con
el consecuente incremento de tiempo de simulación.
36 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
En el método de Euler el mallado es fijo y el dominio (un fluido) se mueve respecto a
este. Es esperable que las grandes distorsiones sean manejadas con relativas facilidad, ya
que en cada punto y tiempo se identifica la variable simulada con el punto del mallado,
pero con la pérdida de precisión en los cambios de fronteras.
La descripción ALE combina ambas descripciones permitiendo que los nodos del ma-
llado se muevan como en la descripción lagrangiana, en la descripción euleriana o de otra
forma especificada, dando capacidad a los elementos de adaptarse al las zonas más distor-
sionadas. Específicamente, en este trabajo la adaptación del mallado se realiza mediante
métodos elípticos y cuasi-elípticos de generación de mallas que permiten describir el mo-
vimiento de las fronteras libres y geometrías complejas. Los trabajos realizados por [3],
[4] y [5] han demostrado ser muy fiables en problemas axisimétricos con fronteras libres,
con grandes deformaciones y con adaptación del dominio con el tiempo. El dominio del
flujo variable con el tiempo se mapea a un dominio constante (figura 3.4), el dominio de
simulación, sobre el cual se calcula la evolución del mallado mediante las ecuaciones de la
sección 3.4.3.1. El dominio computacional coincide con el dominio físico real para t = 0.
Figura 3.4: Mapeo bidimensional del suflé
Este mapeo al dominio computacional conlleva el uso de la regla de la cadena para el
cálculo de las derivadas de las ecuaciones.
3.4.3.1. Ecuaciones del mallado
Los ecuaciones que conforman el problema del mallado de la figura 3.8 son los pro-
puestas por [3] adaptadas a nuestro problema:
∇
(
ε1
√
R2
ξ + Z2
ξ
R2
η + Z2
η
+ (1− ε1)
)
∇ξ = 0 (3.45)
∇ · ∇η = 0 (3.46)
donde ε1 es una constante entre [0, 1] que se ha tomado como 0.1 en nuestro caso. Los
subíndices ξ y η denotan derivadas parciales respecto de ξ y η respectivamente.
La líneas de coordenadas con valores constantes de η serán paralelas a la frontera libre
Víctor Siliang Ruan Zhao 37
deformada teniendo que cumplir la condición de ser ortogonales, razón por la cual la ecua-
ción (3.46) es una ecuación de Laplace. Las líneas de coordenadas con valores constantes
de ξ también tendrán que seguir una condición de ortogonalidad o cuasi-ortogonalidad
dada por la ecuación (3.45). La solución de ambas ecuaciones nos proporciona las coor-
denadas de (r, z) para cada nodo adaptando de esta manera el mallado en cada paso de
tiempo. Se utiliza de nuevo la formulación débil Galerkin con las funciones base (3.35) y
(3.36) para aproximar (r, z). Se obtienen las siguientes ecuaciones:
Ri
mr =
∫
Γ
n̂
[
ψi
(
∇
(
ξ1
√
R2
ξ + Z2
ξ
R2
η + Z2
η
+ (1− ξ1)
)
∇ξ
)]
dS
−
∫
Ω
[
∇ψi ·
(
∇
(
ξ1
√
R2
ξ + Z2
ξ
R2
η + Z2
η
+ (1− ξ1)
)
∇ξ
)]
dV (3.47)
Ri
mz =
∫
Γ
n̂ψi∇ηdS −
∫
Ω
∇ψi · ∇ηdV (3.48)
Los términos de superficie de estas ecuaciones se pueden ignorar dado que se impondrán
condiciones de contorno en las paredes del molde y en la frontera libre. La solución de
estas ecuaciones con sus correspondientes condiciones de contorno darán las coordenadas
R y Z de todos los puntos del mallado.
3.4.3.2. Frontera libre
Como se muestra en la figura 3.8 el cálculo de las variables del problema del flu-
jo precede al problema del mallado. Resuelto este problema, su evolución obedece a las
condiciones de contorno que han determinado la frontera libre, por imposición de la con-
dición cinemática (3.74), tras lo cual se resuelve el problema del mallado constituido por
la ecuaciones (3.47) y (3.48), teniendo como condición de contorno la frontera libre para
esta última, que determinan la posición para los puntos interiores del mallado.
Los puntos de la frontera libre están definidos por un vector posición r = r(ξ, ϕ),
dependiente de la ξ y ϕ pero no de η ya que son mapeado a una altura constante η =
Height, 3.4.
r = R(ξ, ϕ)er + Z(ξ, ϕ)ez (3.49)
Los vectores tangente a la superficie son, [40]:
aξ =
∂r
∂ξ
aϕ =
∂r
∂ϕ
38 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONESDE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
el vector normal unitario:
n̂ =
aξ × aϕ
|aξ × aϕ|
y el diferencial de superficie:
dS = |∂r
∂ξ
× ∂r
∂η
|dξdη
La frontera libre queda definida por las siguientes ecuaciones:
Primer vector tangente:
aξ =
∂R
∂ξ
er +
∂Z
∂ξ
ez√(
∂R
∂ξ
)2
+
(
∂Z
∂ξ
)2
(3.50)
Segundo vector tangente:
aϕ = eθ (3.51)
Vector normal unitario:
n̂ =
−∂Z
∂ξ
er +
∂R
∂ξ
ez√(
∂R
∂ξ
)2
+
(
∂Z
∂ξ
)2
(3.52)
Diferencial de superficie:
dS = R
√(
∂R
∂ξ
)2
+
(
∂Z
∂ξ
)2
dξdϕ (3.53)
Las expresiones del vector normal y el diferencial de superficie son utilizadas en los
cálculos de los términos de superficie de las ecuaciones que necesitan introducir un término
de frontera y la imposición de una condición natural (o de Neumann), por ejemplo los
términos TermA y TermB de las ecuaciones (3.41), (3.42) y (3.43)
3.5. Condiciones iniciales y de contorno
Se establecieron unas condiciones iniciales para todas las simulaciones:
Presión inicial:
P = 105 Pa (3.54)
Temperatura inicial igual a la temperatura ambiente
T = 25oC (3.55)
Víctor Siliang Ruan Zhao 39
Velocidad inicial nula, el suflé se considera en reposo:
ur = 0 uz = 0 (3.56)
Densidades iniciales de las sustancias del modelo: tabla 6.1. Siguiendo las ecuaciones
detallas en la sección 3.3.
La figura 3.5 muestra la división de la frontera del modelo en 4 partes diferentes: (1)
contacto con el plato inferior, B1; (2) contacto con la pared del molde, B2; (3) la frontera
libre, B3; (4) el eje de simetría, B4. En la figura 3.5 se detalla esta división:
Figura 3.5: Fronteras del esquema bidimensional del suflé
Se describen en esta sección las distintas condiciones de contorno que han sido im-
puestas en los problemas del flujo y en el problema del mallado. Antes se describirán las
fenómenos físicos que describen las condiciones.
3.5.1. Condición de simetría axial. Frontera B4
La primera de las condiciones de contorno y que supone uno de los puntos de partido
del modelo, es la condición de simetría axial de la frontera B4 en la figura 3.5. Su
implementación en los términos Term A y Term B de las ecuaciones (3.41) y (3.42),
(3.43) se realiza según se desarrolla a continuación y teniendo en cuenta que el vector
normal es n̂ = −er y que las derivadas de la temperatura y de la velocidad son nulas
respecto de la componente radial r y longitudinal z ,
∂T
∂r
= 0,
∂ur
∂r
= 0,
∂uz
∂r
= 0 (igual
para z).
40 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
TermB =
∫
Γ
ezφ
i((τ · ez)− P ez)dS =
∫
Γ
ezφ
i((τ · ez)dS =
∫
Γ
φiτrzdS =
∫
Γ
φiη
(
∂uz
∂r
+
∂ur
∂z
)
dS = 0 ⇒ TermB = 0 (3.57)
q · n = −k∂T
∂r
= 0 ⇒ TermA = 0 (3.58)
Por lo que ambos términos son nulos: TermA = TermB = 0 en la frontera B4.
3.5.2. Transferencia de calor. Conducción, convección y radiación
El proceso de cocción del suflé ocurre dentro de un horno tipicamente entre 180 -
200 oC en las recetas de los libros de cocina. Se encuentran presentes los 3 fenómenos
principales de la transmisión de calor: conducción, convección y radiación que supondrán
un flujo de calor de hacia el interior del suflé, q.
El término TermA de la ecuación (3.41) introduce el flujo de calor en las fronteras del
suflé y adaptará una u otra expresión según el fenómeno. El molde en el que se introduce
el suflé al horno no se simula explícitamente, es decir, no existe dentro de la simulación de
elementos finitos pero es necesario su espesor, wm, para el flujo de calor dado que existe
conducción de calor. Se representa el flujo de calor en la figura 3.6 para la parte inferior,
B1 en la figura 3.5 donde suceden los fenómenos de conducción y convección. En la parte
lateral, B2 en la figura 3.5, el mecanismo de transmisión de calor sería el mismo.
Figura 3.6: Representación de la conducción y la convección en el modelo del suflé
En una hipótesis cuasi-estacionaria, para un instante de tiempo dt, el flujo de calor
por conducción se puede considerar constante a lo largo de la dirección X pero no en la
dirección ortogonal a esta dado que los puntos interiores tendrán diferente temperatura.
Esta simplificación de la ecuación es posible debido a que en un instante dt, se considera
la variación de ∂ T/∂ x como constante, siendo la temperatura Tsouffle la temperatura a
determinar por en la ecuación de conservación de la energía (3.41). Se expresa con una
expresión general para la conducción (e.g. [41]):
Víctor Siliang Ruan Zhao 41
qconduction = km
∆T
wm
= km
(Tsurface − Tsouffle)
wm
(3.59)
donde Tsurface y Tsouffle son la temperatura de la superficie del molde y la temperatura
del interior del suflé para el tiempo t, km la conductividad del molde y wm el espesor del
molde.
Material k (W/mK) w(mm)
Alúmina 35 6.7
Tabla 3.8: Propiedades de molde de óxido de aluminio a 300K de [41]
De igual manera para el flujo de calor por convección que ocurre desde el aire dentro
del horno a temperatura Toven = 473,15K(200oC) hasta la superficie del suflé:
qconvection = h (Toven − Tsurface) (3.60)
donde h es el coeficiente de película y Toven es la temperatura del aire interior del
horno. Sustituyendo (3.59) en (3.60) quedaría la expresión:
Tsurface =
h (Toven +
km
wm
Tsouffle)
km
wm
+ h
(3.61)
que volviendo a sustituir en la ecuación (3.59):
q =
(Toven − Tsouffle)
1
h
+
wm
km
(3.62)
q = htotal (Toven − Tsouffle) htotal =
1
h
+
wm
km
−1
(3.63)
Es preferible el uso de la ecuación (3.63) por poder ser introducido el flujo como una
condición de contorno tipo Neumann en el término TermA de la ecuación (3.41) para el
cálculo de la temperatura interior Tsouffle. El coeficiente de película htotal representa, en
un analogía eléctrica, la combinación en serie de los efectos de conducción y convección.
El cálculo de los coeficiente de película se ha realizado siguiendo las correlaciones de
[41] tomando como longitud característica L, la pared del moldeHeight para la correlación
de placa vertical (3.65) y As/p = Rcyl, área (As = π R2
cyl) y perímetro (p = π Rcyl) del
molde, para la correlación horizontal, (3.64):
42 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
Correlación de McAdams para superficie inferior de placa fría, respecto al ambiente:
NuL = 0,54Ra
1/4
L (3.64)
Correlación de Churchill y Chu para placa vertical en el supuesto de flujo laminar:
NuL =
{
0,68 +
0,670Ra
1/4
L
(1 + (0,492/Pr)9/16)4/9
}
(3.65)
donde RaL es el número de Rayleigh, NuL el número de Nusselt (en la dirección
de la longitud característica L), calculadas a la temperatura media de película Tmp =
(Tsurface + Toven)
2
y Pr, número de Prandtl.
Se trata de números adimensionales que se definen como:
RaL =
g β(Tsurface − Toven)L3 ρ2
oven
µ2
oven
NuL =
hL
koven
Pr =
Cp oven µoven
koven
(3.66)
el subíndice oven indica propiedades del fluido, en nuestro caso el aire en el interior
del horno. Los parámetros que se introducen en estas expresiones se recogen en la tabla
3.9 las distintas temperaturas.
Fluido T (K) ρoven ( kg
m3 ) Cp oven ( J
kgK
) koven ( W
mK
) µoven ( kg
ms
) Pr β (K−1)
Aire 298.15 0.89 1013.03 0.033 2.28E-05 0.71 2.51E-03
Aire 473.15 0.83 1016.34 0.034 2.39E-05 0.7 2.36E-03
Aire 410.65 0.83 1016.34 0.034 2.39E-05 0.7 2.36E-03
Tabla 3.9: Propiedades del aire. Las propiedades para la temperatura media de película
Tmp = 410,65K se calculan interpolando entre las anteriores para el instante inicial.
La correlaciones dan como parámetros iniciales para t = 0:
Correlación L (m) RaL Pr NuL h ( W
m2K
)
Placa horizontal (3.64)
As
p
= Rcyl = 0,025 8.06E+04 7.07E-01 9.1 12.19
Placa vertical (3.65) Height = 0,03 1.39E+05 7.07E-01 14.19 15.73
Tabla 3.10: Coeficientes de película iniciales
El valor del coeficiente de película se tomó constante durante las simulaciones dado
que según se puedever en la tabla 3.9 los valores de las propiedades del aire no cambian
significativamente en el rango de temperaturas de trabajo. No obstante, la implementación
Víctor Siliang Ruan Zhao 43
para el cálculo del coeficiente del aire se podría hacer fácilmente de manera sencilla dadas
las propiedades del aire y la temperatura media de película de los puntos de las fronteras.
El otro fenómeno de transmisión que ocurre es la radiación en la frontera libre (B3,
en la figura 3.5), a una temperatura T2 provocada por las resistencias en el interior del
horno a T1. La expresión general puede ser encontrada en muchos libros como [42] y es de
la forma:
qradiation = σ(ε1T
4
1 − αsouffleT 4
souffle) (3.67)
donde la ε1 es la emisividad de las resistencias del horno, α2 la absortividad del suflé
y σ la constante de Stefan-Boltzmann. Debido a la dificultad de cálculo de ε1 y α2, la
expresión anterior se transforma teniendo en cuenta un balance de energía entre ambas
superficies en:
qradiation = σF12(T 4
1 − T 4
souffle) (3.68)
que utiliza el concepto de factor de forma, F12 entre dos superficies que en este caso
correspondería a dos placas planas paralelas. Se han tomado los siguientes valores para
las constantes de emisividad y absortividad de [43] y para el factor de forma de [42].
Constante Nomenclatura valor
Emisividad de las resistencias ε1 0.89
Absortividad del suflé αsouffle 0.8
Factor de forma F12 0.7
Tabla 3.11: Constantes para la ecuaciones (3.67) y (3.68)
3.5.3. Condición de deslizamiento
A la hora de cocinar un suflé se engrasan las paredes del molde con una fina capa
de mantequilla con el objetivo de que la velocidad de esos punto no se anula y el suflé
pueda elevarse. En esta pequeña capa la componente tangencial de la velocidad dependerá
tangencialmente de la componente tangencial del tensor de tensiones que actúa en la
superficie [30]:
(τ · n) · t = −β u · t (3.69)
Esta condición es aplicada solamente a la frontera lateral en contacto con el molde,
B2 en la figura 3.5, donde el vector normal toma la dirección exterior del eje r n = er y
el vector tangente es la dirección positiva del eje z t = ez. De manera que: τrz = −β uz, β
es el coeficiente de deslizamiento.
Este tipo de condiciones de contorno pueden ser fácilmente demostradas consideran-
44 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
do una fina película de fluido, en nuestro caso la mantequilla, cerca de la pared con
propiedades reológicas distintas a las del suflé. La implementación de esta condición, se-
gún la ecuación (3.69), ha sido habitualmente utilizada en simulaciones de mecánica de
fluidos como una manera de evitar las singularidades producidas por el no deslizamien-
to [44], [45]. En esta capa la componente τrz única componente del tensor de tensiones
tendrá la forma τrz = ηbutter
∂uz
∂r
, donde la derivada de la velocidad se expresa tal que
∂uz
∂r
=
0− uz souffle
�
��Rcyl −�
��Rcyl + ε
=
−uz souffle
ε
. Por lo que:
τrz = −ηbutter
ε
uz souffle (3.70)
donde
ηbutter
ε
es la constante de deslizamiento β de la ecuación (3.69). La imposición
de condición se realiza para el residuo de la ecuación (3.43) en el TermB dado que para
la dirección r, ur = 0.
Figura 3.7: Esquema de la capa de deslizamiento
3.5.4. Equilibrio de presiones
En el término TermB de las ecuaciones de los residuos de conservación del momento
(3.42) y (3.43),
TermBk =
∫
Γ
n̂ φi((τ · ek)− P ek)dS k = r, z (3.71)
Víctor Siliang Ruan Zhao 45
se tienen dos términos que deben ser desarrollados: (a) τ , tensor de tensiones y (b) P ,
presión interior del suflé. El trabajo de Zhang y Datta [11] demostró que la contribución
del tensor de tensiones no es muy significativa en el equilibrio y que se puede despreciar
en la simulación. La presión debe mantener un equilibrio con la presión interior del horno
en todo momento, es decir, P = Poven y que fue confirmado en nuestras simulaciones.
La condición de contorno de presión se ha impuesto de manera como una condición de
tipo Neumann en las ecuaciones (3.42) y (3.43) sustituyendo el valor de la presión en la
integral por Poven, presión del aire del horno. Este tipo de implementación ha demostrado
ser más estable como se ha demostrado en [27] y en [46].
3.5.5. Ecuación cinemática
El caso del suflé es también un problema de frontera libre. La parte superior, es una
frontera del dominio inicial que evoluciona, se mueve, con las demás variables del pro-
blema. Su movimiento tiene que ser resuelto mediante una ecuación cinemática adicional
que viene expresado por la derivada material:
D ∗
Dt
=
∂ ∗
∂t
+ v · ∇ ∗ (3.72)
Si r es el vector de posición de los puntos de la frontera libre cuyas coordenadas viven
dadas por: r(r, z) = Rer + Z ez, de la figura 3.4. De manera que tiene que seguir el
movimiento:
Dr
Dt
=
∂r
∂t
+ v · ∇r = v |r=r (3.73)
∂r
∂t
+ v · ∇r =
∂(Rer + Z ez)
∂t
+ v · ∇(Rer + Z ez) (3.74)
Por lo que la velocidad de los puntos de la frontera libre estarán dados por las ecuación:
uz =
∂Z
∂t
+ vr
∂Z
∂r
+ uz
∂Z
∂z
(3.75)
vr =
∂R
∂t
+ vr
∂R
∂r
+ uz
∂R
∂z
(3.76)
siendo ambas ecuaciones equivalentes. Utilizando la regla de la cadena para las deriva-
das parciales y teniendo en cuenta que la frontera libre es mapeado a una altura constante
η = Height, fig. 3.4, en el dominio computacional la expresión (3.75) o (3.76), la condición
cinemática es de la forma:
ur
∂Z
∂ξ
− ∂Z
∂ξ
∂R
∂τ
+
∂Z
∂τ
∂R
∂ξ
− uz
∂R
∂ξ
= 0 (3.77)
46 Escuela Técnica Superior de Ingenieros Industriales (UPM)
GASTRONOMÍA MOLECULAR: ESTUDIO DE CONDICIONES DE CONTORNO
DE UN MODELO TERMOMECÁNICO BIDIMENSIONAL DEL SUFLÉ
Aunque esta condición de contorno implique variables del problema de flujo, debe ser
aplicada al problema del mallado en la ecuación (3.48) dado que está relacionado con la
derivación de de la frontera libre por lo que se consiguen los nuevos puntos nodales. De
esta manera la información del problema del flujo es transmitida al problema del mallado
acoplando las ecuaciones.
3.5.6. Resumen de las condiciones de contorno
Explicadas las posibles condiciones de contorno realista aplicables las específicas pa-
ra cada frontera y para cada problema (flujo o mallado) serán (ver figura 3.5 para la
numeración de las fronteras):
1. Plato inferior, B1:
a) Para el problema del flujo:
1) Calor transmitido por conducción y convección según la ecuación (3.63)
y la correlación (3.64) para el coeficiente de película total de los efectos
de conducción y convección. Los valores de los coeficientes de película se
encuentran en la tabla 3.10.
2) Los puntos de esta frontera adquieren la velocidad de la pared y su velo-
cidad es nula ur = uz = 0, como una condición de no deslizamiento.
b) Para el problema del mallado:
1) Se impone que no evolucione la coordenada, Z = 0, y que los nodos estén
equidistribuidos a lo largo de R.
2. Pared del molde, B2:
a) Para el problema del flujo tendremos:
1) Transmisión de calor por conducción y convección según la ecuación (3.63)
y la correlación (3.65) para el coeficiente de película total de los efectos
de conducción y convección. Los valores de los coeficientes de película se
encuentran en la tabla 3.10.
2) La condición de deslizamiento (3.70) para que los nodos tenga una veloci-
dad uz; aplicada a la ecuación (3.43). Velocidad radial ur = 0.
b) Para el problema del mallado:
1) Equidistribución de los nodos en Z y que R = Rcyl.
3. Frontera libre, B3:
a) Para el problema del flujo tendremos:
1) Flujo de calor por radiación dado por la ecuación (3.68)
2) Equilibrio de presiones, P = Poven impuesta en (3.71) como parte de las
ecuaciones (3.42) y (3.43).
b) Para el problema del mallado:
1) Equidistribución de los nodos para la coordenada R.
Víctor Siliang Ruan Zhao 47
2) Condición cinemática (3.73) de la frontera libre a la que aplicada derivación
por medio de la regla de la cadena se obtiene la ecuación (3.77) dará