ES2389098T3 - Método para describir la distribución de la orientación estadística de partículas en una simulación de un proceso de llenado de molde - Google Patents

Abstract

Un método para calcular las estadísticas de la orientación de partículas no esféricas a nivel macroscópico mediante el uso de un modelo de simulación para simular un proceso de moldeo por inyección en el que una cavidad de molde, la que forma al menos parte de la geometría de un dominio de simulación, se llena o parcialmente se llena con una suspensión que se forma por un solvente que contiene un gran número de partículas no esféricas, 5 en donde se proporciona una representación digital o un modelo de computadora de la geometría del dominio de simulación, y en donde se forma una malla con una pluralidad de celdas de cálculo mediante la subdivisión o discretización al menos de una parte del dominio de simulación, el método que comprende: a) especificar las condiciones de frontera; b) establecer las condiciones iniciales; c) resolver las ecuaciones de balance de masa, el momento y la energía de al menos una porción de las celdas del dominio de simulación para obtener un flujo de fluido, un flujo de calor y una transferencia de masa a nivel macroscópico; y d) resolver las ecuaciones dinámicas de orientación de partículas no esféricas basado al menos parcialmente en los resultados de las ecuaciones de balance resueltas para de esta manera determinar los cambios en la orientación de las partículas no esféricas a nivel macroscópico como función del espacio y del tiempo, en donde la orientación de las partículas no esféricas se describe estadísticamente por una función de distribución y en donde la función de distribución se aproxima por tensores de orientación de partículas no esféricas de orden creciente y el tensor de orientación de 4to orden se calcula como una función del tensor de 2do orden usando una aproximación de clausura 20 híbrida estabilizada.

Description

Método para describir la distribución de la orientación estadística de partículas en una simulación de un proceso de llenado de molde.

La presente invención se relaciona con el campo del modelado tridimensional de un flujo de una suspensión que contiene partículas en la cavidad de un molde, y más específicamente con un método y aparato para describir la distribución estadística de la orientación de partículas no esféricas en una simulación de un proceso en donde la cavidad de un molde se llena con una suspensión que contiene un gran número de partículas no esféricas.

Arte anterior

Una simulación tridimensional exacta de un proceso de moldeo por inyección o un proceso de fundición de metales involucra un sistema de hasta varios cientos de miles de ecuaciones. Se han hecho progresos en el pasado para mejorar la eficiencia de los métodos de simulación para hacer frente a estos cálculos complejos. Con el software optimizado y la capacidad de procesamiento de las estaciones de trabajo modernas tales simulaciones pueden realizarse en un puesto de trabajo, es decir los resultados se obtienen lo suficientemente rápido para ser adecuados fuera del área de investigación puramente científica y puede aplicarse por ingenieros en los departamentos de investigación y desarrollo, en las fundiciones y por los fabricantes de artículos moldeados por inyección.

Cuando unas partículas, tales como fibras, se añaden a la composición de un polímero, y la distribución de la orientación de las fibras necesita describirse, la simulación y el conjunto de ecuaciones de pertenencia se hacen significativamente más complejos. La simulación tridimensional exacta de tal proceso no se ha introducido exitosamente hasta ahora en el puesto de trabajo dado que la complejidad de la simulación no permitiría resultados aceptables en las estaciones de trabajo actuales, dado que el tiempo de cálculo sería muy largo, o la precisión de la simulación sería inadecuada.

En las piezas reforzadas con fibras es frecuentemente crucial para los ingenieros de desarrollo tener una descripción de la distribución de la orientación de las fibras para poder predecir aspectos sobre la tensión y la deformación del componente. Típicamente las fibras se usan para mejorar las propiedades mecánicas de piezas plásticas. Pero entonces las propiedades (termo-) mecánicas (como la expansión térmica, resistencia y rigidez) dependen de la orientación de las fibras.

El uso de componentes de plástico moldeado por inyección se ha incrementado sostenidamente en muchas industrias en los últimos años. Los fabricantes de equipos electrónicos, bienes de consumo, equipos médicos, y piezas de automóviles producen cada vez más y más de sus productos y componentes usados en su mercadería de plástico que nunca antes.

Las piezas moldeadas por inyección reforzadas con fibras sustituyen componentes metálicos estructurales debido a que ofrecen una mejor razón de resistencia/peso, durabilidad, integración de componentes y bajos costos.

Al mismo tiempo, las presiones competitivas llevan a los fabricantes en la industria de plásticos moldeados por inyección a encontrar nuevos métodos para optimizar los diseños con el objetivo de adaptar mejor los diseños al proceso de producción. Cuando se descubre tarde que existe necesidad de hacer modificaciones de un componente

o de la configuración de un molde en el proceso de desarrollo del diseño, el retraso y los costos asociados para implementar los cambios necesarios son significativamente mayores que en etapas tempranas de las etapas de desarrollo del diseño.

Las empresas que desean asegurar que sus componentes sean producibles y funcionen óptimamente les gustaría usar técnicas de ingeniería asistida por computadora para simular o modelar los flujos complejos y la orientación resultante de las fibras en un molde de inyección, con el objetivo de comprender mejor el proceso de fabricación e integrar este conocimiento en el diseño de componentes, tempranamente en la fase de diseño.

Hay un número de factores que deberían considerarse cuando se diseña un molde de inyección para un componente reforzado con fibras y para el componente reforzado con fibras que se producirá en el mismo. Los parámetros tales como la geometría general del componente, el grosor mínimo y máximo de las paredes, el número y la localización de las salidas en el molde a través de las que se inyecta el polímero líquido y la suspensión de fibras, el número y la localización de los respiraderos en el molde a través de los que se escapa el gas de la cavidad, la composición y las propiedades del polímero, las propiedades y la cantidad de fibra, la contracción, las tolerancias y la distribución de la orientación de las fibras son unos pocos. Debido a la estrecha relación entre sí, el diseño del componente y el molde no puede basarse confiablemente solo en la forma y función del componente final, sino también debería considerar los efectos del proceso de fabricación.

La simulación en ingeniería asistida por computadora puede usarse ventajosamente para proporcionar a los ingenieros de diseño y fabricación retroalimentación visual y numérica de lo que es probable que suceda dentro de la cavidad del molde durante el proceso de moldeo por inyección, que les permite comprender mejor y predecir el comportamiento de los diseños contemplados de los componentes de manera que se pueda eliminar sustancialmente la tradicional, costosa estrategia de prueba y error para la fabricación. El uso de la simulación en ingeniería asistida por computadora facilita la optimización de los diseños de componentes, los diseños de moldes y de los parámetros de procesamiento de fabricación durante la fase de diseño, donde los cambios necesarios puedan implementarse fácilmente, con el menor costo e impacto en la planificación.

La aplicación de técnicas de simulación CAE dentro del proceso de ingeniería de componentes reforzados con fibras comprende (i) una simulación de un proceso de fabricación de “moldeo por inyección” que incluye el cálculo del flujo de fluidos y de la transferencia de calor y (ii) los cálculos de esfuerzo & resistencia (y posiblemente durabilidad), todos realizados a nivel macroscópico, para estos componentes para determinar sus propiedades mecánicas funcionales bajo cargas externas. Ambos tipos de simulación requieren modelos adecuados de materiales que describen las propiedades del material polímero que contiene las fibras sumergidas en estado líquido así como en estado sólido.

La escala de longitud en el nivel macroscópico se determina por las dimensiones lineales (tamaño total, grosor de paredes, etc.) de la geometría del componente que típicamente varía en el intervalo de unos pocos mm hasta cm. Las dimensiones de las celdas de cálculo tienen que resolver las escalas macroscópicas de longitud con suficiente precisión, por lo tanto ellas son típicamente de hasta un orden de magnitud menor que la menor dimensión macroscópica. Cuando las dimensiones típicas de las fibras sumergidas en piezas reforzadas con fibras cortas son uno o dos órdenes de magnitud por debajo de las dimensiones típicas de las celdas de cálculo macroscópicas, las propiedades de la fibra, las que son relevantes para el modelado del comportamiento del material macroscópico se describen mediante un enfoque estadístico. Para materiales reforzados con fibras cortas las propiedades macroscópicas relevantes son: (a) la concentración volumétrica, la que es típicamente (aproximadamente) constante en toda la pieza, y (b) la distribución local de la orientación de las fibras (FO) dentro de cada celda de cálculo, la que típicamente varía significativamente a través de la geometría de la pieza. (Detalles adicionales de este tema se discuten en las secciones 1.1 y 1.2 de la descripción detallada)

Una descripción simplificada y apropiada para propósitos prácticos de la distribución estadística de la FO local se proporciona por medio de los momentos de bajo orden (es decir 2do y 4to) de la función de distribución correspondiente. Debido a su estructura matemática estos momentos se denotan como tensores de orientación (de 2do y 4to orden respectivamente). Dentro del marco de simulación CAE para componentes reforzados con fibras los tensores de 4to orden se necesitan para predecir las propiedades reológicas, así como mecánicas del material reforzado con fibras en el nivel macroscópico, dado que estas son propiedades del tensor de 4to orden. El tensor de FO de 2do orden es una matriz simétrica 3 × 3 de valores reales que tiene traza unitaria, por lo tanto sólo 5 de sus 9 componentes son independientes. El número de componentes independientes del tensor de FO de 4to orden se reduce de 34 = 81 a 15 por medio de la simetría (total).

El modelo matemático que describe la distribución de la FO en términos de sus momentos se simplifica significativamente por medio de un cálculo aproximado del tensor de orientación de 4to orden en términos del tensor de 2do orden en términos de una relación de clausura. Una relación de clausura proporciona la descripción matemática de tal esquema de cálculo en términos de una función, y el procedimiento de cálculo relacionado se

denota como “aproximación de clausura” si la relación de clausura es sólo aproximadamente válida bajo

suposiciones específicas. El enfoque de usar sólo el tensor de FO de 2do orden junto con cierta aproximación de clausura conduce a un modelo de tipo “Folgar-Tucker” para simular la evolución del tensor de FO de 2do orden en tiempo y espacio durante el proceso de llenado del molde (ver secciones 2.2 a 2.5 en la descripción detallada para más detalles).

El documento “Glass fibre orientation within injection moulded automotive pedal - Simulation and experimental studies, B.R. Whiteside y otros, Plastics, Rubber and Composites, 2000, Volumen 29, No. 1” describe un método

para modelar la distribución de la orientación de las fibras dentro de un artículo termoplástico reforzado usando el análisis y el flujo asimétrico del termoplástico que contiene un algoritmo de predicción de la orientación de las fibras. El software aproxima un modelo tridimensional usando una malla bidimensional de elementos finitos que consiste de elementos triangulares lineales. Los campos de flujo se calculan usando la aproximación generalizada de Hele-Shaw y una variación de la ecuación de Folgar-Tucker para calcular la orientación de las fibras. Los cálculos de la orientación de las fibras, la temperatura y la viscosidad se realizan usando una técnica de diferencias finitas sobre

19 láminas a través del “grosor” de cada elemento con el objetivo de producir una solución tridimensional. Sin

embargo, es importante señalar que este sistema pudiera no describirse como verdaderamente tridimensional porque el modelo no puede simular componentes de velocidad en la dirección hacia fuera del plano (una limitación de la aproximación de Hele-Shaw). El método descrito en este documento, cuando se adapta a una simulación tridimensional exacta, resultaría en una simulación inestable y consumidora de mucha potencia de cálculo que pudiera no usarse en el puesto de trabajo de por ejemplo ingenieros de desarrollo.

En WEN-HSIEN YANG, RONG-YUE CHANG Y OTROS: “Computer simulation of 3d short fiber orientation in injection molding” ANTEC 2003 PLASTCS: ANNUAL TECHNICAL CONFERENCE, VOLUMEN 1: XP-002501225 se

describe un enfoque numérico tridimensional exacto para simular el llenado de colada y la orientación de las fibras en piezas de geometría compleja moldeadas por inyección. Un método de volumen finito eficiente se combina con el método VO para resolver el flujo de colada durante el moldeo. La distribución de la orientación de las fibras se describe por los tensores de orientación de segundo orden. La orientación tridimensional de las fibras se determina mediante la ecuación de Folgar-Tucker, en la que se introduce un coeficiente de interacción para tomar en cuenta la interacción fibra-fibra. Varios casos ejemplares se simulan para verificar el avance predicho de frente de colada y la correspondiente orientación de las fibras.

Descripción de la invención

En este contexto, es un objeto de la presente invención proporcionar un método para determinar la distribución de la orientación de partículas no esféricas en una simulación de un proceso en donde una cavidad de un molde se llena con una suspensión que contiene un gran número de partículas no esféricas que es más estable y menos intensivo en cálculos que los métodos del arte anterior.

Este objeto se logra de acuerdo con la reivindicación 1 proporcionando un método implementado en computadora para calcular las estadísticas de la orientación de partículas no esféricas a nivel macroscópico mediante el uso de un modelo de simulación para simular un proceso de moldeo por inyección en el que una cavidad de un molde, la que forma al menos parte de la geometría de un dominio de simulación, se llena o se llena parcialmente con una suspensión que se forma por un solvente que contiene un gran número de partículas no esféricas, en donde se proporciona una representación digital o modelo de computadora de la geometría del dominio de simulación, y en donde se forma una malla con una pluralidad de celdas de cálculo mediante la subdivisión o discretización de al menos una parte del dominio de simulación, el método comprende los pasos de a) especificar las condiciones de frontera; b) establecer las condiciones iniciales; c) resolver las ecuaciones de balance de masa, de momentos y de energía para al menos una porción de las celdas del dominio de simulación para obtener el flujo de fluido, el flujo de calor y la transferencia de masa a un nivel macroscópico; y d) resolver las ecuaciones dinámicas de la orientación de partículas no esféricas basado al menos parcialmente en los resultados de las ecuaciones de balance resueltas para de esta manera determinar los cambios en la orientación de las partículas no esféricas a nivel macroscópico como función del espacio y el tiempo. Aquí, para el paso d), las ecuaciones de orientación de las partículas sólo pueden resolverse para celdas de cálculo que contienen suspensión.

Se prefiere que el paso c) comprenda además (cc) determinar una superficie libre actualizada o frente de flujo del fluido o suspensión, cuya superficie libre separa las celdas llenas con la suspensión de las celdas vacías de la cavidad, basada al menos parcialmente en los resultados de las ecuaciones de balance resueltas. Preferiblemente, el paso (cc) comprende además actualizar las condiciones de frontera de acuerdo con el frente de flujo actualizado. También se prefiere que el método de la invención comprenda además los pasos e) determinar si el proceso de moldeo por inyección simulado ha finalizado mediante la determinación de si la cavidad del molde está llena por la inyección simulada de la suspensión; y f) repetir los pasos c), cc) y d) hasta que finalice el proceso de moldeo por inyección simulado.

Este objeto de la presente invención también se logra de acuerdo con la reivindicación 6 proporcionando un método para describir la distribución de la orientación estadística de partículas no esféricas en una simulación de un proceso en donde una cavidad de un molde se llena con una suspensión que se forma por un solvente que contiene un gran número de partículas no esféricas, el método comprende (1) proporcionar un modelo tridimensional por computadora que defina la geometría de la cavidad; (2) especificar las condiciones de frontera; (3) discretizar un dominio de solución basado en el modelo para formar una malla con una pluralidad de celdas; (4) resolver las ecuaciones de energía y flujo para al menos una porción del dominio de solución; (5) calcular las condiciones de flujo y de temperatura en las respectivas celdas como una función del tiempo; (6) calcular los cambios en la orientación de partículas no esféricas; y (7) describir la distribución estadística de la orientación de las partículas no esféricas en las celdas respectivas como una función del tiempo.

Los métodos de la invención utilizan el flujo de una suspensión para predecir la distribución de la orientación de las fibras y las propiedades termomecánicas a través de la pieza con un esfuerzo computacional sustancialmente reducido. Los resultados de la simulación pueden usarse por un ingeniero de desarrollo, quien puede de esta manera optimizar los productos en relación con la orientación de las fibras y por lo tanto mejorar la resistencia y la estabilidad de la forma del artículo.

La presente invención se relaciona también con un producto de software para computadora, como se define en la reivindicación 16, almacenado en un medio legible por computadora, el software que incluye un código de software para ejecutar el método de acuerdo la reivindicación 1.

Las modalidades preferidas de la invención se definen en las reivindicaciones dependientes.

Los objetos, las características, las ventajas y las propiedades adicionales del método y el aparato para determinar la distribución de la orientación de las partículas no esféricas en una simulación de un proceso en donde una cavidad de molde se llena con una suspensión que contiene un gran número de partículas no esféricas de acuerdo con la invención se harán evidentes a partir de la descripción detallada.

Breve descripción de los dibujos

En la siguiente parte detallada de la presente descripción, la invención se explicará en más detalles con referencia a las modalidades ejemplares mostradas en los dibujos, en los que:

La Fig. 1 es una vista en sección transversal a través de una representación esquemática de una máquina de moldeo por inyección que incluye un molde, La Fig. 2 es un diagrama de flujo de nivel superior que resume los pasos básicos de proceso de un proceso de moldeo por inyección de acuerdo con una modalidad de la invención, La Fig. 3 es un diagrama de flujo que ilustra con más detalles el paso 5 del diagrama de flujo de la Fig. 2; La Fig. 4 es una micrografía de un artículo plástico reforzado con fibras; La Fig. 5 ilustra un modelo de una fibra usado en una modalidad de la invención; La Fig. 6 es un gráfico que ilustra los términos cuadrático & cúbico y lineal & constante del polinomio característico; La Fig. 7 es un gráfico que ilustra la terna de valores propios de la matriz; La Fig. 8 muestra un ejemplo especial correspondiente a un caso “máximamente simétrico”; La Fig. 9 muestra versiones distorsionadas del volumen en forma de tetraedro; La Fig. 10 muestra un dibujo de la estructura general del conjunto del espacio de fase MFT de las matrices de FO (es decir los tensores de FO de 2do orden) de acuerdo con una modalidad de la invención; La Fig. 11 es un diagrama de flujo que ilustra que el proceso de división de operador está en una forma simplificada de acuerdo con una modalidad de la invención; La Fig. 12 es un diagrama de flujo que ilustra los detalles del método de espaciado del tiempo; y La Fig. 13 es un diagrama de flujo que ilustra un proceso de proyección del espacio de fase que usa reescalado de traza de acuerdo con una modalidad de la invención.

Descripción detallada

La Fig. 1 muestra esquemáticamente una máquina de moldeo por inyección 1. La máquina de moldeo por inyección se proporciona con un tornillo 2 que se alimenta con gránulos de polímero que se disponen en una tolva 3. Los gránulos de polímero se transforman por la acción del tornillo 2 y los elementos de calentamiento 4 en una masa viscosa que se empuja a alta presión en una cavidad de molde 5 en el molde 6. Y la máquina de moldeo y el ciclo de fabricación de moldeo por inyección son bien conocidos en el arte y no se explican en detalle aquí. Con la máquina de moldeo 1 se pueden fabricar piezas plásticas tanto no reforzadas con fibra como reforzadas con fibra.

La simulación numérica del proceso de moldeo por inyección en una computadora puede llevarse a cabo de conformidad con el proceso ilustrado en la Fig. 2.

Los pasos principales identificados de una simulación generalmente son los siguientes:

-

paso 1, proporcionar una representación digital de la geometría del dominio de simulación;

-

paso 2, enmallar, lo que es una subdivisión del dominio de cálculo en muchos elementos pequeños, los que son las bases para discretizar las ecuaciones diferenciales (utilizando diferentes algoritmos de solución) y de esta manera encontrar las soluciones a los fenómenos físicos a ser simulados;

-

paso 3, anexar los datos físicos necesarios a los diferentes dominios de materiales dentro del modelo de simulación (base de datos o banco de datos);

-

paso 4, especificar las condiciones de frontera para el proyecto de simulación,

-

paso 5, simular el proceso de moldeo por inyección (este paso será una reivindicación en mayor detalle más adelante); y

-

paso 6, mostrar los resultados como una presentación gráfica o numérica en la pantalla de una computadora tal como una estación de trabajo.

Los detalles del paso cinco cuando se simula el proceso de moldeo por inyección de un artículo reforzado con fibras se ilustran en el diagrama de flujo de la Fig. 3. En esta parte del proceso las ecuaciones diferenciales para el flujo de calor, el flujo de fluido y los esfuerzos y las tracciones se resuelven usando algoritmos numéricos:

-

en el paso 1 se establecen las condiciones iniciales de las propiedades termofísicas del material y el frente de flujo;

-

en el paso 2 se resuelven las ecuaciones térmicas para todo el dominio y las ecuaciones de flujo en todas las celdas del fluido usando las ecuaciones de conservación de masa, de energía y de momentos;

-

paso 3, en este paso el frente de flujo se mueve y las condiciones de frontera se adoptan de acuerdo con el nuevo frente de flujo;

-

en el paso 4 la orientación y el transporte de las fibras se calcula a partir de la velocidad de flujo obtenida en los pasos anteriores. Este paso se explicará con mayor detalle más adelante. Las condiciones iniciales se aplican a las celdas recién llenas. Sólo se consideran las celdas que contienen material fluido;

-

en el paso 5 se calculan magnitudes adicionales como las reacciones químicas, y se verifica si las celdas se solidifican;

-

en el paso 6 se verifica si el proceso de moldeo por inyección en molde ha finalizado; si no la simulación continúa con el próximo paso de tiempo y el proceso regresa al paso 2;

-

en el paso 7 se calculan las propiedades del artículo después de la extracción del molde.

Las temperaturas del material termofísico fuera del molde es decir la temperatura, contracción, torcedura, etc. se calculan usando la información obtenida a partir de la simulación del moldeo por inyección.

La Fig. 4 es una micrografía de un artículo de plástico reforzado con fibras, en el que puede verse la orientación de las fibras después que el proceso de moldeo por inyección ha finalizado. La orientación de las fibras en el producto final depende en gran medida del patrón de flujo de la masa termoplástica durante el proceso de moldeo por inyección. La orientación de las fibras no se determina exactamente para cada fibra individual, sino más bien descrita por una función de distribución.

Antes de entrar en detalles sobre el cálculo de la orientación de las fibras, se proporciona una revisión del aspecto básico de la mecánica de fluidos de la fundición de un termoplástico con fibras cortas.

1. Aspectos básicos de la mecánica de fluidos de la fundición de un termoplástico con fibras cortas

Una masa de polímero incluye un gran número de fibras cortas dispersas en la misma de manera que el componente producido será fabricado de un termoplástico reforzado con fibras cortas. En el estado fundido, es decir cuando la temperatura es suficientemente alta de manera que la matriz termoplástica es líquida, la mezcla de la masa fundida de plástico y las fibras dispersas constituye un fluido complejo el que comúnmente se denota como una suspensión de partículas en la terminología de la dinámica de fluidos y la física. En general tal suspensión se compone de dos fases diferentes: (i) el solvente, el que en nuestro caso se corresponde con el material plástico fundido sin las fibras dispersas, y (ii) la fase de partículas, la que consiste de todas las fibras sumergidas en el solvente.

1.1 Geometría y propiedades materiales de las fibras

Las partículas esféricas que poseen un eje de simetría rotacional (como se muestra en la Fig. 5) de aquí en lo adelante se denotan como fibras, y los términos partícula y fibra se usan como sinónimos excepto donde se indique explícitamente lo contrario. La geometría de las fibras se caracteriza por su longitud ℓ y su diámetro d así como también la relación de aspecto ra = ℓ/d resultante de ambas cantidades. Los valores típicos de estas cantidades para fibras (por ejemplo de tipo de carbono o de vidrio) dispersas en materiales termoplásticos reforzados con fibras cortas son por ejemplo ℓ ≈ 0.5 mm, d ≈ 5 pm y ra ≈ 100. Típicamente estos valores varían en el rango de ℓ = 0.1 hasta 1 mm, d = 1 hasta 10 µm con los correspondientemente grandes valores de relación de aspecto.

Para la mayoría de los tipos de fibras sumergidas en plásticos reforzados con fibras, la densidad de masa del material de las fibras se compara con la densidad de la masa fundida de plástico que las suspende. El plástico fundido en sí tiene una más bien alta viscosidad a las temperaturas típicas que ocurren durante la fase de llenado del molde. Mediante la combinación de estos dos aspectos los efectos de flotabilidad así como los efectos inerciales se desprecian cuando se describe el movimiento de las fibras dentro de la masa fundida de plástico en el método de acuerdo con la modalidad preferida.

Concerniente a la concentración de fibras, se asume una distribución espacial homogénea de las fibras en toda la suspensión y por lo tanto se toma que la concentración es constante. Los efectos que pueden alterar pueden ser (i) un transporte puramente convectivo de una distribución no homogénea ya presente en la entrada o (ii) la ocurrencia de efectos de migración de partículas inducidos por cizallamiento en el curso de la fase de llenado del molde.

Si bien (i) o (ii) constituyen un efecto significativo, el modelo usado para describir el flujo de la suspensión de las fibras cortas sumergidas en una masa fundida termoplástica necesariamente tiene que ser un modelo de “flujo de dos fases” que acopla mutuamente el flujo de las partículas y la fase de solvente y permite una concentración de

partículas no homogénea. Sin embargo, este no es el caso. El fenómeno (ii) puede comprenderse como un tipo de proceso de difusión, donde la constante de difusión se escala como el cuadrado del tamaño relativo de la partícula, el que se define como la razón entre el tamaño real de las partículas suspendidas (es decir la longitud de las fibras,

o el radio en el caso de partículas esféricas) a una dimensión típica de la cavidad de flujo (por ejemplo el grosor de la pared). Como el tamaño relativo de la partícula es típicamente del orden de 0.1 a 0.01 incluso para piezas de paredes delgadas fabricadas mediante moldeo por inyección de plásticos reforzados con fibras, es perfectamente posible que (ii) sea completamente despreciable en este caso. Concerniente a (i) pudiera imaginarse que pudiera haber un perfil de concentración no homogénea en la entrada debido al proceso de preparación de la masa fundida en el tornillo. Sin embargo, la concentración de las fibras experimentalmente observadas en la práctica para las piezas plásticas reforzadas con fibras cortas típicas muestra solo desviaciones muy pequeñas de un valor constante casi donde quiera dentro de las piezas, lo que es quizás el mejor argumento para asumir una distribución homogénea de fibras [20].

1.2 Influencia de las fibras sobre las propiedades reológicas

Los materiales termoplásticos muestran un comportamiento de flujo no newtoniano. En el caso de una masa fundida de plástico pura sin fibras, las propiedades reológicas del material pueden modelarse por una función escalar de viscosidad que depende de la temperatura así como de las variables de estado del flujo (fluidos newtonianos generalizados). Aunque se conoce que sólo cubren un rango limitado de las propiedades de los flujos no newtonianos (como por ejemplo la fluidificación por cizallamiento), han demostrado ser útiles y satisfactoriamente precisos para los propósitos de la simulación de moldeo por inyección. Los modelos de este tipo se usan también en el método de acuerdo con la modalidad preferida de la invención.

La adición de fibras al material matriz termoplástico le da propiedades mecánicas las que en el estado sólido son fuertemente anisotrópicas y dependen en gran medida de la distribución local de las direcciones en las que están apuntando las fibras embebidas. En principio el comportamiento del material anisotrópico está presente también si el material está en estado líquido (es decir fundido). Para tener en cuenta esta anisotropía se tendría que sustituir la viscosidad escalar mencionada anteriormente por un tensor de viscosidad. Se ha hecho un número de investigaciones que comparan ambas vías de modelar el comportamiento material de las suspensiones de fibras. Para la mayoría de las situaciones de flujo como las encontradas en el llenado de un molde se ha encontrado que el patrón de llenado predicho por ambos tipos de modelos (es decir el escalar y el tensorial) muestra poca diferencia (ver por ejemplo [21], [22]). Por lo tanto los efectos de viscosidad anisotrópicos se desprecian y se usa un simple modelo newtoniano generalizado aquí también cuando la masa fundida termoplástica contiene fibras. Esto es equivalente a despreciar la influencia de la orientación de las fibras dispersas sobre las propiedades reológicas de la masa fundida. La viscosidad depende de la concentración de fibras, pero como esta se supone constante (ver más arriba), este aspecto entra sólo como un parámetro material conocido a priori el que contribuye a las propiedades materiales efectivas de la masa fundida termoplástica.

Debido a este enfoque el método de acuerdo con la modalidad preferida usa un desacople parcial de los cálculos del flujo y los de la orientación de las fibras: mientras que la velocidad de flujo influye localmente en la orientación de las fibras dispersas, hay una influencia despreciable de la orientación de las fibras en el flujo. Por lo tanto el cálculo del flujo se lleva a cabo independientemente del cálculo de la orientación de las fibras. De esta manera la velocidad de flujo local de la masa fundida entra en el modelo usado para el cálculo de la orientación de las fibras como un

conjunto de coeficientes externos.

2. El modelo de Folgar-Tucker 2.1 La ecuación de Jeffery

Aunque las fibras suspendidas en la masa fundida del polímero son partículas delgadas en el sentido de que su relación de aspecto ra = ℓ/d es grande, ellas son lo suficientemente cortas de manera que las fuerzas mecánicas que actúan sobre ellas en el campo de flujo local del solvente no son capaces de provocar ninguna deformación sustancial. Por lo tanto una fibra individual es aquí modelada como un cuerpo rígido simétrico rotacionalmente, delgado cuya orientación se da por un vector unitario p dirigido a lo largo de los ejes de simetría rotacional. Los dos vectores ±p ambos representan el mismo estado de orientación de la partícula.

La Fig. 5 ilustra una partícula rígida de forma elíptica rotacionalmente cuyo movimiento se afecta por el campo del vector de velocidad U en la vecindad de la partícula. En adición a las magnitudes (longitud ℓ y diámetro d) que caracterizan a la geometría de la partícula, se muestran el vector de orientación p. La dirección y el tamaño de la velocidad de flujo U se indica por la dirección y la longitud de las flechas correspondientes. Aunque las direcciones de los vectores de velocidad son todas las mismas, sus longitudes no lo son, lo que indica que la velocidad de flujo en la vecindad de la partícula no es constante. Sin embargo, la cantidad local de cambio de la velocidad de flujo, es decir el gradiente de velocidad, parece ser el mismo, ya que hay un aumento constante de abajo hacia arriba en la longitud de los vectores. En cada punto de la superficie de la partícula, se supone que la partícula se mueve exactamente con la velocidad de flujo local (condición de frontera “sin deslizamiento”). Si la velocidad de flujo fuese constante alrededor de la partícula, esto resultaría en un simple movimiento de traslación, es decir un transporte convectivo puro. De lo contrario, en la presencia de un gradiente de velocidad, la partícula también realiza un movimiento de rotación, como se indica por las flechas de trazos. En total, el tipo más general de movimiento

realizado por la partícula rígida es una traslación con la “velocidad promedio” alrededor de la partícula, combinada

con una rotación alrededor de su centro de masa que se impulsa por el gradiente de velocidad que describe la desviación local de la velocidad de flujo de su valor promedio en la vecindad de la partícula.

La descripción cualitativa dada anteriormente se ha formulado en un modelo matemático por Jeffery [2], asumiendo que las fuerzas viscosas son dominantes, las fuerzas inerciales por lo tanto despreciables y la variación del flujo local es pequeña sobre la región barrida por una sola fibra durante su movimiento. Como se explicó anteriormente, todas estas suposiciones se cumplen en el caso de las fibras cortas dispersas en una masa fundida termoplástica. La asumida pequeñez de la variación de la velocidad de flujo alrededor de la fibra implica que el tensor de gradiente

U essuficiente para describir estavariación con precisión. Estetensor es una matriz 3×3
ٔ׏ local de velocidad

cuyos elementos se calculan a partir de las derivadas parciales de las componentes Uj del vector de velocidad de

=ijU)
ٔ׏ en el espacio, es decir: flujo con respecto a las coordenadas cartesianasxide un punto r = (x1, x2, x3)T

∂Uj/∂xi. De manera que si r0 es la coordenada espacial del centro de masa de la partícula, los valores de la velocidad de flujo U en los puntos r cerca de la partícula y el tiempo t son bien aproximados por un desarrollo de Taylor de

)conerror despreciable, y el gradiente de velocidad puede r0,1))T · (r -r0U(
ٔ׏ ,t)+(r0t) ≈ U(,rprimer orden U() set,0r(U
ٔ ׏ , t) ≈r(U
ٔ ׏ , es decirpuede asumirsequelocalmente constante considerarse como una magnitud

cumple para todos los puntos r cerca de r0 .

En este caso, el estado de orientación transitoria de una fibra, el que se da por su vector de orientación p(r, t) como una función de las coordenadas del espacio y el tiempo, se calcula usando la ecuación de Jeffery,

la que se escribe aquí en la forma Euleriana compacta. Mientras que la derivada convectiva (o material)

en el lado izquierdo (l.h.s.) de la ec. (1) describe el transporte de la orientación de las fibras, FO, puro en el campo de velocidad local U(r, t) del flujo, el lado derecho (r.h.s.) modela el movimiento de rotación de la partícula impulsada por el tensor

de gradiente de velocidad local efectiva . En el caso ideal especial de fibras infinitamente delgadas (es decir la relación de aspecto ra→∞) el tensor

es idéntico al tensor gradiente de velocidad

<∞) se define por el término ar, en el caso general de fibras que tienen una relaciónde aspecto finita (< U
ٔ ׏

que contiene el parámetro de geometría de la fibra

Este parámetro codifica cómo la geometría de la partícula afecta al movimiento de rotación de las fibras en el caso de partículas elipsoidales de simetría rotacional.

Vías alternativas para escribir la ecuación de Jeffery hacen uso de la descomposición definida única del tensor de gradiente de velocidad en sus partes simétricas y antisimétricas dadas por los tensores de tasa de cizallamiento y de vorticidad dados por

Los tensores de tasa de cizallamiento y de vorticidad se relacionan con el tensor de gradiente de velocidad efectiva mediante las identidades

Sustituyendo estas identidades en (1) y tomando en cuenta que pT ·

· p = 0 se obtiene la ecuación de Jeffery en

la forma alternativa

Aún otra forma de la ecuación de Jeffery puede obtenerse reescribiendo el primer término de la r.h.s de (1d) usando la identidad - · p =

T · p = ½ rot U × p, la que relaciona el tensor de vorticidad con el rotacional rot U del campo vectorial de velocidad de flujo (ver por ejemplo [3]).

En el caso de partículas cuya geometría es simétrica rotacionalmente pero no elipsoidal, la forma de la ec. (1) así como la definición del gradiente de velocidad efectiva se mantiene igual, pero la fórmula del parámetro de geometría

λ debe ser alterada, de manera que la relación de aspecto “geométrica” ra = ℓ/d como se definió anteriormente se

sustituye por una relación de aspecto efectiva

el valor apropiado de la que ha de determinarse experimentalmente. En este sentido, el parámetro de geometría λ puede considerarse como un parámetro del material. En el caso de partículas asféricas asimétricas la ec. (1) ha de sustituirse por un sistema acoplado de tres ecuaciones de tipo Jeffery que contienen tres parámetros de geometría [3].

2.2 Distribución macroscópica de la orientación de las fibras

En el caso de una simulación de moldeo por inyección para plásticos reforzados con fibras cortas realizada a nivel macroscópico, el número de fibras contenidas en una sola celda de cálculo (es decir un elemento de volumen del dominio de cálculo) es muy grande debido a las pequeñas dimensiones de las fibras individuales. Este punto se ilustra por la micrografía mostrada en la Fig. 4, la que representa una muestra típica tomada de una pieza hecha de plástico reforzado con fibras cortas. Por tanto el estado de FO local del material contenido en una celda de cálculo se describe por medio de un modelo macroscópico, es decir una función de distribución Ψ(p,r,t) del vector de FO p (ver [5]). Notamos que aunque la concentración de fibras se asume que es homogénea, la función de distribución Ψ aún depende paramétricamente de las coordenadas (r, t) mediante el campo de velocidad de flujo local U(r, t) y su

U. Comoun puntoadicional, un término que tiene encuenta la interacción mutua de las fibras, que
ٔ ׏ gradiente

influye en su orientación, tiene que incluirse en el modelo a nivel macro.

La transición de la ecuación de Jeffery como un modelo micro de fibras individuales, que no interactúan a un modelo macro el que da las estadísticas de FO de muchas fibras que interactúan se logra mediante la ecuación correspondiente de Fokker-Planck

con la función de distribución de la FO Ψ como variable dependiente. Aquí

sirve como un atajo para la r.h.s. de

esfera ·... simbolizan los operadores gradiente y divergencia definidos en lap
׏ ... yp
׏ la ecuación de Jeffrey (1), y

unitaria S2 del espacio euclidiano tridimensional

De acuerdo con el enfoque de Folgar y Tucker [4], la introducción de un coeficiente de difusión r = Clγeff proporcional a la tasa de cizallamiento efectivo

del gradiente de velocidad local da un modelo simple de interacción fibra-fibra en suspensiones concentradas. Aquí las Gij son las componentes del tensor de tasa de cizallamiento

como se define en (1b), y el argumento de la raíz cuadrada - usando la convención de suma de Einstein - es un atajo para su autocontracción

El coeficiente de interacción no negativo, adimensional CI es un parámetro del material de la suspensión. Típicamente tiene un valor pequeño (positivo) en el rango de 10-3 ... 10-2 para plásticos reforzados con fibras cortas. Notamos que en flujos incompresibles (así como casi incompresibles), la debilidad relativa del término de difusión “estocástica” -comparado con el término que representa la dinámica de Jeffery “determinista” -puede ser la fuente de problemas de estabilidad.

2.3 Tensores de orientación de las fibras y la ecuación de Folgar-Tucker

Calcular la distribución de la FO local por medio de la ecuación de Fokker-Planck (2) requiere la solución numérica de una PDE definida en la esfera unitaria S2 para cada celda de cálculo del dominio de simulación de flujo, lo que es

una tarea prohibitivamente cara para problemas tridimensionales de “tamaño industrial”. Por lo tanto, Advani y

Tucker [6] propusieron el uso de tensores de orientación de las fibras, los que se definen como momentos de la función de distribución, y por lo tanto para sustituir la ecuación de Fokker-Planck por una jerarquía de ecuaciones de momentos para los tensores de FO. Debido a la simetría de inversión (-p,...) = Ψ(p,...) de la distribución de la FO con respecto a la variable p, la que refleja el hecho de que las direcciones ±p corresponden al mismo estado de orientación, todos los momentos de orden impar se anulan idénticamente, de manera que la expansión de los momentos de Ψ contiene sólo elementos de orden par. El primer momento no trivial de esta expansión es por lo tanto el segundo dado por

(en notación de índice:

). La notación

significa que la integral se calcula sobre la superficie de la esfera S2, donde dS(p) es la medida de integración. El segundo momento â(2)(r, t) se llama tensor de FO de 2do orden (o matriz de FO) y es, por definición, una matriz 3×3 simétrica real. Como la distribución de la FO (también por definición) se normaliza de acuerdo con

â(2)

la matriz de FO tiene la propiedad obvia de que su traza as

El siguiente momento no trivial en la jerarquía de expansión es el tensor de FO de 4to orden â(4) definido como

información completa acerca de â(2).

La jerarquía de las ecuaciones de momentos mencionadas anteriormente se obtiene para los momentos de cada orden mediante el intercambio de las operaciones de derivación e integración, es decir

reemplazar por los términos en el r.h.s. de la ecuación de Fokker-Planck (2) y evaluar analíticamente las integrales correspondientes. Si este procedimiento se aplica a la matriz de FO â(2), la así llamada ecuación de Folgar-Tucker (FTE) se obtiene como la primera ecuación no trivial en la serie organizada jerárquicamente de ecuaciones de los momentos:

Como en la ecuación de Jeffery, la derivada convectiva en el l.h.s de (5) describe el transporte puro del estado de FO local (representado por la matriz de FO â(2) en este orden de la expansión de los momentos) debido al movimiento de traslación de las fibras en un sistema de referencia Euleriano, mientras que los primeros dos términos en el r.h.s de (5) modelan su movimiento de rotación impulsado por el gradiente de velocidad efectiva local . El 3er término del r.h.s. de (5) resulta de la presencia del término de difusión en la ecuación de Fokker-Planck (2). Al nivel de la FTE esto produce un efecto amortiguador que arrastra la matriz de FO hacia el estado isotrópico â(2) = ⅓ î.

En su forma Euleriana (5) la FTE es un sistema acoplado de PDE de 1er orden de tipo convección-reacción. Desde el punto de vista Lagrangeano (5) es un sistema acoplado de ODE para las componentes de â(2).

2.4 Propiedades básicas de la matriz de FO

Como una matriz 3×3 simétrica real la matriz de FO â(2) posee 3 valores propios μk reales con vectores propios

correspondientes Ek que forman una base ortonormal de

es decir: â(2) . Ek = µk Ek con Ek · El = δkl, donde δkl es el símbolo de Kronecker el que es igual a 1 para k = l y cero en caso contrario. Resulta directo mostrar que, por construcción, los valores propios μk de la matriz de FO â(2) son no negativos y – dado que se cumple que p2=1 - su suma, la que es idéntica a la traza de â(2) (ver más arriba), siempre es igual a uno (es decir:

Para cualquier vector unitario p0, la matriz de FO especial â0 = representa un así llamado estado de

0p
ٔ 0p

orientación uniaxial donde el 100% de las fibras se orientan en la dirección dada por ±p0. Por lo tanto, el signo de p0 se quita del producto diádico que define esta matriz de FO especial. El vector p0 es un vector propio de â0 con valor propio 1, y los restantes dos valores propios son cero con vectores propios correspondientes que yacen en el plano ortogonal a p0 y por lo demás son arbitrarios.

La matriz de FO se puede escribir en términos de su representación espectral

la que

formada por los vectores kE
ٔ kEtiene la forma de una suma ponderada de los estados de orientación uniaxiales

propios de la matriz de FO, con los pesos dados por los valores propios. Esto permite una interpretación de un valor

5 propio μk como la fracción local de fibras orientadas a lo largo de la dirección del correspondiente vector propio Ek. En este sentido los datos espectrales {µk, Ek} k=1.2..3 de la matriz de FO representan el estado de orientación macroscópica local de las fibras dentro de un pequeño volumen de la suspensión.

Esto sirve como una motivación para la siguiente definición (matemáticamente formal) de una matriz de FO y el espacio de fase de la FTE [12]:

Definición: Una matriz 3×3 simétrica real â es una matriz de FO si y sólo si es semidefinida positiva y su traza es igual a 1. El espacio de fase MFT de la FTE es el conjunto de todas las matrices de FO.

15 Se puede mostrar que el conjunto MFT es equivalente al conjunto de todas las matrices 3×3 simétricas, reales que resultan de las integrales de los momentos del tipo (3) usando funciones de distribución normalizadas apropiadamente, pero arbitrarias por lo demás. Una caracterización matemática del espacio de fase MFT por medio de sus propiedades topológicas y geométricas se ha dado recientemente por uno de los autores (ver [12]). Es de gran importancia práctica conocer con exactitud si una matriz pertenece al conjunto MFT, dado que la interpretabilidad de los resultados de una integración numérica de la FTE requiere que la variable dependiente â(2) tenga propiedades espectrales especiales como se describió anteriormente durante todos los pasos (es decir también en los intermedios) del cálculo de FO.

25 2.5 El problema de la clausura

El así llamado problema de la clausura se origina del hecho de que en cada orden de la expansión de los momentos, la DE para el momento â(2n) contiene el momento â(2n+2) del orden superior siguiente como una variable. Mientras que â(2n) puede expresarse en términos de â(2n+2) por medio de una simple identidad algebraica (es decir una suma sobre un par de dos índices iguales, ver más arriba), lo contrario no es posible, de manera que â(2n+2) ha de tratarse como una incógnita. En la FTE este problema de la clausura se manifiesta el mismo por la aparición de â(4) en el

r.h.s. de (5), lo que evita que el sistema sea soluble a menos que se cierre expresando â(4) como una función de â(2) por medio de una aproximación de clausura. Aplicar una aproximación de clausura a la FTE significa sustituir el tensor de FO de 4to orden exacto - pero desconocido - â(4) en el r.h.s. de (5) por alguna función tensorial apropiada

35 (en general no lineal) Â(4)[â(2)] de la matriz de FO [8].

Un ejemplo conocido es la clausura híbrida [7] definida por la fórmula

la que, a pesar de algunos inconvenientes bien conocidos, es una elección aceptada debido a su simplicidad algebraica y

su robustez numérica [7]. La clausura

se define como una interpolación (convexa) entre dos clausuras de 45 tipo más simple: la clausura cuadrática definida como

la que da resultados exactos en el caso especial de una distribución de orientación uniaxial, y la clausura lineal dada por

la que es exacta en el caso isotrópico. El peso de interpolación entre estos casos extremos se proporciona por el factor de orientación escalar

Como el determinante det(â(2)) es un invariante de la matriz de FO, lo mismo es cierto para el factor de orientación escalar.

Denotamos la variante especial de la FTE (5) combinada con la aproximación de clausura híbrida (6) por el acrónimo FTE-hyb. Esta variante es de gran interés práctico, la modalidad preferida de la invención también se basa en la variante FTE-hyb del modelo de Folgar-Tucker.

2.6 La FTE como un sistema algebraico diferencial

El r.h.s. de (5) es formalmente bien definido para matrices simétricas reales arbitrarias. La cuestión de si las propiedades del r.h.s. de (5) confinan automáticamente las trayectorias de la solución al dominio de MFT puede responderse positivamente al menos para la FTE en su forma exacta (es decir sin una aproximación de clausura): Si el tensor de FO de 4to orden â(4) exacto se inserta en el r.h.s. de (5) y entonces â(2) se calcula como una solución de (5), el resultado es idéntico al 2do momento exacto de la distribución de la FO el que se obtiene mediante la evaluación de la integral de los momentos (3) usando la solución Ψ de la ecuación de Fokker -Planck (2) con los mismos parámetros

y Dr. Debido a este argumento puede concluirse que la solución de (5), usando el tensor de FO de 4to orden â(4) exacto, puede escribirse como una integral de los momentos (3), y su trayectoria por lo tanto está necesariamente confinada al dominio de MFT MFT.

Sin embargo, este argumento ya no es válido si se aplica una aproximación de clausura, lo que siempre es necesario si se tiene que resolver numéricamente la FTE sin ningún conocimiento previo de la función de distribución completa. De esta manera, el problema del confinamiento de las soluciones de la FTE al espacio de fase MFT siempre surge en todos los problemas de interés práctico.

La necesidad de confinar las trayectorias de solución de la FTE al dominio de MFT impone restricciones adicionales sobre la variable dependiente â(2) las que pueden formularse en términos de desigualdades algebraicas de sus invariantes (ver [12]). Estas restricciones sobre la variable dependiente convierten la FTE en un sistema algebraico diferencial (DAS) y ha de tenerse cuidado en el procedimiento usado para la integración numérica.

3. Caracterización matemática del espacio de fase de la FTE

Las propiedades topológicas básicas del espacio de fase MFT puede deducirse directamente de las propiedades especiales de las matrices de FO como se resume en

Teorema 1: El espacio de fase MFT es un subconjunto convexo acotado del espacio vectorial de las matrices 3×3 simétricas reales confinado al hiperplano de cinco dimensiones definido por la condición de traza Tr(â)=1.

La convexidad de MFT permite la definición de una proyección de mapeo sobre MFT la que puede combinarse con cualquier integrador de ODE adecuado para construir un método de integración apropiado para la FTE (ver el capítulo IV.4 de [9]). Una caracterización algebraica invariante de las matrices de FO puede obtenerse por un análisis del polinomio característico

de una matriz 3×3 simétrica real.

Los coeficientes

y Da = det(â) son invariantes de la matriz. (En la literatura estos invariantes se denotan ocasionalmente por I1 = Sa ' I2 = Ka y I3 = Da). La caracterización algebraica de las matrices de FO puede formularse en términos de estos invariantes de acuerdo con

Teorema 2: Una matriz 3×3 simétrica real â es una matriz de FO si y sólo si su traza Sa es igual a 1 y sus invariantes Ka y Da son no negativos.

La Fig. 6 ilustra esta afirmación mediante un análisis separado de los términos cuadrático y cúbico y los términos lineal y constante del polinomio característico Pa(μ): una traza positiva Sa obviamente proporciona la existencia de valores propios positivos μk, mientras que los valores no negativos de los dos invariantes Ka y Da evitan la existencia de negativos. Aún más, las condiciones Sa > 0, Ka ≥ 0 y Da ≥ 0 son no sólo necesarias sino también suficientes para que todos los valores propios de â sean no negativos, dado que la matriz â siempre tiene tres valores propios reales. Combinado con la condición de traza Sa=1 esto completa la prueba del teorema anterior.

Una noción de la geometría del espacio de fase MFT, el que es un objeto de 5 dimensiones de acuerdo con el teorema 1, puede obtenerse observando por separado a los elementos de la diagonal y fuera de la diagonal de las

matrices de FO. Primero notamos que los elementos de la diagonal dados por

son siempre no negativos y satisfacen la condición de traza Σkakk=1. Por lo tanto, independiente de la elección del sistema de coordenadas, ellos se confinan al conjunto triángulo (ver Fig. 7)

Esto proporciona una condición necesaria la que caracteriza las matrices de FO con respecto a sus elementos de la diagonal y da una representación invariante de la “parte diagonal” del conjunto de espacio de fase MFT.

Una descripción formal de la parte “fuera de la diagonal” del dominio de MFT puede obtenerse mediante la introducción de la notación (x, y, z) y(u, v, w) para las ternas de los elementos de la diagonal y fuera de la diagonal

ay que formalmente Δ
א )x, y, zde una matriz 3×3 simétrica real, tomando una terna diagonal arbitraria (pero fija) (

define el conjunto

de todas las ternas fuera de la diagonal “admisibles” que pertenecen a una terna fija de la diagonal. Algebraicamente

el conjunto N(x,y,z) puede caracterizarse como el conjunto de todas las ternas fuera de la diagonal (u, v, w) que satisfacen simultáneamente el siguiente par de desigualdades, como se explicó anteriormente con respecto al Teorema 2:

Un ejemplo especial correspondiente al caso “máximamente simétrico” x = y = z=1/3 se muestra en la Fig. 8, otros casos (ver la Fig. 9) corresponden a versiones distorsionadas de este volumen en forma de tetraedro.

Para complementar las restricciones (9) y (10) de los invariantes Ka y Da a valores no negativos, notamos que para cualquier matriz 3×3 simétrica real con traza positiva Sa >0 estos invariantes se restringen también a partir de lo anterior

por

de manera que para las matrices de FO (Sa=1) sus valores siempre se confinan por 0 ≤ Ka ≤ ⅓ y D y a ≤ 1/27. La fibración (formal) de MFT = y N(x,y,z)

aobtenida por la introducción de los conjuntos Δ )z,y,x(N) ×z,y,x(Δa
א )z,y,x(U

ayuda a obtener una imagen de la estructura general de MFT (ver la Fig. 10): Las “fibras” individuales pueden visualizarse mediante un procedimiento que adjunta localmente el conjunto N(x,y,z) de ternas admisibles fuera de la

, y la variación subsecuente de este punto base a través de todo el conjuntoaΔ
א diagonal a su “punto base” (x, y, z)

de triángulo completo Δa cubre el espacio de fase completo.

Los resultados resumidos en esta sección muestran claramente que el espacio de fase MFT es un objeto matemático complejo. De acuerdo con el resultado matemático riguroso establecido en el Teorema 2, el confinamiento de cualquier trayectoria de solución t հ â(2) (t) de la FTE, la que naturalmente se define como una ecuación de evolución en el espacio vectorial real de las matrices 3×3 simétricas, al dominio del espacio de fase necesariamente implica el cumplimiento de las desigualdades invariantes (9) y (10) combinadas con la condición de traza unitaria. Este hecho inevitablemente convierte la FTE en un DAS. A diferencia de la condición de traza unitaria, cuyo cumplimiento está “integrado en” la FTE exacta en la ausencia de cualesquiera aproximaciones de clausura para el tensor de FO de 4to orden y la que es aún válida sobre prerequisitos más generales para una clase grande de aproximaciones de clausura (ver la sección 4), la validez de las desigualdades invariantes (9) y (10) generalmente no se conserva en la presencia de cualesquiera de las aproximaciones de clausura conocidas actualmente.

(También, la validez de (9) y de (10) en el caso de la FTE “exacta” puede deducirse más bien por un razonamiento

indirecto como se da en el sección 2.6 que por medio de la estructura algebraica de la propia FTE). Estos hechos matemáticos comúnmente se pasan por alto en los ejemplos simples usualmente mostrados en el trabajo académico.

Sin embargo, una aproximación de clausura se contiene necesariamente como un elemento esencial de cualquier procedimiento que resuelve la FTE numéricamente. Por lo tanto cualquier procedimiento de simulación adecuado para manipular situaciones más complejas similares a las que ocurren típicamente en las aplicaciones industriales debe tener esto en cuenta por medio de procedimientos de control adecuados que confinen la trayectoria de solución a su dominio admisible teóricamente. Por lo tanto un tratamiento apropiado de la FTE como un sistema algebraico diferencial es obligatorio para las aplicaciones industriales serias.

En la sección 7.7 se presenta un procedimiento que proporciona un tipo adecuado de control invariante. Este procedimiento se ha implementado en una modalidad preferida.

4. El problema de la conservación de la traza y la estabilidad

La aplicación de una aproximación de clausura a la FTE significa sustituir el tensor de FO de 4to orden â(4) exacto (pero desconocido) en el r.h.s de (5) por alguna (en general no lineal) función tensorial Â(4)[â(2)] de la matriz de FO (es decir: â(4) ← Â(4)[â(2)] en notación matemática, ver la sección anterior sobre el problema de la clausura para detalles adicionales). En dependencia de la elección de la clausura, el tensor de 4to orden Â(4) posee una cantidad menor o mayor de las propiedades del tensor exacto â(4).

Un requisito especial concerniente a las propiedades de simetría de Â(4) el que se cumple por todas las aproximaciones de clausuras razonables es la así llamada simetría ortotrópica como se da por el conjunto de identidades

es decir se requiere que el tensor de 4to orden Â(4) sea simétrico con respecto al primer y segundo pares de índice (ij) ↔ (kl) así como al intercambio (ij) ↔ (kl) de ambos pares. Si se asume que â(2) es una solución de (5) aumentada

â(4)

por una aproximación de clausura ← Â(4)[â(2)] con propiedades de simetría ortotrópica (11), podemos formalmente derivar la siguiente ecuación diferencial (abreviado: DE) para su traza (ver [17]):

La condición de traza Tr(â(2)) = 1 no se usó aún en la derivación de (12), dado que la estabilidad y la conservación de la traza para las soluciones de la FTE se investigarán mediante un análisis de esta ecuación. Si la aproximación de clausura se aplica a la FTE, la última adicionalmente satisface las condiciones de normalización

y la DE (12) para la traza se simplifica a

Como el parámetro de difusión Dr es por definición no negativo, puede concluirse que o bien (en el caso de Dr > 0) el hiperplano definido por la condición de traza es una variedad integral estable o que (en el caso de Dr = 0) la traza es una primera integral de la FTE, siempre que la aproximación de clausura â(4) ← Â(4)[â(2)] satisfaga ambas condiciones (11) y (13). (Debería notarse que el tensor de FO de 4to orden exacto â(4) es por definición totalmente

simétrico y satisface ∑ y que la matriz de FO exacta â(2) es por definición simétrica y satisface la condición de traza)

Las consideraciones anteriores muestran cómo la aplicación de una aproximación de clausura a la FTE puede influir en la validez de la condición de traza, la que es la más simple de las condiciones impuestas a los invariantes de una matriz de FO, en dependencia de la cantidad de propiedades que la función tensorial de aproximación Â(4)[â(2)] y el tensor exacto â(4) tienen en común.

La estabilidad de la traza tiene que considerarse cuidadosamente en el caso especial de la clausura híbrida (6), la que posee sólo las propiedades de simetría (11) pero no satisface las condiciones de normalización (13). En el caso de esta aproximación de clausura, la DE para la traza asume la forma (ver también [17])

con un prefactor

Este prefactor se desvía de la expresión más simple 6Dr válida para todas las clausuras que satisfacen (13), lo que refleja el hecho de que la clausura híbrida no lo hace.

De aquí en adelante se muestra que bajo ciertas condiciones los términos adicionales que aparecen en el prefactor φa pueden provocar que este término se haga negativo. Esto significa que aunque el hiperplano Tr(â(2) = 1 es aún una variedad integral de la FTE con clausura híbrida, puede convertirse en inestable localmente.

La posibilidad de que el prefactor øa se haga negativo puede deducirse por el siguiente procedimiento de cuatro etapas de argumentos:

(i) En el caso especial de un campo de flujo incompresible (es decir la Tr( ) = div U = 0) el prefactor asume la

: â(2)

forma simplificada φa =6Dr + fs 2λ(

.

(ii) En el caso de fibras delgadas siempre puede asumirse que λ=1, y si el estado de orientación local descrito por â(2) es aproximadamente uniaxial (es decir uno de sus valores propios μk es cercano a 1), el valor del factor de orientación escalar fs =1-27det(â(2) también es aproximadamente 1 dado que el determinante de â(2) es cercano a cero. De esta manera se obtiene la expresión simplificada øa ≈6Dr +

2(

: â(2)), la que aproxima el valor exacto de øa bajo las suposiciones especificadas.

kE
ٔ kEkµk= Σ(2)â y la inserción de la representación espectral effγl = C rDMediante la sustitución de (iii)

en la contracción se obtiene una forma modificada de la expresión aproximada para el prefactor:

(iv)

En el último paso se normalizan las componentes del tensor de la tasa de cizallamiento a la tasa efectiva de cizallamiento

El tensor de la tasa de cizallamiento normalizado tiene los mismos vectores propios que , traza cero y por lo tanto – justo los mismos de -valores propios negativos. Dado que las componentes de son del orden de 1 por construcción, el tensor necesariamente tiene al menos un valor propio negativo con un valor absoluto del orden de 1.

Este hecho es la clave para el estimado deseado de φa bajo las varias circunstancias asumidas hasta ahora en la derivación de la expresión aproximada para el prefactor φ a. Usando el tensor de la tasa de cizallamiento normalizado

esta expresión aproximada puede escribirse en la forma

Tomando en cuenta que en las simulaciones prácticas el coeficiente de interacción Cl es un número positivo pequeño (típicamente 0< Cl < 10-2, ver también la sec. 1.2), puede verse que φa se hará negativo si se evalúa con una matriz de FO la que tiene un valor propio dominante μj ≈ 1 con vector propio correspondiente Ek que apunta a lo largo de la dirección propia del mayor valor propio negativo del tensor de tasa de cizallamiento, como en este caso puede ser

estimado

y por lo tanto se obtiene φa ≈6γeff [Cl - 2] < 0.

Los argumentos anteriores, pueden formularse de manera rigurosa. De esta manera es posible dar una demostración matemática de que para cualquier flujo incompresible y valores 0 Cl < λ del coeficiente de interacción la traza de una solución de la FTE con clausura híbrida se hace localmente inestable en ciertas regiones del espacio de fase MFT, dado que la expresión (15b) para el prefactor øa se hace inevitablemente negativa en estas regiones. Por lo tanto, es claramente necesario tener cuidado con el problema de la estabilidad de la traza durante la integración numérica de la FTE.

5. La integración numérica de la FTE: aspectos generales

La elección de un método adecuado para la integración numérica de la FTE (incluyendo cualquier tipo de aproximación de clausura) depende de la clasificación matemática del tipo de ecuación así como los aspectos relacionados con la forma algebraica específica de la FTE. La estructura general de la FTE cerrada muestra que constituye un sistema acoplado de ecuaciones diferenciales parciales (PDE) hiperbólicas de tipo convección

reacción para los elementos

de la matriz de FO como variables dependientes, como se explica en detalle en el párrafo siguiente.

Comenzando a partir de la forma de la FTE dada en la ec. (5), esta ecuación puede reescribirse sustituyendo

explícitamente la derivada material

en el l.h.s. y la notación matemática â(4) ← Â(4)[â(2)] que simboliza una aproximación de clausura en el r.h.s. de (5), lo que da una formulación equivalente de la FTE cerrada:

Esta forma de la FTE ya revela la mayoría de su estructura matemática: el l.h.s. consiste de un simple operador diferencial parcial de primer orden de tipo de transporte o de convección con la velocidad de flujo local U(r, t) que

gobierna el transporte desacoplado, puramente convectivo de los elementos

de la matriz de FO que son las variables dependientes de la ec. (16). Aún más la estructura algebraica de los varios términos del r.h.s. de (16) muestra que el r.h.s. – justo como la propia matriz de FO, y como se requiere por la consistencia matemática constituye una función tensorial simétrica de 2do orden, la que se denota como (2)(...) en la siguiente forma resumida (pero equivalente) de la ec. (16):

Aunque la dependencia explícita del r.h.s. de la matriz de FO â(2) es lineal, la naturaleza de su dependencia implícita, la que se determina por la forma funcional â(4) ← Â(4)[â(2)] de la aproximación de clausura, en general es no lineal. Por lo tanto la función

(2)(...)

ha de considerarse como una función no lineal de la variable dependiente â(2) a menos que se elija una aproximación de clausura lineal. Además puede reconocerse que la función

(2)(...)

en el r.h.s. conduce inevitablemente a un acoplamiento mutuo de las ecuaciones individuales

para los elementos de FO

a menos que el gradiente de velocidad efectiva

sea una matriz diagonal y la función de clausura Â(4)[â(2)] tenga una estructura algebraica muy especial.

y el parámetro de U
ٔ׏ del tensor de gradiente de velocidad existente

Después de resolver la dependencia degeometría de la fibra λ por la ec. (1a) y resustituyendo la expresión que define Dr = Clγeff para el parámetro de difusión de rotación en términos de la tasa de cizallamiento efectiva, la ec. (17a) puede reescribirse en la forma

Esta forma de la FTE cerrada revela explícitamente la dependencia del r.h.s. de los parámetros constantes del modelo λ y Cl, pero aún diferencia entre la dependencia explícita e implícita de la matriz de FO así como entre la

inducida por la U
ٔ y la dependencia implícita de V U
ٔ dependencia explícita del tensor de gradiente de velocidad V

fórmula de la tasa de cizallamiento efectivo

(ver la sección de la ecuación de Fokker-Planck para más detalles). En la literatura la ec. (17b) a veces se escribe en la forma de componentes usando notación de índice, es decir

Frecuentemente se encuentra una versión muy simplificada

de esta ecuación, usando una notación condensada que oculta todas las dependencias de los parámetros y que no discrimina entre las dependencias explícita e

implícita de las funciones componentes (...) de la matriz de FO o del gradiente de velocidad. Revirtiendo de la notación de índice a la notación tensorial completa usada en las ec. (17a/b), se obtiene la versión simplificada

de estas ecuaciones.

6. La integración numérica de la FTE: “División de operador”

de enfoques para la integración numérica de un sistema de PDE acoplado

de tipo convección-reacción. Un enfoque que ha probado ser robusto y flexible especialmente en el caso de funciones en el r.h.s. no lineales F(...) comienza a partir de una reformulación equivalente de la PDE

y procede mediante un tratamiento separado de los dos términos del r.h.s. dados por

(...). En la literatura matemática este F... y el “operador” determinado por la función
׏ ·Uel operador diferencial lineal

enfoque se conoce como [23-37] “método de pasos fraccionados”, “método de división” o simplemente “división de operador” (ver [28-31] para enfoques alternativos). La integración numérica de una ecuación que tiene la forma

mediante la “división de operador” hace uso de las soluciones numéricas (o en algunos

casos incluso analíticas) de las dos ecuaciones separadas

y

obtenidas de la ecuación completa estableciendo cualquiera de los dos operadores en el r.h.s. a cero. La primera ecuación

es equivalente a la ecuación de convección homogénea,

la segunda constituye un sistema de ecuaciones diferenciales ordinarias (ODE) el que es en general acoplado y no lineal. Este tipo de enfoque se usa en la modalidad preferida. Usando la notación de las ecs. (17a/b), la aplicación del enfoque de división a la FTE cerrada conduce a las ecuaciones parciales

( “paso de convección” con una matriz nula en el r.h.s)

las que se tratan como subproblemas separados dentro del marco del método “división de operador”. El efecto físico modelado por la ec. (18a) es la redistribución espacial de las estadísticas de FO (como se codificada por los elementos de la matriz de FO) dentro del dominio lleno de la cavidad del molde debido a un transporte puramente convectivo de la masa fluida de acuerdo con la velocidad de flujo, mientras que la ec. (18b) desprecia completamente los efectos del transporte convectivo y solamente considera la tasa local de cambio de la distribución de la FO debido a la cinemática de rotación de las fibras impulsadas por el gradiente de velocidad local así como la interacción mutua entre las fibras. De aquí en lo adelante se describen varias variantes que muestran cómo las soluciones (numéricas) de ambos subproblemas pueden combinarse para dar una solución aproximada de la ecuación completa (16) (o su equivalente (17a/b) en notación abstracta).

6.1 “División simple de operador”

Usando la manera más sencilla (usualmente denotada como “división simple de operador”) para obtener una solución aproximada de la ecuación a ser resuelta en realidad a partir de las dos ecuaciones parciales, la última se resuelve en orden consecutivo, tomando la solución intermedia resultante de la primera como entrada (es decir valor inicial) para la segunda. Para explicar esto en detalle, simplificamos la notación usando la ecuación modelo

con las identificaciones y ↔ â(2) and F(y) ↔

(2)(â(2) ,...).

El solucionador de flujo, es decir el software que modela el flujo de la suspensión, da las variables de estado del

), medianteuna solución de last,r(U
ٔ ׏ ) y su tensor gradiente t,r(Ufluido, entre ellas la velocidad de flujo

ecuaciones de transporte discretizadas para la masa, el momento y la energía, en instantes discretos de tiempo tn en todas las celdas de cálculo localizadas alrededor de los puntos discretos rm en el espacio los que se contienen en la parte llena Ω(n) del dominio de cálculo total que cubre la cavidad del molde, el sistema de respiraderos y la entrada. El solucionador de flujo comienza en el instante t0 = 0 y procede desde el instante tn hasta el instante tn+1, usando un tamaño de paso tn+1 := tn+1 -tn. Los cálculos de flujo implican un cálculo del dominio Ω(n+1) en el instante tn+1 que depende de su estado anterior Ω(n) y el nuevo campo de velocidad U(n+1). Esto puede hacerse usando un método de volumen de fluido (VOF). Calcular el nuevo frente de flujo y el nuevo dominio Ω(n+1) da nuevos valores de todas las variables de estado en las celdas de Ω (n+1) en el instante tn+1. En un paso correspondiente, los nuevos valores

en el instante tn+1 en todas las celdas de cálculo localizadas alrededor de los puntos rm del dominio recién lleno Ω(n+1) deben calcularse mediante una solución numérica de la PDE

comenzando a partir de los valores

definidos en el dominio “viejo” Ω(n) los que ya se conocen desde el paso de cálculo anterior

Esto se logra mediante el siguiente procedimiento de tres pasos, el que constituye una posible variante de “división simple de operador”:

1. Paso de extensión: Comenzando desde

en las celdas del dominio Ω(n), calcula una

inicialización en todas las celdas del nuevo dominio Ω(n+1).

2. Paso de convección: Comenzando desde los valores iniciales

proporcionados por el paso de

extensión, calcular los valores intermedios

mediante una solución numérica de la ecuación de

convección

con un tamaño de paso Δtn+1, usando la velocidad de flujo U(rm,tn+1) calculada por el solucionador de flujo en el dominio Ω(n+1).

3. Paso de rotación: Comenzando desde los valores iniciales

proporcionados por el paso de

convección anterior, calcular los valores finales

mediante una solución numérica del

)+1nt,mr(U
ٔ ׏ , usando el gradiente de velocidad +1ntcon tamaño de paso Δ

sistema de ODE proporcionado por el solucionador de flujo en el dominio Ω(n+1).

La Fig. 11 ilustra esta variante de división de operador.

La variante de “división simple” co mo se describió anteriormente da una solución aproximada de la ecuación

completa

con un error de discretización de tiempo de 1er orden -es decir O(Δtn+1) -en el tamaño del paso y se usa en el método de acuerdo con la modalidad preferida. Los detalles sobre los tres pasos de este procedimiento se explicarán en las secciones subsecuentes.

Una variante alternativa de “división simple de operador”

,nt en el instante U
ٔ ׏ usando el gradiente de velocidad “viejo” )n(comienza con un paso derotación en Δ (i)

entonces

(ii) procede con un paso “intermedio” para extender los resultados de (i) desde Ω(n) hasta Ω(n+1) y

(iii) termina con un paso de convección usando la velocidad de flujo U en el instante tn+1.

En esta variante, los pasos de convección y rotación se ejecutan en orden inverso. Aunque el error teórico de discretización de tiempo es del mismo orden que en el caso de la primera variante descrita anteriormente, esta variante alternativa tiende a ser menos consistente, dado que se usan la velocidad de flujo y el gradiente de velocidad de diferentes instantes de tiempo. Esto puede evitarse mediante otra variante

(i)

comenzando con un paso de extensión de los valores conocidos

desde Ω(n) hasta Ω(n+1) como se describió anteriormente, entonces

(ii)

realizar un paso de rotación seguido por

(iv)

un paso de convección

+1)n(en Ωn+1tU proporcionados por el solucionador de flujo en el instante
ٔ ׏ usando los valores “nuevos” de U y

para ambos pasos(ii) y (iii).

6.2 “División simétrica de operador”

Teóricamente puede obtenerse un orden de precisión mayor con respecto al error de discretización de tiempo por medio de un enfoque de “división simétrica de operador”. La idea básica de este enfoque es intercalar un paso intermedio de una de las ecuaciones parciales con el tamaño de paso completo entre los pasos de las otras ecuaciones con la mitad del tamaño del paso. Una variante posible de este enfoque resulta en el siguiente procedimiento* de cuatro pasos:

1. Paso de extensión: Comenzando a partir de

en las celdas del dominio Ω(n), calcular una

inicialización

en todas las celdas del nuevo dominio Ω(n+1).

2. Paso de prerotación: Comenzando a partir de la extensión

de los valores conocidos para Ω(n+1) proporcionados por el paso de extensión, calcular los valores pre rotados mediante una solución

numérica del sistema de ODE

con la mitad del tamaño del paso Δtn+1/2, usando el gradiente

.+1)U(rm,tn+1) proporcionado por el solucionador de flujo en el dominio Ω (n
ٔ ׏ de velocidad

3. Paso de convección: Comenzando a partir de

como valores iniciales proporcionados por el paso de

prerotación, calcular los valores intermedios

mediante una solución numérica de la ecuación de convección

con tamaño de paso completo Δtn+1, usando la velocidad de flujo U(rm, tn+1) calculada por el solucionador de flujo en el dominio Ω(n+1).

4. Paso de postrotación: Comenzando a partir de los valores iniciales

proporcionados por el paso de

convección precedente, calcular los valores finales mediante una solución numérica del sistema de

)+1nt,mr(U
ٔ /2, usando el gradiente de velocidad +1ntcon la mitad del tamaño del paso Δ

ODE

proporcionado por el solucionador de flujo en el dominio Ω(n+1).

La variante de “división simétrica” como se describió anteriormente da una solución aproximada de la ecuación

completa

la que tiene un error de discretización de tiempo de 2do orden -es decir

- en el tamaño del paso.

Otra variante de “división simétrica” que es teóricamente posible intercala un paso de rotación entre dos pasos de convección con la mitad del tamaño del paso. Esta comienza

(i)

con una extensión de los valores conocidos

al nuevo dominio Ω(n+1),

(ii)

procede con un paso de “pre convección” de la mitad del tamaño del paso en el nuevo dominio usando la velocidad de flujo U(rm,tn+1), entonces

) y+1nt,mr(U
ٔ ׏ usandoel gradiente de velocidad +1)n(realiza un paso de rotación completo sobre Ω (iii)

(v)

termina con otro paso “post convección” de la mitad del tamaño del paso usando la

-

velocidad de flujo U(rm, tn+1),

la que finalmente resulta en una aproximación de

que también tiene una precisión de 2do orden. Como en el caso de “división simple” aún hay una tercera variante que intercambia el paso de extensión inicial y el

paso de prerotación, comenzando con el último y de esa manera usando los valores “viejos” del gradiente de

.(n)) en ΩU(rm, tn
ٔ ׏ velocidad

Como las situaciones de flujo típicas encontradas en el llenado de moldes frecuentemente involucran un flujo de

cizallamiento el que conduce a una rotación significativa de las fibras, la parte “operador de rotación” siempre juega

un rol significativo - y dominante en la mayoría de las situaciones -comparado con los efectos de un transporte puramente convectivo. Por lo tanto el intercalamiento simétrico del paso de convección entre los pasos de prerotación y postrotación como se describió en detalle para la primera de las variantes de división simétrica

constituye la variante preferida del enfoque de “división simétrica” [19]. Sin embargo, en la mayoría de los casos

prácticos, el tamaño del paso de tiempo proporcionado por los cálculos de flujo de fluido es suficientemente

pequeño, y la variante de “divisón simple” como se presentó en la sección 6.1 trabaja con suficiente precisión en

estos casos.

6.3 Variantes alternativas de división

En dependencia de la estructura algebraica de la función F(...) en el r.h.s., son posibles una variedad de enfoques alternativos de división. Si la función consiste en una suma de términos F(...) = ΣkFk(...), es posible dividir el r.h.s. en

un nuevo conjunto de subecuaciones

determinadas por las funciones parciales Fk...). Posiblemente algunas de las funciones parciales involucren operadores lineales, es decir Fk(y) =

k·y. En este caso se puede considerar algunos (o todos) estos términos lineales en el r.h.s de la ecuación de convección parcial, es decir

la que de esta manera se mantiene lineal, y dejar algunas (o todas) de las funciones

genuinamente no lineales Fj(...) para el tratamiento dentro de una sola ODE o una lista:

de ODE separadas como se explicó anteriormente.

Observando el r.h.s. de la ec. (16) se reconoce una suma de tres términos, siendo el primero y tercero lineales, y el del medio que involucra una aproximación de clausura que posiblemente tiene una dependencia no lineal de la matriz de FO. (Dependiendo de la elección de la clausura este término del medio se pudiera dividir en otra suma de subtérminos). De acuerdo con los comentarios anteriores, la estructura de términos en la función del r.h.s. de (16) abre la posibilidad de construir una variedad de esquemas adicionales de división. En contraste con las variantes basadas en la división de las ecs. (18a/b), a estas alternativas adicionales no le corresponden una distinción clara de

efectos físicos mutuamente separables (como en el caso de “transporte convectivo” frente a “dinámica de la rotación” local). Por lo tanto pueden considerarse meramente como “divisiones matemáticas artificiales”, las que son

difícilmente útiles para una simulación apropiada del fenómeno físico de la FO. La posibilidad de estas variantes de división “artificial” se menciona aquí solo para completitud y no se discutirán más aquí.

6.4 Enfoque numérico para el “paso de convección”

Esta sección explica el método numérico usado para integrar la ecuación de convección (18a), la que a de

resolverse como un subproblema dentro del enfoque de “división de operador” para la FTE cerrada. Basado en la

velocidad U(n+1) calculada para el nuevo paso de tiempo se ensambla la matriz del sistema de ecuaciones de transporte usando un esquema al viento de primer orden [32], y el sistema de ecuaciones resultante se resuelve para cada componente de la matriz de FO en el dominio de Ω(n+1).

6.5 Enfoque numérico para el “paso de extensión”

Se usa un método numérico para extender los valores conocidos de la matriz de FO dentro de las celdas de cálculo del dominio “viejo” Ω(n) en las celdas del dominio “nuevo” Ω(n+1). El cálculo del dominio “nuevo” Ω(n+1) se realiza por un método VOF modificado como una parte del algoritmo aplicado en el método de acuerdo con la modalidad preferida para la solución numérica de las ecuaciones de flujo (es decir las ecuaciones de conservación de la masa, del momento y de la energía). Los valores de la matriz de FO se quedan sin cambiar durante el paso de extensión en

aquellas celdas que son comunes a dos dominios. En general el dominio Ω(n+1) contiene un número de “celdas recién llenas” donde se necesita una asignación inicial de valores razonables para los elementos de la matriz de FO. De

acuerdo con la modalidad preferida esto se logra por medio de un promedio ponderado correspondiente al transporte de masa dentro de una celda de cálculo desde/hacia sus celdas vecinas de acuerdo con el flujo de masa calculado previamente.

Como la concentración de partículas se asume que es homogénea y la matriz de FO resulta de un procedimiento de promedio de volumen a un nivel macroscópico, se asume que una celda obtendrá una contribución de FO de sus celdas vecinas proporcional a la cantidad neta de masa de fluido transportada hacia dentro de ella. De acuerdo con la caracterización matemática resumida previamente del espacio de fase, es decir el conjunto de todas las matrices de FO posibles, de la FTE, lo último es un subconjunto acotado y convexo del espacio vectorial de seis dimensiones de las matrices 3×3 simétricas, reales confinado al hiperplano (de cinco dimensiones)

. Como un promedio ponderado de acuerdo con el transporte de masas se constituye en una combinación convexa de matrices de FO y por lo tanto siempre resulta en una matriz de FO válida, el procedimiento de inicialización es compatible con la estructura topológica (teóricamente requerida) del espacio de fase del modelo de FT, la que es una propiedad importante de este procedimiento.

6.6 Tratamiento del efecto “Flujo de fuente”

El término “flujo de fuente” caracteriza el patrón de flujo macroscópico total cerca de la superficie libre en el frente en

el caso de una gran clase de fluidos viscosos, no newtonianos. En un “flujo de fuente” las partículas aguas arriba

cerca del frente de flujo se transfieren desde la región central a las fronteras de las paredes. Este efecto realmente

sucede “automáticamente” debido a las propiedades materiales del fluido y no requiere ningún modelado adicional,

es decir: el “flujo de fuente” es un fenómeno emergente, al menos en los cálculos de flujo tridimensional basados en las ecuaciones de Navier-Stokes con leyes constitutivas no newtonianas apropiadas. Sin embargo, este no es el caso si los cálculos de flujo se basan en modelos simplificados (es decir del tipo Hele-Shaw como se discute en la sección 11.3 del artículo de Crochet en [1]). En simulaciones de un canal de flujo con el método de acuerdo con la

modalidad preferida se puede reconocer claramente un patrón de “flujo de fuente” en las trazas de las partículas, lo

que muestra que este fenómeno se tiene en cuenta por el solucionador de flujo.

Se necesita hacer algunas acciones especiales con el objetivo de asegurar que la componente normal de FO a la

superficie libre sea cero (como se requiere para la consistencia de los efectos observados de “flujo de fuente” en la

orientación de las fibras) para el procedimiento de inicialización de la orientación de las fibras para las celdas recién llenas. La anulación “requerida" de la componente normal de FO no se refuerza (intencionalmente) mediante el procedimiento “paso de extensión” descrito anteriormente por las siguientes razones: desde el punto de vista matemático, el modelo de FTE es un sistema de ecuaciones de transporte hiperbólicas acopladas mediante sus

r.h.s. (es decir un sistema de tipo convección-reacción con el término de reacción no lineal, ver más arriba). Por lo tanto no es apropiado - es decir matemáticamente hablando: imposible - prescribir la anulación de la componente de la matriz de FO en la dirección de la normal a la superficie libre como una condición de frontera para la matriz de FO en las celdas de superficie libre. La anulación de esta componente de FO debería de hecho emerger automáticamente a partir del comportamiento calculado de “flujo de fuente” del flujo.

Sin embargo, debido a las imprecisiones numéricas introducidas por los errores de discretización temporal y/o espacial así como al procedimiento usado para calcular el movimiento de la superficie libre durante el llenado del

molde, posiblemente se tienen que adicionar algunas “correcciones” para asegurar que el procedimiento de

inicialización para las celdas recién llenas sea consistente con el comportamiento de la FO en la superficie libre como se requiere por el fenómeno de “flujo de fuente”. En todas las celdas recién llenas se chequea que la traza del tensor de orientación de fibras difiera en no más de 1% del valor teórico de 1. Si la diferencia es demasiado grande el tensor de orientación de las fibras se corrige mediante una proyección ortogonal al triángulo admisible de los vectores propios en su espacio propio.

6.7 Inicialización de la matriz de FO al inicio de la simulación

La estructura hiperbólica de la FTE (cerrada) requiere la inicialización de la matriz de FO en las celdas de cálculo en (así como cerca de) la entrada en cada paso de tiempo de la simulación del llenado del molde como una condición de frontera. Por lo tanto tiene que hacerse una elección apropiada de este estado inicial de FO. En el método de acuerdo con la modalidad preferida se elige un estado de FO isotrópico â(2)=⅓

(correspondiente a una distribución aleatoria de FO) para este propósito como se explica más adelante.

En un flujo de cizallamiento altamente viscoso, el estado de FO se afecta fuertemente por una alta tasa de cizallamiento, cuando se conduce rápidamente a su estado “final” de cuasiequilibrio local. Por lo tanto el estado de FO observado en las entradas de colada – estas son donde realmente entra la masa fundida a la pieza – se determina completamente por su historia de flujo a través del sistema de canal de colada y en gran medida independiente de su estado inicial asumido en la entrada. Por otra parte un análisis de la estructura analítica de la FTE así como de su cinemática (es decir los posibles estados en el espacio de fase) muestra que la suposición de una orientación aleatoria en la entrada es la elección óptima en vista de una adaptación subsecuente por el campo de flujo en el sistema de canal de colada cerca de la entrada, debido a que un estado de orientación aleatoria da un acoplamiento completo de todas las componentes del gradiente de velocidad.

El proceso de división de operador descrito anteriormente se describe en forma simplificada con referencia a la Fig.

11:

Paso 1: establecer las condiciones iniciales y las condiciones de frontera en las celdas recién llenas,

Paso 2: mover la orientación de las fibras de acuerdo con el campo de flujo, y

Paso 3: calcular la rotación de las fibras de acuerdo con el gradiente de velocidad.

7. Enfoque numérico para el sistema de ODE del “paso de rotación”

Se usa un método numérico para integrar el sistema acoplado de ODE (18b), el que ha de resolverse como una

parte del enfoque de “división de operador” para la FTE cerrada. Este método se ha desarrollado específicamente

para la implementación del módulo de FO del método de acuerdo con la modalidad preferida. El mismo abarca varios aspectos que son relevantes dentro del contexto de las aplicaciones industriales típicas de una simulación de moldeado por inyección tridimensional para materiales termoplásticos reforzados con fibras cortas, tales como:

La formulación de la FTE usando el conjunto completo de 6 elementos independientes de la matriz de FO (opuesto al “procedimiento estándar” de usar la condición de traza para eliminar uno de los elementos de la diagonal y reducir en uno el número de variables dependientes de la FTE);

La adición de un término de control de “penalización” al r.h.s. de la FTE que transforma el hiperplano definido por la condición de traza en una variedad integral estable de la FTE;

La explotación del comportamiento de escalado específico de la función del r.h.s. de la FTE con respecto a las componentes del gradiente de velocidad para construir un método de integración el que selecciona el esquema de tiempo de integración de acuerdo con el tamaño de la tasa de cizallamiento local (máx. componente del gradiente de velocidad).

La implementación de una “clausura híbrida estabilizada” que contiene un factor de orientación escalar que se confina al intervalo [0,1].

La implementación de un procedimiento eficiente para la evaluación de la función especial del r.h.s de la FTE con clausura híbrida con una cantidad mínima de operaciones.

Estos aspectos se discuten en detalle más adelante.

7.1 Estabilización de la clausura híbrida

Las aproximaciones de clausura híbridas como se definen por las ecs. (6) y (6a/b) experimentan problemas de estabilidad especialmente en casos donde el gradiente de velocidad tiene una estructura complicada - es decir no es de tipo simple de cizallamiento/estiramiento sino que refleja un patrón complejo de flujo tridimensional - y el tamaño del paso de tiempo (como se determina por el solucionador de flujo) se hace más bien grande. La desviación del factor de orientación escalar fs (â(2)) a valores negativos es una fuente primaria de estas inestabilidades: una vez que este factor se hace negativo, la solución numérica se hace rápidamente inestable y diverge exponencialmente hacia valores muy lejos del conjunto de espacio de fase admisible para las matrices de FO.

Por lo tanto un control de los valores de fs(â(2)) mejora la estabilidad numérica del procedimiento de cálculo de FO.

Como el determinante de una matriz de FO está siempre confinado al intervalo

los valores admisibles teóricamente del factor de orientación como se definen por la ec. (6c) también se restringen mediante 0 ≤ fs(â(2)) ≤ 1. Por lo tanto, la definición estándar (6) de la clausura híbrida se modifica de la siguiente manera:

Las definiciones de los términos lineal y cuadrático de la clausura como se dan por las ecs. (6a/b) se mantienen igual, y el factor de orientación fs(â(2)) como se define por la ec. (6c) se sustituye por el factor de orientación confinado s(â(2)) como se define en la ec. (19) anteriormente. Esta definición de s(â(2)) implica que s(â(2)) = fs(â(2)) si el término 1-27 det(â(2)) se evalúa a un valor dentro del intervalo [0,1] y de lo contrario toma 0 ó 1 como valor mínimo

o máximo.

La definición (19) se denota como aproximación de clausura híbrida estabilizada. Los experimentos numéricos han mostrado que el confinamiento de los valores de s(â(2)) al intervalo [0,1] evita el desarrollo de inestabilidades severas en los ensayo ejemplares considerados. La clausura híbrida estabilizada se ha probado exitosamente en una variedad de simulación de moldeo por inyección tridimensional.

7.2 Reducción del número de variables mediante la condición de traza

El tensor de FO de 2do orden a^^((2) ), el que es la variable dependiente de la FTE cerrada como se presenta en la

ec. (16), es una matriz 3×3 simétrica. Dado que sus elementos

constituyen a priori un conjunto de seis variables mutuamente independientes, la FTE es un sistema acoplado de seis PDE. Aunque la estructura algebraica de la FTE admite formalmente matrices 3×3 simétricas arbitrarias, las matrices tienen que cumplir un número de condiciones adicionales formuladas como restricciones en sus invariantes que las califican como matrices de FO.

La condición de traza

, la que es la más simple de estas condiciones invariantes, abre una

posibilidad obvia para eliminar uno de los elementos de la diagonal

de la matriz de FO y de esta manera reducir en uno el número de variables. La mayoría (si no todos) de los artículos de investigación publicados que investigan el flujo de una suspensión de fibras bajo más bien situaciones académicas de flujo (como por ejemplo los flujos simples de cizallamiento o de estiramiento) usan este enfoque tomando una elección fija para el elemento de la diagonal específico a ser eliminado. Esta eliminación se logra entonces sustituyendo por ejemplo

en el r.h.s. de la ec. (16), el que de esta manera se hace independiente de

Como la FTE es formalmente compatible con la condición de traza para todas las aproximaciones de clausuras razonables por medio de las ecs. (14) o (15a), tanto la condición de traza y la PDE para el elemento de la diagonal

se cumplen idénticamente (es decir como identidades algebraicas formales), y esto es suficiente para considerar en lo adelante sólo el conjunto restante de las cinco componentes de la matriz de FO y el sistema acoplado de sus PDE.

Este procedimiento es puramente formal. Si la elección del elemento de la diagonal a ser eliminado de esta manera

.U
ٔ׏ del flujo y de su gradiente Ues buena depende de las propiedades espaciales del campo de velocidad local

Para los tipos simples de flujos que sólo visualizan las variaciones planas de la velocidad de flujo (es decir en el

plano x1-x2) la eliminación de como se explicó anteriormente es una elección razonable, como en este caso la dinámica de rotación dominante de las fibras se restringe al plano de flujo, como la forma especial del gradiente de velocidad induce sólo un “acoplamiento débil” a la dirección ortogonal de x3 en este caso. Por este argumento

también es plausible que elegir ya sea

o

para la eliminación conduciría más probablemente a problemas numéricos (por ejemplo concernientes a la estabilidad de la solución) y por lo tanto a resultados menos satisfactorios. Sin embargo, la situación es muy diferente en las simulaciones de moldeo por inyección tridimensional

pueden variar U
ٔ ׏ y su gradiente Utanto la velocidad de flujo de geometrías complejas de molde, ya que dentro

de una manera compleja en el espacio así como en el tiempo y por lo tanto no puede suponerse que tienen alguna forma simple específicamente, de manera que una elección fija (como la anterior) es muy probable que no sea óptima. Por lo tanto la eliminación de uno de los elementos de la diagonal de la matriz de FO por medio de la condición de traza se considera que es inapropiada para la simulación de flujos de suspensión industriales.

Otro argumento - relacionado con la estabilidad numérica - contra la eliminación de uno de los elementos de la diagonal mediante la condición de traza unitaria se presenta en una serie [34-36] de artículos de Shampine, donde se investigan los sistemas de ODE con restricciones algebraicas especiales (“leyes de conservación”). En [34-36] se considera un sistema acoplado de ODE ∂/∂t c = F(c) donde el vector n-dimensional c = (c1,..., cn)T de “concentraciones de especies” se restringe a un hiperplano de dimensión n-1 por medio de la condición

de conservación de la masa. Aunque esta condición formalmente da la posibilidad de eliminar una de

las concentraciones, por ejemplo sustituyendo la relación

en la función F(...) en el r.h.s., y de esta manera cumplir la ley de conservación exactamente al mismo tiempo, Shampine argumenta que usar este “truco” dentro de un procedimiento de integración numérica conduce meramente a una acumulación del error numérico de las concentraciones c1,...,cn-1 las que se calculan mediante la integración numérica del sistema de ODE (reducido) dentro de la concentración cn restante obtenida algebraicamente mediante la ley de conservación Aunque la ley de conservación por construcción se satisface siempre exactamente, la solución numérica calculada de esta manera no se garantiza que sea exacta por ningún medio. De hecho se conoce que especialmente para sistemas de ODE rígidos la solución numérica puede apartarse arbitrariamente lejos de la solución real si se usa descuidadamente este “enfoque de eliminación directa”. De acuerdo con [34 -36] estos argumentos se mantienen válidos también en el

caso de las leyes de conservación en general lineales que involucran un vector w = (w1,..., wn)T de los factores de peso y conducen a las mismas conclusiones incluso si se consideran restricciones no lineales más generales de la forma g(c) = const. .

Como la ecuación de Folgar-Tucker aumentada por la restricción de “traza unitaria”, la que es una “ley de conservación” lineal, es un ejemplo específico de las ecuaciones matemáticas presentadas en [34-36], el análisis matemático de Shampine muestra que eliminar uno de los elementos de la diagonal mediante la condición de traza unitaria muy probablemente conducirá a inestabilidades numéricas en tanto que se introducen inexactitudes incontrolables en los elementos de la diagonal de la matriz de FO.

Debido a los argumentos dados anteriormente, el método de acuerdo con la modalidad preferida - desviándose del “enfoque estándar” de eliminar uno de los elementos de la diagonal mediante la condición de traza unitaria - usa todos los seis elementos de la matriz de FO y las ecuaciones del sistema de FTE para el cálculo de la FO.

7.3 Estabilización dinámica de la traza mediante un término de penalización

Si el sistema de FTE cerrada y su parte “paso de rotación” se consideran como dados por las ecs. (16) y (18b) y la

condición de traza no se usa para eliminar uno de los elementos de la diagonal de la matriz de FO del conjunto de variables dependientes, se necesita algún método adicional con el objetivo de mantener la traza de la matriz de FO calculada numéricamente tan cerca como sea posible de su valor deseado 1.

En los artículos [34-36] se recomienda evitar cualquier tipo de enfoque de eliminación (como se discutió en la sección anterior 7.2) y tratar el sistema completo mediante una combinación de métodos usuales de integración de ODE y algún procedimiento de proyección que tiene en cuenta (ya sea exacta o de manera aproximada) el cumplimiento de las leyes de conservación o restricciones algebraicas requeridas. Para cada paso de tiempo la solución numérica se calcula primero sin considerar las restricciones algebraicas y subsecuentemente se corrige mediante un mapeo de proyección sobre el punto más cercano de la variedad algebraica definida por la ecuación de restricción (por ejemplo un hiperplano en el caso de una restricción lineal).

El principio alternativo para este enfoque trata de tomar en cuenta las restricciones mediante la adición de términos de control apropiados al r.h.s. de las ecuaciones del modelo, la traza de la matriz de FO calculada numéricamente puede mantenerse tan cerca como sea posible de su valor deseado 1. La presencia de los términos de control da un cumplimiento aproximado de las restricciones o leyes de conservación sin ninguna medida adicional. En muchos casos puede mostrarse que un control “duro” resulta efectivamente en un tipo de proyección, mientras que un control “suave” permite pequeñas desviaciones de las restricciones requeridas.

En el método de acuerdo con la modalidad preferida el enfoque de control se usa especialmente para mantener la traza de la matriz de FO calculada numéricamente aproximadamente en su valor unitario requerido. La estabilización se basa en el concepto de la transformación de un hiperplano de 5 dimensiones definido por la condición de traza Tr(â) = 1 en una variedad integral estable del sistema de FTE cerrada mediante la adición de un término de control apropiado a la función del r.h.s. de la ec. (16) o (18b), respectivamente. Desde un punto de vista diferente esto es equivalente a la situación donde la Tr(â) = 1 se convierte en un punto fijo estable de la ecuación diferencial (12) de la traza.

Se ha determinado que los requisitos generales que se deben cumplir por cualquier término de control razonable son los siguientes:

El término de control tiene que anularse idénticamente si la condición de traza ya se cumple.

Si la solución numérica no cumple la condición de traza, el término de control debe conducir a la solución numérica de regreso al hiperplano Tr(â) =1.

Aparte de estos dos requisitos principales, un término de control razonable debería ser compatible con el comportamiento de escalado de la función en el r.h.s. de la FTE con respecto al reescalado del gradiente de velocidad (ver la sección 7.4 sobre la integración en el tiempo controlada por flujo para más detalles sobre este punto), y debería contener un parámetro de ajuste que permite un ajuste de un control “suave” o “duro”, el último que se asemeja efectivamente a una proyección en el límite del control “infinitamente duro”.

La elección especial de un término de control apropiado ha de ser adaptada a la elección de una aproximación de clausura específica, en este caso de la aproximación de clausura híbrida (tanto la variante estándar y estabilizada).

(2) (...) la DE para la traza se obtiene mediante las identidades

En el caso especial de la clausura híbrida (estándar) la traza de la función del r.h.s. puede evaluarse a lo que da la forma simple (15a) de la DE de la traza con el prefactor variable øa

fs(â(2))

φa como dado por la (15b). En el caso de la clausura híbrida estabilizada el factor de orientación escalar de

se sustituye por su variante restringida (estabilizada) definida en la ec. (19), de manera que formalmente se obtiene exactamente la misma DE para la traza con un prefactor ǿa que contiene s(â(2)) en lugar de fs(â(2)).

La elección específica del término de penalización a ser adicionado a la función

del r.h.s. en el caso especial de la clausura híbrida (estabilizada) es:

El parámetro α se requiere que sea positivo, pero no necesariamente constante en el tiempo. La matriz 3×3

se requiere que sea simétrica con traza unitaria, pero arbitraria por lo demás. La matriz

define la “dirección” del término de control (20). Puede ser constante o variable en el tiempo. Lo último incluye la posibilidad para definir como una función de â(2) . Diferentes elecciones de la matriz

corresponden a diferentes variantes del término de control La variante usada en el método de acuerdo con la modalidad preferida usa

=

⅓ 1^ (es decir bij = ⅓ δij), lo que corresponde a un término de control en la dirección ortogonal al hiperplano Tr(â) = I. Una variante alternativa se da mediante la elección de

=

â(2)/Tr(â(2)) (ver [18]). de la adición del término de control (20) la variante generalizada

de la FTE cerrada con clausura híbrida estabilizada debe también considerarse. Similarmente el término de control entra en el r.h.s. de la ODE (18b), la que de esta manera asume la forma de la ec. (21) con la derivada material D/Dt sustituida por la derivada parcial ∂/∂t. Como la traza del término de control (20) se evalúa a Tr( hyb(...))= (α-ǿa)·[1-Tr(

(2))], la forma específica del término de control se adapta para modificar (15a/b) de acuerdo con

Como α > 0, el valor Tr(â(2)) =1 se corresponde con un punto fijo estable de la DE (22), y el hiperplano correspondiente por lo tanto se convierte en una variedad integral estable de la forma modificada (“controlada”) (21) de la FTE con clausura híbrida. Como una consecuencia, la traza de cualquier solución numérica de (21) se mantiene cerca del valor requerido 1 por el término de control (20), dado que este término obliga a todas las soluciones que tienen Tr(â(2)) > 1 en la dirección de

hacia valores menores de traza, y aquellas soluciones que tienen Tr(â(2)) < 1 son correspondientemente obligadas hacia valores mayores de traza (también en la dirección de

).

La “fuerza” del término de control (20) es ajustable mediante el ajuste apropiado del parámetro α: mientras que valores pequeños del parámetro α conducen a correcciones relativamente “suaves” de los valores de la traza hacia 1, valores grandes de α resultan en un fuerte impulso (cerca del efecto de una proyección instantánea) hacia el hiperplano α resultan en un fuerte impulso (cerca del efecto de una proyección instantánea) hacia el hiperplano Tr(â(2)) =1. En la práctica valores grandes de α pueden provocar problemas cuando se usan con métodos de integración de ODE explícitos, dado que en este caso las correcciones grandes de la traza empujan la solución numérica a cualquier lado del hiperplano de una manera alternante y de esa manera inducen oscilaciones espurias. Se ha mostrado [18] que una elección de α dentro de un relativamente amplio rango de valores intermedios conduce a resultados numéricos apropiados sin artefactos indeseados inducidos por el término de control.

Para cualquier elección específica de una aproximación de clausura que resulta en una función correspondiente en el r.h.s. de la FTE cerrada, la definición genérica

del término de control conduce idénticamente a la traza de DE (22) y por lo tanto a una estabilización del hiperplano Tr(â(2)) =1 como se describió anteriormente para el caso especial de la clausura híbrida (estabilizada). Si permitimos que la matriz 3×3 simétrica

de traza unitaria sea una función arbitraria de â(2), la ec. (20a) constituye la forma más general del término control.

7.4 Tiempo de integración controlado por flujo

Usando el enfoque "división de operador” ha de calcularse una solución numérica del sistema de ODE (18b) “paso de rotación” en el intervalo de tiempo [tn, tn+1] con el tamaño del paso global Δtn+1 = tn+1 -tn proporcionado por el solucionador de flujo dentro de todas las celdas de cálculo (etiquetadas por el vector espacial rm) del dominio lleno

)entra como un conjunto de +1nt,mr(U
ٔ ׏ . En cada celda las componentes del gradiente de velocidad local +1)n(Ω

coeficientes externos en el r.h.s. de (18b). Si éstos son grandes, puede esperarse que los valores de los elementos de la matriz de FO cambien significativamente durante el intervalo de tiempo, mientras que un gradiente de velocidad muy pequeño puede resultar en cambios despreciables de los elementos de la matriz de FO desde sus valores iniciales como se proporcionaron por el paso de cálculo precedente.

En todos los casos prácticos de una simulación de moldeo por inyección el gradiente de velocidad varía significativamente sobre el dominio lleno dentro de la cavidad del molde durante la fase de llenado. Mientras que los valores de la tasa de cizallamiento son típicamente altos cerca de las paredes de la cavidad, se han observado valores considerablemente más bajos en la región central de la separación entre las paredes. Esto resulta en una orientación más bien rápida de las fibras sujetas al flujo de cizallamiento fuerte cerca de las paredes en la dirección de la velocidad de flujo, mientras que las fibras en la región central toman un tiempo notablemente más largo para

un cambio en su estado de orientación, lo que conduce a la bien conocida “estructura de sándwich” (que consiste en una región central de “núcleo” con un grado relativamente bajo de la orientación de las fibras entre las regiones

altamente orientadas cerca de las paredes) lo que se observa típicamente en las piezas moldeadas por inyección fabricadas de termoplásticos reforzados con fibras cortas.

Los esquemas numéricos para la integración de los sistemas de ODE tienen en cuenta la “fuerza” de la función del

r.h.s. mediante una elección apropiada (preferiblemente adaptativa) del tamaño del paso. De esta manera se intenta evitar las inexactitudes que ocurren debido a un espaciado grueso a través de intervalos de valores grandes o a un cambio fuerte en la función del r.h.s., así como también tomar pasos muy pequeños a través de regiones donde los valores (absolutos) de la función del r.h.s. son más bien pequeños. Mediante una elección apropiada del tamaño del paso así como del esquema de integración es posible controlar el error numérico de la integración de la ODE hasta un límite deseado.

En el caso especial de un sistema ODE de “paso de rotación” debe proporcionarse una integración numérica sobre el intervalo de tiempo [tn, tn+1] con tamaño de paso global Δtn+1 así como un error numérico global εFO que debe ser uniforme sobre todas las celdas del dominio de cálculo, a pesar de la variación de la función en el r.h.s. debida a la variación del gradiente de velocidad Para este propósito se ha diseñado un esquema numérico especial que logra este objetivo con un número mínimo de evaluaciones de la función en el r.h.s. mediante la elección de

La "regla de integración" (Euler, punto medio o Runge Kutta) 4to orden y el tamaño del paso local(adimensional) así como el número local de pasos de este tamaño con la regla escogida

adaptada a la “fuerza” de la función local del r.h.s. como se determine por los valores del gradiente de velocidad local. El esquema es diferente de los enfoques conocidos de “control adaptativo del tamaño del paso” para la

integración numérica de ODE, ya que aprovecha el comportamiento de escalado de la función en el r.h.s de las ecs.

(16) o (18b) con respecto a las componentes del gradiente de velocidad y se adapta específicamente para el sistema de la FTE. El procedimiento de elección de la regla de integración y del número y tamaño de los pasos de tiempo locales inherentes al esquema es un enfoque heurístico que se basa en el error de discretización teórico de las reglas de integración usadas. Ambos aspectos del esquema se explican en detalle más adelante.

La forma específica de la función F(2) (...) en el r.h.s. de las ecs. (17b) o (18b) puede obtenerse a partir de la ec. (16) como

U explícitamentese tienen que
ٔ ׏ Con el objetivo de resolver todas las dependenciasdel gradiente de velocidad

sustituir las definiciones (1a/b) del tensor de gradiente de velocidad efectiva y el de tasa de cizallamiento, es decir

así como la identidad + T = 2λ

y las fórmulas

que definen el parámetro de difusión de rotación y la tasa de cizallamiento efectivo. Para las siguientes consideraciones será útil suprimir las dependencias indirectas de la función (2)( ...) de la matriz de FO mediante la función de clausura Â(4)[â(2)] y del gradiente de velocidad mediante la tasa de cizallamiento efectivo así como la dependencia de los parámetros constantes del modelo C1 y λ y por lo tanto reescribir la expresión algebraica de la función del r.h.s. (como ya se hizo para el r.h.s. de la ec. (17c)) en una forma simplificada como

(2)Uen la lista de argumentos deF
ٔ dejando sólo la dependencia de la matriz de â(2)y el gradiente de velocidad

(...). Una reformulación equivalente del sistema de ODE (18b) (que corresponde a (17c)) usando la función del r.h.s.

(23) puede escribirse en la forma

U por algún factor σ > 0 induce las
ٔ׏ U)ij de los componentes de
ٔ׏ (σհU)ij
ٔ׏ (Una multiplicación

multiplicaciones correspondientes

հ σ

հ

,

հ σ

, γeff հ σγeff que, debido a la estructura algebraica especial

del r.h.s. de (23), conduce a la identidad algebraica

que expresa el comportamiento de escalado de la función (2)(...) con respecto a su 2do argumento.

Como se mencionó antes, la tarea general es una integración numérica del sistema acoplado (18c) de ODE en el intervalo de tiempo [tn, tn+1] con un error numérico global εFO que es uniforme sobre todas las celdas del dominio de celdas llenas. Esto puede lograrse explotando la propiedad especial (24) de la función

(2)(...) del r.h.s. para un

escalado del gradiente de velocidad, la función del r.h.s. y el sistema de ODE (18c) de la siguiente manera:

Determinar el valor de la componente del gradiente de velocidad con el mayor valor absoluto, es decir:

depende del vector espacial max γ también el parámetro U
ٔ ׏ (Note que con

rm que etiqueta la celda de cálculo local, así como de la coordenada de tiempo tn o tn+1 en la que el gradiente de velocidad entra a la función del r.h.s.).

Introducir el gradiente de velocidad escalado

y la coordenada de tiempo escalada τ: = γmax t - note que ambas magnitudes son adimensionales, y que este escalado se aplica sólo localmente en la celda

usando las magnitudes escaladas y la ec. (24). (24).

La ODE original (18c) ha de ser integrada sobre el intervalo de tiempo “global” [tn, tn+1] con tamaño de paso “físico” Δtn+1 = tn+1 (medidos por ejemplo en [s]), los que son los mismos para todas las celdas de cálculo del dominio

,U
ٔ , ylas funciones del r.h.s. cuya “fuerza” varía de acuerdo con la variación del gradiente de velocidad V+1)n(Ω

-

tn

como se indica por los diferentes valores del parámetro γmax para las diferentes celdas en el dominio. Diferente de que, la ODE transformada (25) ha de integrarse sobre un intervalo de tiempo escalado “localmente” [τn, τn+1] de tamaño adimensional

pero con las funciones del r.h.s. que, como el gradiente de velocidad escalado

una “fuerza unitaria” aproximadamente uniforme debido a los siguientes hechos:

U entra en el r.h.s. de (25), tiene
׏ ׏

i. El gradiente de velocidad escalado

es por construcción una magnitud adimensional, con los valores U
׏ ׏

absolutos de sus componentes iguales o menores que 1.

ii.

Los valores de los elementos de la matriz de FO así como aquellos elementos del tensor de 4to orden suministrados por la función de clausura Â(4)[â(2)] también se limitan por 1.

ii.

Por lo tanto todas las componentes de la función F(2)(...) del r.h.s. son de orden O(1) cuando se evalúan

(2) (...) tienen valores absolutos comparables con γmax

,mientras que las mayores componentes deU
׏ ׏ con en su lugar. U
ٔ ׏ cuando se evalúan con

para medir la “fuerza” de
צ ...
צ Usando una norma matricial arbitraria

una forma matemática precisa:

(2) (â(2)

Como las funciones escaladas del r.h.s.

) de la ODE transformada (25) de acuerdo con (27) tienenU
׏ ׏ ,

“fuerza unitaria” uniforme para todas las celdas en el dominio Ω(n+1), se puede lograr la integración numérica con

error global uniforme εFO sobre intervalos de tamaño variable localmente mediante la elección de

i. una “regla de integración" con error de discretización suficientemente pequeño así como

ii. un número adecuado de subpasos para cubrir el intervalo de tiempo escalado.

El esquema de integración usado en el método de acuerdo con la modalidad preferida usa un conjunto de reglas de integración simples que pertenecen a la clase de métodos de un paso (ver las secciones 16.1 de [38] y 7.2.1-7.2.3 de [39]) incluidas en la presente como referencia. Las reglas de integración específicas usadas dentro del módulo de FO del método de acuerdo con la modalidad preferida son:

1er

la regla Euler hacia adelante simple, a que es un método de precisión de orden y requiere 1 evaluación de la función del r.h.s. por paso,,

la regla del "punto medio" o Runge-Kutta ,2do orden (RK2) la que es un método de precisión de 2do orden y requiere 2 evaluaciones de la función del r.h.s por paso, y

la regla de Runge-Kutta de 4to orden (RK4) la que es un método de precisión de 4to orden y requiere 4 evaluaciones de la función del r.h.s por paso

Las siguientes propiedades de estas reglas de integración que dependen de su orden p se usan en la construcción del esquema especial usado en una modalidad de la presente invención:

a) Un método de orden p requiere p evaluaciones de la función del r.h.s. por subpaso. b) Si un ODE (sistema) se integra numéricamente sobre un intervalo de tamaño total Δτ con N subpasos equidistantes de tamaño h = Δτ/N, el error total εtot acumulado en el subpaso final del intervalo puede estimarse como: εtot ~ Δτ·hp. c) Si se requiere que este error total sea menor que un valor umbral preseleccionado εmax, el tamaño de los subpasos se acota por h < (εmax/Δτ)1/p. Similarmente el número de subpasos tiene que ser mayor que N > Δτ·(Δτ/εmax)1/p. d) La integración sobre el intervalo total con un error total menor que εmax puede realizarse tomando N

· Np)1/(p+1)

subpasos con un método de orden p siempre que el tamaño del intervalo se limite por Δτ < (εmax.

Esto implica que sea suficiente un subpaso único (es decir N =1) de tamaño h = Δτ si se cumple

Tomando esto en cuenta, se puede construir un esquema de integrac ión “híbrido” que da una integración numérica

de la ec. (25) sobre el intervalo de tiempo escalado de tamaño

definido en (26) con un error global menor que εmax = εFO mediante la elección de

(i)

ya sea un paso único con una de las tres reglas de integración o

(ii) múltiples subpasos de tamaño apropiado usando la regla de RK4

donde la elección depende del tamaño relativo de comparado con el umbral de error requerido εFO (ver más adelante). Como este límite de error es el mismo para todas las celdas de cálculo del dominio Ω (n+1), la integración

se realiza con un error uniforme εFO independiente de los valores locales de (o similarmente

El método propuesto elige específicamente la regla de integración que alcanza el límite de error requerido con la cantidad mínima de evaluaciones de la función del r.h.s. como se estime de acuerdo con las relaciones presentadas en los puntos a) hasta d) anteriormente. Aunque el tamaño del paso de tiempo Δtn+1 se determina por el flujo es típicamente más bien pequeño (Δtn+1 ≈10-2 ...10-4 s) comparado con el tiempo total de llenado del molde (el que es

del orden de segundos), la magnitud adimensional puede ser del orden O(1), dado que los valores de γmax pueden llegar a ser más bien grandes (del orden de 10...100s-1) debido al tamaño pequeño de la separación y la gran viscosidad de la masa fundida de polímero. Como los valores actualizados de la matriz de FO son también del

, de manera que -4 ... 10 -210
؄ FO(1), una elección razonable del límite del error requerido varía en el rango εOorden

puede asumirse un tamaño de subpaso de h ≤ 0.1 en aplicaciones típicas. Puede mostrarse que para tamaños de paso h < ½ un paso único de tamaño h con la regla RK4es más preciso que dos pasos de la mitad del tamaño h/2 usando la regla del "punto medio" (RK2) o cuatro pasos de Euler de un cuarto del tamaño h/4,aunque el esfuerzo numérico (4 evaluaciones de la función del r.h.s.) es el mismo para todas las variantes.

La elección de la regla de integración así como el número de (sub)pasos se basa en el tamaño de

relativo a la

{1, 2, 4} que define el orden de
א secuencia de los límites del error calculados para los valores p

precisión de las varias reglas de integración. Como 0 < εFO < 1 (de hecho típicamente se cumple que εFO << 1 el

tamaño de los límites del error crece monótonamente con el valor creciente del parámetro de orden p (siempre tenemos 0 < εp < εp+1 < 1), de manera que los límites relevantes del error para la formulación del algoritmo dado más adelante están siempre ordenados de acuerdo con 0 < εFO < ε1 < ε2 < ε4 < 1. (Para una elección típica de por ejemplo εFO = 0.001 se obtienen ε1 ≈ 0.032, ε2 = 0.1 y ε4 ≈ 0.25 los valores numéricos.) Similarmente existe un valor

de umbral mínimo εmin para el tamaño del intervalo de tiempo escalado adimensional por debajo del que

incluso un paso único de Euler de tamaño meramente da una actualización de los elementos de la matriz de FO al valor previo que es despreciablemente pequeño. Esto ocurre típicamente en todas las celdas del dominio donde el gradiente de velocidad es muy pequeño, por ejemplo en la región central de núcleo en el caso de un tamaño de separación relativamente grande entre dos paredes adyacentes. Un umbral mínimo de por ejemplo εmin = 10-6 es una elección razonable para aplicaciones típicas de llenado de molde.

Tomando en cuenta todos los aspectos como se discutieron anteriormente, el algoritmo para realizar dentro de cada celda del dominio e ilustrado en el diagrama de flujo de la Fig. 12, puede formularse de la siguiente manera:

y calcular gradiente velocidad escalado ∞
צ U
ٔ׏צ :=maxγPrimero evaluar 1.

elde

y el tamaño de intervalo escalado

2. Chequear si en el paso 12.1.

a.

En caso afirmativo, entonces saltar la integración de la ODE (es decir no hacer nada), mantener el valor anterior de la matriz de FO en la celda como el nuevo valor (actualizado) y salir del procedimiento.

b.

Si no, entonces ir al paso 12.2.

3. Chequear si

a. En caso afirmativo, entonces actualizar el valor anterior de la matriz de FO en la celda mediante un

para la evaluación de la función del r.h.s. y salir U
׏ ׏ paso único de Euler de tamañousando

del procedimiento.

b. Si no, entonces ir al paso 12.3.

4. Chequear si

a. En caso afirmativo, entonces actualizar el valor previo de la matriz de FO en la celda mediante un paso

U ysalir del procedimiento.
׏ ׏ usando

único de “punto medio” (RK2) de tamaño

b. Si no, entonces ir al paso 12.4.

5. Chequear si

a. En caso afirmativo, entonces actualizar el valor previo de la matriz de FO en la celda mediante un paso

U ysalir del procedimiento.
׏ ׏ usando

único de "Runge-Kutta de 4to orden (RK4)" de tamaño

b. Si no, entonces ir al paso 12.5.

6. En el caso general la actualización de la matriz de FO se calcula realizando varios pasos de RK4 con un tamaño de paso apropiado.

a.

El menor número de pasos para lograr la precisión requerida εFO es el menor entero N que satisface la desigualdad N > Δτ · (Δτ/εFO)1/p. Usar las funciones int(...), que devuelve la parte entera de un número de punto flotante, y max(..., ...), que devuelve el mayor de dos números, el entero N ≥ 1 puede calcularse mediante la fórmula:

b.

El tamaño de paso requerido se calcula entonces mediante

c.

Una vez que N y h se calculan, actualizar el valor previo de la matriz de FO en la celda mediante N

en el r.h.s.) y salir del procedimiento U
ٔ׏ (usando hcon tamaño de paso RK4pasos de la regla

tiempo [tn,tn+1] con un error global uniforme menor que εFO y un número mínimo de evaluaciones de la función del

r.h.s. de (18b). La aplicación de las reglas de integración a la versión escalada (25) de (18b) se hace “en el lugar” usando el gradiente de velocidad escalado en la evaluación de la función del r.h.s. y el cálculo del tamaño del paso h

basado en el tamaño de la longitud del intervalo escalado

En aplicaciones típicas de este algoritmo (con parámetros de límite de error εFO =0.001 y εmin =10-6) se observa que en la mayoría de las celdas (es decir alrededor del 80%) en el dominio Ωn+1 la actualización de "Paso de Rotación" de la matriz de FO se realiza mediante un único paso de Euler, un menor número de celdas son meramente

“pasadas por alto” debido a que es decir valores pequeños de γmax ), y un número razonable de celdas, las que presumiblemente se localizan cerca de las paredes de la cavidad del molde y por lo tanto tienen un alto valor de tasa de cizallamiento, se actualizan con uno o varios (típicamente no más de 20) pasos de " Runge-Kutta de 4to orden"

7.5 Integración de la ODE “Paso de Rotación” con estabilización dinámica de la traza

Si la función del r.h.s. contiene un término de control - ya sea del tipo especial (20) en el caso de la clausura híbrida

o del tipo general como se da por la ec. (20a) - para lograr la estabilización dinámica de la traza, el comportamiento de escalado de la función del r.h.s. sigue siendo compatible con el algoritmo de integración “híbrido” descrito en esta sección siempre que el parámetro de control de las ecs. (20) o (20a) se defina como α = α0 · γmax con un parámetro de control adimensional de tamaño típico α0 ~ O(1).

Ilustramos esto mediante la observación del término de control en su forma general (20a). Usando la propiedad de

escalado (24) de la función del r.h.s. (...) que contiene cualquier aproximación de clausura, podemos reescribir su traza de acuerdo con la secuencia de identidades

y, usando el factor de control α = α0 · γmax, obtenemos la expresión equivalente

para el término de control (20a). La expresión en el r.h.s. contiene ahora γmax como un factor separado y un término

, pero que U
ٔ Uen lugar de V
׏ ׏ (...) que se evalúa con el gradiente de velocidad escalado que depende de

es por lo demás formalmente idéntico a (20a). Esto puede expresarse por la identidad

la que también muestra que el término de control en su forma general (20a) tiene exactamente el mismo comportamiento de escalado que la función del r.h.s. “no controlada”

(...) con respecto al escalado del gradiente de velocidad mediante γmax si el parámetro de control se elige como α = α0 · γmax.

En el caso especial de la clausura híbrida (estabilizada) puede extraerse γmax del prefactor

a de la misma manera reescribiendo el término (20) como

׏ en lugar de U
ٔ׏ usando aǿ se calcula mediante la misma fórmula que maxγ/aǿ= aødonde el prefactor “reescalado” . Estomuestra explícitamenteque (28) es válida también para la forma especial del término de control asumido U
ٔ

en el caso de clausura híbrida

Las consideraciones anteriores muestran que el algoritmo de integración “híbrido” para el sistema de ODE “Paso de Rotación” pueden aplicarse sin cambio también en la presencia de un término de control, siempre que α0 se use como parámetro de control y los términos que dependen de la función del r.h.s. se evalúen con el gradiente de velocidad escalado. La versión escalada de la ODE correspondiente a la ec. (21) se obtiene simplemente añadiendo el término de control escalado al r.h.s. de la ec. (25), es decir la ODE

ha de considerarse en lugar de (25) como una base para los pasos de tiempo discreto del algoritmo de integración.

El esquema de integración en el tiempo implementado en el módulo de FO del método de acuerdo con la modalidad

preferida realiza la integración numérica del sistema de ODE “Paso de Rotación” mediante la aplicación del algoritmo de integración “híbrido” como se define en esta sección al sistema de ODE (29) que incluye la

estabilización dinámica de la traza por medio del término de control (20).

7.6 Evaluación eficiente de la FTE con la clausura híbrida estabilizada

La evaluación de la función del r.h.s. es la parte más costosa dentro del procedimiento de cálculo de FO, por lo tanto el método evalúa en la manera más eficiente.

El aspecto de la eficiencia se ha abordado en parte por el enfoque algorítmico de la “integración en el tiempo controlada por el flujo”(FCTI) como se propuso en la sección 7.4. Usando este enfoque, se puede lograr una precisión uniforme con un mínimo de evaluaciones de la función del r.h.s. de la FTE (posiblemente incluyendo un término de control).

El factor más importante el que directamente afecta el coste computacional del cálculo de FO es la elección de la aproximación de clausura. La clausura híbrida estabilizada da una precisión razonable a bajo costo computacional y se usa por el método.

Reformulación algebraica de la FTE con clausura híbrida

El primer paso para mejorar la eficiencia es la sustitución de la definición algebraica (6) de la clausura híbrida en la función del r.h.s. (23) de la FTE que da (después de un par de transformaciones algebraicas) la función del r.h.s. en una forma algebraica específica:

â(2)

La suma de los dos primeros términos del r.h.s. de (30) es formalmente idéntica a la expresión ·

+

T · â(2),

pero contiene la matriz

en lugar de la matriz

del gradiente de velocidad efectiva dada por 1 (a). Los tres términos restantes en el r.h.s. tienen todos la forma øxX de un producto de un prefactor escalar øx que multiplica a una cantidad matricial . . Los prefactores se dan por el conjunto de fórmulas

Las fórmulas (32) son válidas en el caso de cualesquiera campos de velocidad de flujo compresible o incompresible. Aquí se hace énfasis en que aunque teóricamente div U = 0 debe cumplirse cuando se asume que el flujo es incompresible, esto no se usa explícitamente, sino más bien todos los términos que contienen Tr ( ) = div U se mantienen en las fórmulas (32). Si Tn(

M) ≠ 0 debido a errores numéricos (lo que inevitablemente ocurre durante el cálculo del flujo), se encuentra que (32) aún conduce a un comportamiento estable, mientras que ocurren inestabilidades en el caso que los términos proporcionales a Tr(

M) se asume a priori que sean cero en las fórmulas de los prefactores (32).

El uso de la función del r.h.s. en la forma (30) con las magnitudes (31) y (32) es el paso básico hacia la evaluación eficiente, dado que la estructura algebraica de la función del r.h.s. (30) se organiza de una manera tal que una variedad de cálculos se tiene que hacer sólo una vez (cuando se calculan los prefactores), lo que ahorra muchas operaciones.

Varias operaciones algebraicas surgen dentro de las subexpresiones comunes que aparecen frecuentemente (por ejemplo los términos 2Dr o (1- fs)/7 en (31) y (32)), y calcularlos sólo una vez y guardar sus valores en variables auxiliares para su uso posterior reduce significativamente el trabajo computacional. Unas pocas operaciones más se ahorran mediante el siguiente paso de redefinición de la matriz

de acuerdo con

y usar la matriz

a en lugar de

para la evaluación de la función del r.h.s.. Haciendo uso de la identidad

algebraica

el cálculo explícito del término Ԅa â(2) se omite, y la función del r.h.s. se redefine de acuerdo con

Las operaciones adicionales para calcular

a mediante la ec.. (31a) son (i) 1 multiplicación para obtener ½øa y (ii) 3 "adiciones” para sustraer ½øa de los elementos de la diagonal de la matriz . El número de operaciones ahorradas asciende a una multiplicación de todas las componentes de la matriz de FO por el prefactor øa y las operaciones de sustracción para - øa â(2) . Por lo tanto el número total de operaciones es menor si las ecs. (30a) & (31a) se usan para la evaluación de la función del r.h.s. en lugar de (30) & (31).

En el caso de la clausura híbrida estabilizada el factor de orientación escalar restringido fs se calcula de acuerdo con (19). Sustituir

por fs en (32) resulta entonces en “prefactores estabilizados” ǿa, ǿ1 y ǿg sin afectar el costo computacional.

La adición de un término de control (ver la ec. (20))

conduce a operaciones adicionales necesarias para la evaluación del r.h.s.. La estructura general de este término de control se da por

hyb(...) = ǿb

con un prefactor ǿb =(α-ǿa)·[1-Tr( (2))] y por lo tanto es de la misma forma que ǿx como se discutió anteriormente. El prefactor ǿb se absorbe dentro de uno de los prefactores ǿ1 o ǿa si cualquiera de

= ⅓

o = â(2)/Tr(â(2)) se eligen para la matriz de proyección. En el primer caso el término

se adiciona al prefactor 1, lo que conduce a una función del r.h.s. modificada

La segunda elección conduce a un término adicional φa = (α -ǿa) · [1-Tr( (2)]/Tr( (2)) que ha de que adicionarse al prefactor ǿa.

Algoritmo eficiente para la evaluación del r.h.s.

Haciendo uso de los pasos anteriores, se presenta un esquema que evalúa (30b) con un número mínimo de operaciones algebraicas de acuerdo con la siguiente secuencia de cálculos:

α,U) & Parámetros: λ , C1
׏ ׏ ,(2)Entrada: (â1.

Variables auxiliares: , a,

M, λ+, λ-, fs, s, ǿa, ǿ1, ǿg, φ1, ζ, ξ1,...,ξ6

2. Inicializar a con la matriz de gradiente de velocidad efectiva calculada mediante (1a):

3.

Calcular la matriz simétrica

M ← Na + aT y su traza ξ1 ← Tr(

M).

4.

Calcular ξ2 = 2Dr a partir de la contracción ξ ←

M :

M mediante

א , 0))fs← min (1, max ( fs),su versión con tilde (2)-27 det( ←1Calcular el factor de orientación escalar fs 5.

[0,1] y el término ξ3 ← (1-fs)/7.

M : â(2)

6. Calcular la contracción ξ4 ←

.

7.

Calcular los prefactores de acuerdo con (32) usando los valores almacenados en las variables auxiliares ξk mediante las fórmulas siguientes:

8.

Calcular el término φ1 ← ⅓(α-ǿa)·[1-Tr(â(2))], entonces calcular ǿ1 ← ǿ1 + ǿ1.

9. Calcular la matriz a usando las variables auxiliares ξ5 ← 2ξ3 y ξ6 ← ½ ǿa. El cálculo se supone que se hace mediante la siguiente secuencia de operaciones:

← â(2)

10.

Inicializar la matriz resultado del r.h.s. mediante el cálculo de

.

a + aT · â(2). Entonces

sucesivamente completar la evaluación de la función del r.h.s. de (30b):

(2) (â(2)

) ←U
׏ ׏ ,

11. Resultado: Reformulación de la FTE en forma vectorial usando "notación contraída" (CN)

Otro elemento en la construcción de un procedimiento de evaluación eficiente para la función del r.h.s. se basa en el hecho de que la matriz de FO â(2) así como la función del r.h.s. (2) (...) ambas son matrices simétricas y poseen sólo 6 elementos matriciales independientes. Por lo tanto no es práctico almacenarlos como matrices 3×3 generales y es ineficiente implementar las operaciones que han de realizarse para la evaluación de (2) (...) como operaciones sobre matrices. De acuerdo con la invención las matrices 3×3 simétricas se tratan como vectores de 6 componentes

{1,..., 6} ylos pares de
א ) de losíndices del vector µ ijes usando el siguiente esquema de identificación μ ↔ (real

índices (ij) = (ji) de las matrices simétricas conocido como "notación contraída" (CN, ver por ejemplo [22]) (En mecánica estructural, esto se conoce también como "notación de Voigt"):

µ

1 2 3 4 5 6

(ij)=(ji)

(11) (22) (33) (23)=(32) (13)=(31) (12)=(21)

Esta elección especial de un esquema μ ↔ (ij) es por supuesto sólo uno de los 6!=720 esquemas equivalentes. La tabla anterior muestra la elección usual que también se adopta en el módulo de FO de la modalidad preferida. (Variantes alternativas elijen usualmente una asignación diferente de los elementos de la matriz fuera de la diagonal

{4, 5, 6} .)
א para los índices vectoriales μ

Usar el esquema CN da un mapeo biyectivo aij = aji հ aµ de los elementos de una matriz 3×3 simétrica sobre las componentes de un 6-vector correspondiente (es decir un vector de (es decir un vector de R6) como se ilustra por De esta manera la matriz de FO, la matriz

M y la matriz

que almacena el resultado de la evaluación de la

función del r.h.s.

) se asigna a las magnitudes vectoriales correspondientes. De la misma manera lasU
׏ ׏ ,

componentes de un tensor de 4to orden Â(4) se asignan a los elementos Aµv de una matriz 6×6 si el tensor es

5 simétrico con respecto a los pares de índices primero y segundo (es decir si se cumple

). Si en adición a que el tensor es también simétrico con respecto a un intercambio (ij) ↔ (kl) de ambos pares de índices y de esta manera posee simetría ortotrópica (ver la ec. (11)), Aµv = Avµ se cumple, es decir la matriz

es simétrica

ella misma y posee elementos independientes. El número de elementos independientes puede reducirse por la presencia de simetrías adicionales y relaciones algebraicas entre las componentes del tensor. En el

10 caso de un tensor de 4to orden completamente simétrico existen las 6 identidades adicionales

que reduce el número de elementos independientes de a 15. Las condiciones de normalización

15 y la condición de traza dan un conjunto de relaciones algebraicas que reduce el número de elementos independientes de la matriz a un valor aún menor (ver [22] para más detalles). Introducir el enfoque CN conduce a la FTE en “forma vectorial” como se da por

Para una matriz arbitraria 3 × 3 la matriz

el mapeo â հ â· + T·â define un mapeo lineal en el espacio vectorial de seis dimensiones de las matrices 3 × 3 simétricas reales, el que puede escribirse formalmente como un producto matriz vector ·a en R6 descrito por una matriz 6 × 6 real , cuyos elementos son cero o dados por términos

25 algebraicos simples de los elementos de .

En este sentido el esquema CN da la identificación

Similarmente la operación de contracción Â(4):

M se define elemento a elemento mediante la fórmula

usando las asignaciones de índices s µ ↔ (ij), v ↔ (kl) como se da mediante el esquema CN. Esta operación también puede escribirse como un producto matriz vector ·GM, donde los elementos de la matriz  e relacionan con 40 los elementos de  mediante

45 Esto revela que la matriz  (a diferencia de Â) es no simétrica. La matriz  resulta de asignar los elementos del tensor Â(4)[â(2)], el que él mismo se calcula a partir de la matriz de FO mediante una aproximación de clausura, de acuerdo con el esquema CN (es decir

[a] representa la aproximación de clausura en CN.

Costo computacional del algoritmo propuesto

La implementación más eficiente de las operaciones de cálculo realizadas en los pasos 2 al 10 del procedimiento propuesto para una evaluación eficiente de la función del r.h.s. resulta de aplicar el esquema CN a este

de: número
ٔ de operaciones (#La siguiente tabla da una visión general sobre el númeroprocedimiento. & sustracciones) y llamadas a funciones que se necesitan: número de adiciones
ْ divisiones, #multiplicaciones, &

10 para realizar los pasos individuales del procedimiento usando el enfoque CN:

Paso

Cantidad ٔ# ْ# Función

2

λ±, a ← 17 8 -

3

M, Tr( M) - 8 -

4

ζ ← M : M, ξ2 10 5

5

det(â(2)), fs, s, ξ3 13 7 min(...) max(...)

6

ξ4 ← M : â(2) 7 5 -

7

ǿa, ǿl, Ԅg 6 4 -

8

φ1, ǿ1, ← Ԅ1+ φ1 2 5 -

9

ξ5, ξ6, a ← a -... 8 12 -

10

F ← â(2)· a + a T ·â(2) ← +... 33 30 -

Σ

(2)(â(2), )←U׏ٔ 96 84 3

Los puntos siguientes comentan sobre los aspectos de la implementación así como sobre el conteo de operaciones dado en la tabla para varios pasos del algoritmo presentado anteriormente:

Como al matriz a es no se almacena como una matriz 3×3. Un almacenamiento extra de

no es necesario si a se inicializa con el cálculo de . La matriz simétrica

M se necesita en forma matricial para el cálculo de a (→ paso 9), pero también aparece como un 6-vector en el paso final 10.

Las operaciones de almacenamiento y asignación no se cuentan, dado que su contribución al costo computacional total es despreciable.

El número de variables auxiliares se reduce mediante la reutilización de variables ya existentes una vez que sus valores no se necesitan más. Esto no se ha hecho en la versión del algoritmo propuesto anteriormente, dado que esto reduciría la claridad de la presentación del algoritmo.

La operación de contracción de un par de matrices simétricas (como ocurre en los pasos 4 y 6) se

:

-usando CN.. (En notación matricial se define mediante
ْ - y 5 operaciones
ٔ implementa con 7 operaciones

- y 5
ٔ El determinante de una matriz 3×3 simétrica real (→ paso 5) se calcula con 11 operaciones -usando CN.
ْ operaciones

-sise
ْ - y 12 operaciones
ٔ El ensamble de la matriz a realizado en el paso 9 requiere 6 operaciones

usa una variable auxiliar que almacena

M como una matriz 3×3 simétrica, además 2 multiplicaciones para calcular las variables ξ 5 y ξ6.

← â(2)

- usando CN.
ْ -y 21 operaciones
ٔ a + a T·â(2) requiere 27 operaciones·

La operación matricial

Cualquier operación del tipo ... + βî que adiciona (o sustrae) un múltiplo de la matriz unidad requiere 35 solamente tres adiciones de β a los elementos de la diagonal.

La última fila de la tabla muestra que una única evaluación de la función del r.h.s. se obtiene al costo computacional

de un númeroraíz cuadrada tpara calcular la-y 3 llamadas a función
ْ -84 operaciones 96de 96 operaciones

de punto flotante (doble precisión) así como el mínimo y máximo de un par de tales números. En computadoras con

40 hardware moderno el costo computacional de las operaciones de adición y multiplicación es aproximadamente el mismo, una llamada a una función min o a una función max cuesta aproximadamente de 1.5-2 operaciones, y el costo del cálculo de una raíz cuadrada asciende aproximadamente a 6-10 operaciones (es decir un promedio de 8 operaciones, en dependencia de la precisión requerida).

En total una única evaluación de (30b) se logra mediante el algoritmo propuesto anteriormente a un costo computacional de aproximadamente 190 operaciones. Este es el enfoque más eficiente para evaluar el r.h.s. de la FTE con clausura híbrida (estabilizada) incluyendo la estabilización dinámica de la traza.

Evaluación de la eficiencia computacional

Una evaluación de la eficiencia computacional del procedimiento propuesto resulta a partir de una comparación con un enfoque “estándar” que evalúa la función del r.h.s. como se da en forma vectorial mediante la ec. (35). Tal 10 enfoque tendría que ser adoptado por ejemplo si se pretende que el código tenga cierta estructura modular tal que

(i) todas las operaciones que son independientes de la aproximación de clausura se realicen en una rutina “genérica” la que toma la matriz 6×6 real  como entrada y (ii)los elementos de la matriz A se calculan en una rutina separada de acuerdo con la elección específica de una función de clausura Â(4)[â(2)] usando el esquema CN y (35b) como se explicó anteriormente. El costo computacional de tal enfoque “estándar” puede evaluarse mediante las

15 siguientes consideraciones:

El procedimiento “estándar” comenzaría con el cálculo de la matriz , el 6 vector

y el parámetro escalar 2Dr (ver pasos 2 al 4), el que en la suma requiere 46 operaciones además del cálculo de una raíz cuadrada.

20 La operación âհâ· + T·â que corresponde a ·a requiere 48 operaciones si el producto matriz vector se lleva a cabo analíticamente por adelantado y sólo se evalúan numéricamente las fórmulas resultantes para las seis componentes.

El cálculo del producto matriz vector · La adición de 2Dr ·1 y la sustracción de 6Dr ·a (usando 6Dr = ¾·2Dr) requieren otras 3 + 12 + 1 = 16

M y su sustracción del r.h.s. requiere 66 + 6 = 72 operaciones.

25 operaciones.

Contando 8 operaciones para el cálculo de la raíz cuadrada, el costo computacional de este enfoque para evaluar el

r.h.s. de (35) adiciona hasta este punto. 190 operaciones La inclusión de un término de control para la estabilización de la traza adiciona algunas operaciones más (ver el paso 8), pero aproximadamente la misma cantidad de

30 operaciones pudiera ahorrarse sustituyendo la sustracción explícita de 6Dr·del r.h.s. mediante una sustracción del término 3Dr =3/2·2Drde la diagonal de

Estos "costos fijos" son independientes de la aproximación de clausura requerida para calcular los elementos de la matriz

, el que ha de hacerse en un procedimiento separado.

En el caso de una clausura híbrida la matriz  tiene 21 elementos independientes, como la función tensor

definida por (6) y (6a/b/c) tiene simetría ortotrópica pero es no totalmente simétrica debido al término

cuadrático (6a). Trasladar la definición de a CN da el siguiente conjunto de expresiones las que han de usarse para calcular los elementos de la matriz Aµv:

Las matrices simétricas

y

ambas contienen las componentes de las expresiones de un tensor de 4to orden totalmente simétrico las que tienen que calcularse analíticamente por adelantado.

Aunque los elementos de matriz Dµv son constantes, los elementos de la matriz Eµv son expresiones lineales de las componentes del vector de aµ. El cálculo analítico de los elementos de las matrices Dµv y Eµv da el siguiente par de matrices:

Usando estas matrices, el costo computacional para ensamblar la matriz  se determina de acuerdo con el siguiente procedimiento:

El cálculo del factor de orientación escalar restringido s y la variable auxiliar ζ2se realizan en el paso 5 del algoritmo efectivo y requieren 23 operaciones (contando 3 operaciones para las llamadas anidadas a min(0, max(..., 1))). La variable auxiliar ζ1 =ζ2/5 se calcula mediante una división adicional, de manera que la cantidad de operaciones requeridas para calcular s y ζ1/2 es 24.

La matriz A se inicializa mediante Aµv ← s aµav, lo que toma 2 multiplicaciones para cada uno de los 21 elementos matriciales independientes del triángulo superior, ó 42 operaciones en total. (Los restantes elementos matriciales se asignan por simetría).

Los elementos de la matriz ζ1

se calculan mediante una única operación (es decir el cálculo de 3ζ1) y una cantidad de asignaciones de ζ1 o 3ζ1 a los correspondientes elementos matriciales no nulos. Similarmente laoperación de actualización completa  ←  -ζ1 D requiere sólo la sustracción de cualquiera de ζ1 or 3ζ1de las entradas correspondientes de A . Como es necesario realizar estas sustracciones solo en los elementos triangulares superiores de Â, el conteo total de operaciones para la operación de actualización se reduce a 10 (es decir 9 sustracciones más 1 operación para calcular 3ζ 1).

La matriz

contiene 12 expresiones diferentes, de las que sólo 9 han de calcularse mediante una única adición o multiplicación. (El resto se obtiene mediante asignaciones). Como ninguna de las 12 expresiones puede asumirse a priori que sean cero, el cálculo de ζ2Ê necesita 12 multiplicaciones, de manera que la cantidad total de operaciones para calcular ζ2Ê es de 21 para los elementos independientes de esta matriz. Los elementos matriciales restantes se obtienen subsecuentemente mediante asignación.

El paso final en el procedimiento de ensamble de  consiste de la operación de actualización  ←  + ζ 2

, la que requiere 21 operaciones que adicionan los triángulos superiores de las matrices  y ζ2

y las asignaciones para obtener los elementos matriciales restantes del triángulo inferior.

El procedimiento completo para ensamblar la matriz  de acuerdo con (36) realiza en secuencia los pasos anteriormente descritos, lo que requiere 24 + 42 + 10 + 21 + 21 = 118 operaciones. (¡Se debe notar que este procedimiento se configura también de una manera muy eficiente!). Como la matriz no simétrica A se calcula a partir de la matriz simétrica  de acuerdo con (35b) mediante 18 multiplicaciones, la cantidad total de operaciones para que el procedimiento anterior calcule la matriz  es de 118 + 18 =136.

En total el “procedimiento estándar” esbozado anteriormente requiere alrededor de 326 operaciones cuando se aplica a una clausura híbrida. Esto excede el costo computacional del procedimiento eficiente propuesto aproximadamente en un factor de 1.72 (es decir aproximadamente el 70%). Usar el procedimiento eficiente da el resultado de una única evaluación de la función del r.h.s. de la FTE con clausura híbrida estabilizada (incluyendo la

estabilización dinámica de la traza) en aproximadamente el 60% del costo computacional del enfoque “estándar”.

Esta mejora en la eficiencia computacional se logra principalmente mediante una reorganización inteligente de la de la función del r.h.s. basada en el cálculo analítico precedente de

, donde el nivel óptimo de eficiencia se alcanza finalmente por una reducción de las ecuaciones mediante el esquema CN.

7.7 Reescalado de la traza, monitoreo de los invariantes y proyección del espacio de fase

Cualquier matriz de FO se requiere que tenga valores propios no negativos y traza unitaria (→ sección 2.4). Esto implica un conjunto de restricciones algebraicas a las soluciones de la FTE, la que de esta manera se convierte en un DAS (→ sección 2.6). Aparte de la igualdad Sa = Tr(â) =1 para la traza, las restricciones de no negatividad para los valores propios pueden formularse equivalentemente en términos de un par de desigualdades Ka ≥ 0 y Da ≥ 0 (→

ecs. (9), (10)) para los invariantes 2do y 3ro y Da = det(â) de la matriz.

La condición de la traza se tiene en cuenta en una manera “activa” y económica computacionalmente mediante la estabilización dinámica de la traza (DTS) por un término de penalización adecuado adicionado al r.h.s. de la FTE (→ sección 7.3) que automáticamente mantiene la traza de la solución numérica cerca del valor requerido Sa =1. Usando esta técnica se obtiene una solución numérica de la FTE que cumple la condición de la traza sólo aproximadamente (es decir Sa ≈ 1). Desde un punto de vista práctico esto difícilmente afecta la corrección de la solución, dado que siempre es posible reescalar los elementos matriciales mediante una multiplicación con un factor de 1/Sa. Esta operación de reescalado de la traza se describe matemáticamente mediante

y tiene varias propiedades favorables:

Si µk son los valores propios

, los valores propios de la matriz reescalada se dan por µ'k=µk/Sa, so

1.

(i) de manera que no cambian su signo si Sa > 0 , y

2.

(ii) sus valores numéricos son sólo ligeramente reescalados si Sa ≈ 1 (como se garantiza por la

Los correspondientes vectores propios no se afectan por la operación de reescalado.

Por lo tanto la operación de reescalado de la traza no induce cambios cualitativos a la solución numérica calculada, sino corrige ligeramente la solución mientras preserva todas las propiedades esenciales de la matriz Desde este punto de vista, cualquier solución numérica de la FTE que tiene valores propios no negativos µk (o equivalentemente: invariantes 2do y 3ro Ka y Da no negativos)y Sa ≈ 1 puede interpretarse en la práctica como una matriz de FO.

Aspectos generales del control de invariantes para la FTE

De acuerdo con los resultados de [12] (→ sección 3), la no negatividad de los valores propios de una matriz 3 × 3

simetría real se garantiza si su traza es positiva y los invariantes 2do y 3ro Ka y Da son no negativos.

Como la traza de una matriz es una función lineal simple de los elementos de la matriz y el operador derivada

material

conmuta con la operación de traza, es posible derivar una ecuación de evolución para la traza de una manera directa (→ sección 4). Debido a la estructura algebraica especial del r.h.s. de la FTE esta ecuación de evolución toma una forma cerrada, en el caso “exacto” así como en los casos de las aproximaciones de clausura ortotrópica híbrida o general (ver las ecs. (14) y (15a/b)). Debido a estas propiedades algebraicas especiales de la FTE es posible adicionar un término de penalización al r.h.s. de la FTE misma sin afectar la ecuación de evolución para la traza de manera especial (intencionada) (→ DTS).

Sin embargo, a diferencia de la traza los invariantes 2do y 3ro Ka y Da son funciones no lineales de los elementos matriciales. Por lo tanto es imposible aplicar el mismo procedimiento (o uno similar) con el objetivo de obtener las ecuaciones de evolución en forma cerrada para estos invariantes (en general el r.h.s. de las ecuaciones de evolución de los invariantes es una función de todos los elementos matriciales, no sólo las combinaciones invariantes de ellos, y no es independiente de las propiedades direccionales de la matriz que se codifican en los vectores propios). Similarmente es más bien difícil (sino imposible) adicionar cualesquiera términos específicos al

r.h.s. de la FTE que afecten la evolución de los invariantes no lineales de una manera predeterminada. Por lo tanto no se dispone de un control “activo” de los invariantes 2do y 3ro por medio de unos términos de penalización adecuados en el r.h.s. de la FTE para mantener los invariantes dentro del dominio no negativo.

Un procedimiento alternativo consiste en una combinación de un procedimiento de integración “no controlado” (o

parcialmente controlado) con una operación de proyección que mapea la solución numérica de regreso a su dominio admisible. Tal procedimiento “pasivo” resulta en un esquema de integración válido que preserva la precisión del esquema “no controlado” si las operaciones de integración numérica y la proyección subsecuente se aplican en cada

paso de tiempo (ver el cap. IV.4 de [9] sobre los métodos de proyección para una discusión detallada).

Aplicado a la FTE este procedimiento combinado consiste de

1.

un paso de integración usando los métodos descritos en las secciones 5, 6 y 7 (hasta la subsección precedente 7.6) seguido por

2.

un paso de monitoreo de invariantes para determinar si la matriz calculada en el paso de integración precedente es aún una matriz de FO, y

3.

un paso de proyección en el espacio de fase el que mapea la matriz de regreso al conjunto de fase MFT si el resultado del paso de monitoreo ha sido negativo (es decir las condiciones de los invariantes se han violado).

Como el espacio de fase MFT es un conjunto convexo (→ Teorema 1 en la sección 3), para cada matriz 3 × 3 simétrica, real existe una única matriz en MFT a una “distancia mínima” (medida con una métrica adecuada, ver más adelante). Esto define el mapeo de proyección requerido para el tercer paso del procedimiento anterior. Detalles adicionales de la proyección del espacio de fase se discuten más adelante.

Monitoreo eficiente de los invariantes y cálculo de los valores propios

De acuerdo con el procedimiento descrito anteriormente, un monitoreo de los invariantes ha de realizarse en cada paso de tiempo y en cada celda (llena) del dominio de cálculo. Por lo tanto el cálculo de los invariantes – ya sea el invariante 2do o 3ro o los valores propios, los que ellos mismos son magnitudes invariantes también – ha de llevarse a cabo de la manera más eficiente para evitar cualquier sobrecarga computacional innecesaria.

Usando la notación (x, y, z) y(u, v, w) para los elementos de la diagonal y fuera de la diagonal de la matriz â como se introdujo en la sección 3, los invariantes se dan como expresiones polinómicas de los elementos matriciales mediante las fórmulas

las que en conjunto pueden evaluarse con 11 adiciones y 13 multiplicaciones. Como la traza Sa se calcula en cada paso de tiempo durante la integración numérica de la FTE, un monitoreo de los invariantes mediante el chequeo de las desigualdades Ka ≥ 0 y Da ≥ 0 puede hacerse a un costo de 22 operaciones adicionales por paso de tiempo para cada celda..

Una verificación de los invariantes Ka y Da es la manera más eficiente para determinar si una matriz 3×3 simétrica, real pertenece al conjunto MFT de matrices de FO, dado que una verificación de las condiciones µk ≥ 0 de los valores propios requiere un cálculo explícito de estos últimos.

Un análisis detallado muestra que un cálculo numérico de los valores propios de la matriz a puede hacerse más eficientemente mediante el cálculo de las raíces de su polinomio característico Pa (µ), lo que cuesta aproximadamente 100 operaciones sobre las 24 operaciones necesarias para calcular los invariantes como coeficientes de Pa (µ) y por lo tanto es aproximadamente 5 veces más caro. Esto cuantifica la sobrecarga computacional de usar los valores propios en lugar de los invariantes para monitorear las propiedades de la matriz de FO y enfatiza que usar los invariantes (38) es el método a elegir para un control eficiente de las propiedades espectrales de a . Una diagonalización numérica de la matriz es el método más caro para calcular los valores propios. Por lo tanto tal procedimiento no se recomienda a menos que los vectores propios también sean de interés (por ejemplo como parte de un procedimiento de proyección de espacio de fase discutido a continuación).

Proyección del espacio de fase

Si el monitoreo de los invariantes indica que la solución numérica de la FTE se mueve fuera del conjunto de espacio de fase MFT (es decir si una o ambas de las desigualdades Ka ≥ 0 y Da ≥ 0 se violan) dentro del paso de tiempo actual, la solución ha de corregirse por medio de una proyección sobre el espacio de fase MFT.

Para cualquier matriz 3×3 real

F de dos matrices 3 × 3 reales se mide en la norma de Frobenius
צ

-

צ ) =,(donde la distancia dF ) en el espacio vectorial de las ·T=Tr(
ۄ |
ޛ la que es autoinducida por el producto escalar

matrices 3 × 3 reales. La restricción de estas magnitudes al subespacio de seis dimensiones de las matrices

א â*existe a una única matriz

simétricas es directa. Puede mostrarse [18] que para cualquier matriz 3 × 3 real MFT tal que la distancia dF( , â) entre

y los elementos â de MFT se minimiza precisamente por â = â*.

que se define únicamente mediante las siguientes FTM
א â*Esta solución de distancia mínima se da por una matriz

propiedades (ver [18] para una prueba formal):

Los vectores propios de â* son los mismos que los de la matriz .

La terna de valores propios de â* es el único punto del conjunto del triángulo (ver la ec. (7) en la sección 3 y la Fig. 7)

,FTM
א âste mapeo de proyección da por definición la matriz “más cercana” e

que tiene la menor distancia euclidiana en R3 hasta el punto (µ1, µ2, µ3)dada por la terna de valores propios de la matriz .

El mapeo de proyección definido en el espacio vectorial real de seis dimensiones

de las matrices 3 × 3 simétricas y que tiene valores en MFT se da en notación matemática compacta mediante la fórmula

y el algoritmo para calcular realmente el valor del mapeo de proyección consiste de los pasos siguientes:

1. Dada una matriz 3 × 3 simétrica real como entrada, calcular sus invariantes Sa, Ka y Da.

yFTM
א ≥ 0. Si todas las condiciones se cumplen, entonces ya â aD≥ 0 yaK =1, aS2. Chequear las condiciones

el mapeo de proyección es justo la identidad:

En este caso asignar la matriz de entrada a a la salida y salir. De lo contrario (es decir si ya sea que Ka < 0 o Da < 0) ejecutar los siguientes pasos:

3. Calcular los valores µk y los correspondientes vectores propios Ek o de la matriz (es decir â·Ek = µk Ek). Este paso puede representarse formalmente como:

4. A continuación encontrar el único punto que tiene distancia euclidiana mínima hasta µ = (µ1, µ2, µ3) en R3, es decir:

5. Finalmente componer el valor â* = PFT[â] del mapeo de proyección mediante

, de manera que (39d) coincide con (39a) en este caso. En este sentidoFTM
א Nótese que (39c) da µ* = µ si â

(39b/c/d) ya codifican la definición general del mapeo de proyección. Los pasos individuales para calcular los

, aunque (39a) formalmente completa la definiciónFTM
ב resultados de (39b/c/d) han de llevarse a cabo sólo para â .FTM
א âen el caso que

Aunque la descripción del algoritmo para calcular

FT[â]dado anteriormente es matemáticamente preciso, una implementación práctica de este algoritmo tomando matrices de entrada con traza Sa ≈ 1 (como se suministra por un esquema de integración en el tiempo de la FTE con DTS) omitiría una evaluación adicional de la traza en el Paso 1 del procedimiento. Aún más una implementación práctica incorporaría las siguientes variaciones y/o especificaciones de los varios pasos algorítmicos, ver la Fig. 13:

En Paso 1 los invariantes Ka y Da de la matriz de entrada se calculan mediante (38).

En el Paso 2. se chequean sólo las condiciones Ka ≥ 0 y Da ≥ 0 y la matriz de entrada se mantiene sin cambios y se asigna a la salida si ambas condiciones se cumplen.

Los resultados de Paso 3. se obtienen mediante una diagonalización matricial numérica, por ejemplo mediante rotaciones cíclicas de Jacobi iterativas (ver [38], cap. 11.1 o [39], cap. 6.5.2) o un único paso de la reducción de Givens/Householder y las iteraciones QR subsecuentes (ver [42] o [38], cap. 11.2). 11.2).

Realmente el Paso 4. se lleva a cabo sólo si al menos uno de los valores propios de entrada µk es negativo, de manera que el punto de solución

se localiza siempre en los bordes (incluyendo los vértices) del conjunto triángulo Δa.

Usar los valores propios de la “traza reescalada k = µk/Sa (ver la ec. y los comentarios subsecuentes sobre este tema) confina el problema de minimización que ha de resolverse para determinar μ* µ* de acuerdo con (39c) al plano de Sa = 1 fijado por los tres vértices (1,0,0), (0,1,0) y (0,0,1) del triángulo Δa. Este problema de minimización “plana” puede resolverse equivalentemente en R2.por ejemplo en el plano (µ1, µ2)-encontrando el único par

localizado en los bordes del triángulo proyectado

(es decir la proyección ortogonal de Δa sobre el plano (µ 1, µ2) -que tiene la menor distancia euclidiana hasta el par

( 1, 2). La tercera coordenada de µ* entonces se da por

Dado que Sa ≈ 1 mediante DTS, el algoritmo que resuelve el problema de minimización “plana” da una muy buena aproximación a la solución exacta de (39c), aunque es mucho más sencillo y puede implementarse más eficientemente que un algoritmo que calcula μ* directamente de (39c) en R3 .

En el Paso 5.la fórmula para calcular los elementos individuales de la matriz â*se dan por

en términos de los elementos de la matriz Rij de la matriz de rotación

= (E1, E2, E3) que se ha calculado previamente en el Paso 3.

8. Implementación

El resultado del proceso, el que es una función de distribución que describe la probabilidad de la orientación de las fibras en un área dada del artículo se presenta en forma gráfica o numérica en una pantalla de una computadora de una estación de trabajo (no mostrada). Esta pudiera ser la pantalla de la computadora o de la estación de trabajo en la que se realizan los cálculos, o pudiera ser en una pantalla de una computadora que está conectada en red a la computadora en la que se realiza la simulación.

Un diseñador de moldes o de productos usará los resultados de la simulación para mejorar la calidad del artículo que resulta del proceso de moldeo por inyección. Los resultados de la simulación también pueden usarse por los ingenieros para reducir el costo de fabricación del artículo en cuestión. La ventaja de una información confiable sobre la información de las fibras permite a los ingenieros tener una mejor comprensión e información antes de determinar las características de resistencia, rigidez y estabilidad del artículo, dado que la orientación de las fibras, cuyas fibras típicamente tienen características de mucho más alta resistencia que el material polímero, tienen una influencia decisiva en las características de resistencia, rigidez y estabilidad del artículo. Además, la orientación de las fibras tiene una influencia en los efectos de deformación que pueden ocurrir cuando se moldea por inyección una masa de polímero con fibras suspendidas en la misma. Conociendo la orientación de las fibras, la deformación y otros esfuerzos que crea la comisión al efecto pueden predecirse mejor o evitarse los cambios al diseño.

Los resultados de la simulación también pueden transmitirse a otras aplicaciones tales como un software CAD. Por lo tanto, los resultados de la simulación pueden usarse en un proceso interactivo entre el software de diseño y el software de simulación.

La invención tiene numerosas ventajas. Las diferentes modalidades o implementaciones pueden dar una o más de las siguientes ventajas. Debe notarse que esta no es una lista exhaustiva y pueden existir otras ventajas las que no se describen en la presente. Una ventaja de la invención es que la distribución de la orientación de las fibras de los artículos reforzados con fibras de fabricación puede determinarse a un esfuerzo computacional significativamente reducido. Una ventaja de la presente invención es que la distribución de la orientación de las fibras en los artículos reforzados con fibras de fabricación puede determinarse con una mayor precisión. Otra ventaja de la presente invención es que la distribución de la orientación de las fibras en artículos reforzados con fibras de fabricación puede determinarse con mayor rapidez.

Aunque la presente invención ha sido descrita en detalle para propósitos de ilustración, se entiende que tal detalle es exclusivamente para este propósito, y pueden hacerse variaciones de la misma por los expertos en la materia sin apartarse del alcance de la invención. El término “que comprende" como se usa en las reivindicaciones no excluye otros elementos o pasos. El término “un” o “una” como se usa en las reivindicaciones no excluye una pluralidad.

Referencias

[1] S.G. Advani (Ed.): Flow and Rheology in Polymer Composites Manufacturing, Elsevier, Amsterdam (1994)

[2] G.B. Jeffery: The motion of ellipsoidal particles immersed in a viscous fluid, Proc. R. Soc. A 102, 161-179 (1922)

[3] M. Junk y R. Illner: A new derivation of Jeffery's equation, J. Math. Fluid Mech. 8, 1-34 (2006)

[4] F. Folgar y C.L. Tucker III: Orientation behaviour of fibers in concentrated suspensions, J. Reinf. Plast. Compos. 3, 98-119 (1984)

[5] C.L. Tucker y S.G. Advani: Processing of short-fiber systems, ch. 6, pp. 147 en S.G. Advani (Ed.): Flow and Rheology in Polymer Composites Manufacturing (ver referencia [1] anterior)

[6] S.G. Advani y C.L. Tucker III: The use of tensors to describe and predict fiber orientation in short fiber composites, J. Rheol. , 751-784 (1987)

[7] S.G. Advani y C.L. Tucker III: Closure approximations for three-dimensional structure tensors, J. Rheol., 367386 (1990)

[8] J.S. Cintra y C.L. Tucker III: Orthotropic closure approximations for flow-induced fiber orientation, J. Rheol. 39, 1095-1122 (1995)

[9] E. Hairer, C. Lubich y G. Wanner: Geometric numerical integration, Springer, Berlin (2002)

[10] J. Linn, J. Steinbach, A. Reinhardt: Calculation of the 3D fiber orientation in the simulation of the injection molding process for short-fiber reinforced thermoplasts, Conferencia de ECMI 2000, Palermo (2000)

[11] J. Linn: Exploring the phase space of the Folgar-Tucker equation, Conferencia de SIAM-EMS, Berlin (2001)

[12] J. Linn: On the frame-invariant description of the phase space of the Folgar-Tucker equation, p. 27-332 in A.

Buikis, R. iegis, A.D. Fitt (Eds.): Progress in Industrial Mathematics at ECMI 2002, Springer (2004)

[13] S. Hess: Fokker-Planck equation approach to flow alignment in liquid crystals, Z. Naturforsch. 31A, 1034 ff. (1976)

[14] M. Doi: Molecular dynamics and rheological properties of concentrated solutions of rodlike polymers in isotropic and liquid cristalline phases, J. Polym. Sci., Polym. Phys. Ed. 19, 229-243 (1981)

[15] M. Grosso, P.L. Maffetone, F. Dupret: A closure approximation for nematic liquid crystals based on the canonical distribution subspace theory, Rheol. Acta 39, 301-310 (2000)

[16] M. Kröger: Simple models for complex nonequilibrium fluids, Phys. Rep. 390, 453-551 (2004)

[17] J. Linn: The Folgar-Tucker Model as a Differential Algebraic System for Fiber Orientation Calculation, pp. 215-224 en Y. Wang, K. Hutter: Trends in Applications of Mathematics to Mechanics, ponencias de la STAMM 2004 en Seeheim (Alemania), Shaker (2005)

[18] U. Strautins: Investigation of fiber orientation dynamics within the Folgar-Tucker model with hybrid closure, Tesis de Maestría, Departamento de Matemáticas, Universidad de Kaiserslautern (2004)

[19] J. Linn: Dreidimensionale Vorausberechnung der anisotropen Materialeigenschaften in thermoplastischen Spritzgusserzeugnissen (Anlso-3D), Projektbericht für die MAGMA GmbH, Teil IIa (2002)

[20] J. Linn: Dreidimensionale Vorausberechnung der anisotropen Materialeigenschaften in thermoplastischen Spritzgusserzeugnissen (AnIso-3D), Projektbericht für die MAGMA GmbH, Teil IIb (2002)

[21] B.E. VerWeyst, C.L. Tucker, P.H. Foss, J.F. O'Gara: Fiber orientation in 3D injection moulded features: prediction and experiment, Internat. Polymer Processing 14, 409-420 (1999);

[22] B.E. VerWeyst: Numerical predictions of flow-induced fiber orientation in three-dimensional geometries, Tesis Ph.D, Univ. de Illinois en Campaña Urbana (1998)

[23] G.I. Marchuk: Splitting and Alternating Direction Methods, pp. 197-462 en P.G. Ciaret & J.L Lions (Eds.): Handbook of Numerical Analysis, Volumen I, North-Holland, Elsevier (1990)

[24] K.W. Morton: Numerical solution of convection-diffusion problems, Chapman & Hall, London (1996)

[25] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser (1992)

[26] G. Strang: "On the construction and comparison of difference schemes", Revista SIAM Num. Anal. 5, 506517 (1968)

[27] M.G. Crandall y A. Majda: "The method of fractional steps for conservation laws", Math. Comp. 34, 285-314 (1980)

[28] H.V. Kojouharov, B.M. Chen: "Nonstandard methods for the convective transport equation with nonlinear reactions", Numer. Meth. Partial Diff. Eq. 14, 467-485 (1998); "Nonstandard methods for the convectivedispersive transport equation with nonlinear reactions" en R.E. Mickens (ed.): Applications of non-standard finite difference schemes, minisymposium on non-standard finite difference schemes: theory and applications, Reunión Anual de SIAM , Atlanta GA, USA 1999, publicado por Singapore: World Scientific (2000)

[29] H. Wang, X. Shi y R.E. Ewing: "An ELLEM scheme for multidimensional advection-reaction equations and its optimal-order error estimate", Revista SIAM Numer. Anal. 38, 1846-1885 (2001)

[30] P.J. van der Houwen: "Note on the time integration of 3D advection-reaction equations", J. Comput. Appl. Math. 116, 275-278 (2000)

[31] W. Hunsdorfer, J.G. Verwer: "A note on splitting errors for advection-reaction equations", Appl. Numer. Math. 18, 191-199 (1995)

[32] S.V. Patankar: Numerical heat transfer and fluid flow, Hemisphere Publ. Corp. (1980)

[33] C.A.J. Fletcher: Computational techniques for Fluid Dynamics, Volumen I: Fundamental and General Techniques (2da edición), Springer (1991)

[34] L.F. Shampine: "Conservation laws and the numerical solution of ODEs", Comp. Math. Applic. 12B, 12871296 (1986)

[35] L.F. Shampine: "Linear conservation laws for ODEs", Comp. Math. Applic. 35, 45-53 (1998)

[36] L.F. Shampine: "Conservation laws and the numerical solution of ODEs, parte II", Comp. Math. Applic. 38, 5 61-72 (1999)

[37] J. Linn: "Entwicklung eines Software-Moduls zur Berechnung der Faserorientierung in der Spritzgießsimulation mit SIGMASOFT", Reporte Técnico para GmbH MAGMA (2001)

[38] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery: Numerical Recipes in Fortran 77: The Art of

Scientific Computing (2da Edición), Universidad de Cambridge (1992) 10 [39] J. Stoer, R. Bulirsch: Introduction to Numerical Analysis (3ra Edición), Springer (2002)

[40] D.H. Chung, T.H. Kwon: "An invariant-based optimal fitting closure approximation for the numerical prediction of flow-induced fibre orientation", J. Rheol. 46, 169-194 (2002)

[41] F. Dupret, V. Verleye, B. Languilier: "Numerical prediction of moulding of short-fibre composite parts", Conf. Proc. 1er ESAFORM, 291-294 (1998)

15 [42] G.H. Golub, H.A. van der Vorst: "Eigenvalue Computation in the 20th Century", J. Comp. Appl. Math. 123, 35-65 (2000)

Claims (16)

1. REIVINDICACIONES

   1.

   Un método para calcular las estadísticas de la orientación de partículas no esféricas a nivel macroscópico mediante el uso de un modelo de simulación para simular un proceso de moldeo por inyección en el que una cavidad de molde, la que forma al menos parte de la geometría de un dominio de simulación, se llena o parcialmente se llena con una suspensión que se forma por un solvente que contiene un gran número de partículas no esféricas, en donde se proporciona una representación digital o un modelo de computadora de la geometría del dominio de simulación, y en donde se forma una malla con una pluralidad de celdas de cálculo mediante la subdivisión o discretización al menos de una parte del dominio de simulación, el método que comprende:

   a) especificar las condiciones de frontera; b) establecer las condiciones iniciales; c) resolver las ecuaciones de balance de masa, el momento y la energía de al menos una porción de las celdas del dominio de simulación para obtener un flujo de fluido, un flujo de calor y una transferencia de masa a nivel macroscópico; y d) resolver las ecuaciones dinámicas de orientación de partículas no esféricas basado al menos parcialmente en los resultados de las ecuaciones de balance resueltas para de esta manera determinar los cambios en la orientación de las partículas no esféricas a nivel macroscópico como función del espacio y del tiempo, en donde la orientación de las partículas no esféricas se describe estadísticamente por una función de distribución y en donde la función de distribución se aproxima por tensores de orientación de partículas no esféricas de orden creciente y el tensor de orientación de 4to orden se calcula como una función del tensor de 2do orden usando una aproximación de clausura híbrida estabilizada.

2. 2.

   Un método de acuerdo con la reivindicación 1, en donde la función de distribución es un sistema acoplado no lineal de ecuaciones diferenciales parciales hiperbólicas de reacción convección (16) que se reformula en forma vectorial usando notación contraída (CN).

3. 3.

   Un método de acuerdo con cualquiera de las reivindicaciones 2, en donde la estructura algebraica del lado derecho de un sistema de ecuaciones diferenciales parciales hiperbólicas de reacción convección (16) se explota para lograr una clausura híbrida estabilizada con pocas operaciones algebraicas.

4. 4.

   Un método de acuerdo con la reivindicación 3, en donde el número de operaciones algebraicas requeridas para determinar la orientación de las partículas no esféricas se minimiza mediante la identificación de subexpresiones comunes en el lado derecho del sistema de ecuaciones diferenciales parciales hiperbólicas (16).

5. 5.

   Un método de acuerdo con cualquiera de las reivindicaciones 3 ó 4, en donde la ecuación de Folgar-Tucker (16) se reformula mediante la sustitución de fórmulas de definición para la clausura híbrida en la función del r.h.s. para cálculos analíticos y la reorganización subsecuente de los términos resultantes para lograr una estructura de términos la que es óptima para una evaluación eficiente..

6. 6.

   Un método de acuerdo con cualquiera de las reivindicaciones 3 hasta la 5, en donde un término de control se incluye en la ecuación de Folgar-Tucker reformulada (16) para la estabilización dinámica de la traza.

7. 7.

   Un método de acuerdo con cualquiera de las reivindicaciones 3 hasta la 6, en donde la función del r.h.s. reformulada se evalúa con un número mínimo de operaciones usando CN para una eficiencia óptima.

8. 8.

   Un método de acuerdo con cualquiera de las reivindicaciones 3 hasta la 7 que además incluye chequear si la solución numérica de la FTE está en el dominio admisible mediante el control de los invariantes 2do y 3ro.

9. 9.

   Un método de acuerdo con cualquiera de las reivindicaciones 3 hasta la 8, que además incluye el uso de un

   “reescalado de traza” para interpretar matrices de FO que tienen por medio de DTS sólo aproximadamente traza

   unitaria.

10. 10.

    Un método de acuerdo con cualquiera de las reivindicaciones 3 hasta la 9, que incluye una descripción algorítmica de la proyección del espacio de fase.

11. 11.

    Un método de acuerdo con la reivindicación 10, que además incluye una combinación de la proyección del espacio de fase con las otras técnicas propuestas para una integración numérica de la FTE para confinar la solución numérica al dominio admisible.

12. 12.

    Un método de acuerdo con la reivindicación 1, en donde se usa una solución aproximada para el tensor de orientación de 4to orden, dicho tensor de orientación de 4to orden aproximado que se calcula a partir del tensor de orientación de 2do orden mediante el uso de una aproximación de clausura.

13. 13.

    Un método de acuerdo con la reivindicación 12, en donde la aproximación de clausura es una aproximación de clausura híbrida.

14. 14.

    Un método de acuerdo con la reivindicación 13, en donde la aproximación de clausura híbrida se ha estabilizado

    5 para lograr la estabilidad numérica y en donde la estabilización de la aproximación de clausura híbrida incluye el uso de un factor de orientación escalar fSC, el que se confina a un rango de valores numéricos de [0 1].
15. 15. Un método de acuerdo con la reivindicación 14, en donde el valor del factor de orientación escalar se confina al

    rango de valores numéricos [0 1], estableciendo el valor igual a 0 cuando el valor calculado de fSC es menor que 0, y 10 estableciendo el valor igual a 1 cuando el valor calculado de fSC es superior a 1.
16. 16. Un producto de software de computadora en un medio legible por computadora que incluye un código de software para ejecutar un método de acuerdo con cualquiera de las reivindicaciones 1 hasta la 15.