CoMPhy Lab blogs

Non-Newtonian flows with Basilisk

20 min read

PDF version

Download this page as PDF

TL;DR

This assignment explores computational approaches to viscoelastic fluid dynamics using Basilisk. We investigate four advanced multiphase flow phenomena involving complex rheology: bursting bubbles in polymer solutions, viscoplastic drop impact, soft material impacts, and Taylor-Culick retraction of viscoelastic filaments. Each simulation incorporates non-Newtonian behavior characterized by dimensionless parameters including Weissenberg, Deborah, Ohnesorge, and elasto-capillary numbers. Through these simulations, we develop practical skills in implementing log-conformation formulations for numerical stability, handling complex interfaces with adaptive mesh refinement, and visualizing rheologically complex flows. This assignment represents the culmination of our computational fluid dynamics progression, building on earlier work with heat conduction, single-phase flows, and simple multiphase systems.

Introduction

In this final assignment, we explore the fascinating world of viscoelastic fluid dynamics using computational simulations. Viscoelastic materials—from polymer solutions to biological fluids, food products, and advanced materials—exhibit both viscous (fluid-like) and elastic (solid-like) properties, leading to complex and often counterintuitive behaviors. Understanding these behaviors is crucial for applications ranging from inkjet printing and food processing to biological systems and advanced manufacturing.

Building on our previous work with heat conduction (1st-workingAssignment), single-phase flows (2nd-workingAssignment), and simple multiphase flows (3rd-workingAssignment), we now tackle the additional complexity introduced by non-Newtonian rheology, viscoelasticity, and viscoplasticity. These phenomena represent some of the most challenging aspects of computational fluid dynamics, requiring specialized numerical techniques for stable and accurate simulations.

Learning Objectives

  • Understand key mathematical models for viscoelastic and viscoplastic fluid behavior
  • Implement log-conformation formulations for numerical stability in viscoelastic simulations
  • Analyze how elasticity affects interfacial dynamics in multiphase flows
  • Apply adaptive mesh refinement strategies to resolve complex interface evolution
  • Develop skills in visualizing and interpreting non-Newtonian flow phenomena
  • Relate dimensionless parameters to physical behaviors in complex fluids

Throughout this assignment, we’ll explore four complex multiphase flow simulations:

  1. Bursting bubbles in viscoelastic media
  2. Viscoplastic drop impact on a solid surface
  3. Soft impact dynamics with viscoelastic fluids
  4. Taylor-Culick retraction of viscoelastic filaments

Each case presents unique physical phenomena and numerical challenges, providing a comprehensive overview of computational approaches to viscoelastic fluid dynamics.

Part 1: Viscoelastic Fluid Fundamentals

Before diving into specific simulations, let’s establish the mathematical framework for understanding viscoelastic fluid behavior.

1.1 Constitutive Equations for Complex Fluids

Traditional Newtonian fluids are characterized by a linear relationship between stress and strain rate, governed by a constant viscosity:

where is the stress tensor, is the viscosity, and is the strain rate tensor.

For non-Newtonian fluids, this relationship becomes more complex:

  1. Viscoplastic fluids (e.g., Bingham plastics, Herschel-Bulkley fluids) exhibit a yield stress that must be exceeded before flow occurs:

  1. Viscoelastic fluids exhibit memory effects, with stress depending on the deformation history. A common model is the Oldroyd-B model:

where is the relaxation time, is the retardation time, and denotes the upper-convected derivative of the stress tensor.

1.2 Key Dimensionless Parameters

Several dimensionless parameters characterize the behavior of complex fluids:

  1. Weissenberg number (Wi): Ratio of elastic to viscous forces where is the relaxation time and is a characteristic strain rate

  2. Deborah number (De): Ratio of material relaxation time to observation time

  3. Elasticity number (El): Ratio of elastic to inertial effects

  4. Ohnesorge number (Oh): Ratio of viscous forces to inertial and surface tension forces

  5. Plasto-capillary number (J): Ratio of yield stress to capillary pressure

Numerical Challenge Simulating viscoelastic flows presents significant numerical challenges, particularly at high Weissenberg numbers (known as the "high Weissenberg number problem"). Traditional numerical methods often fail due to the loss of positive-definiteness of the conformation tensor. The log-conformation approach we'll use in these simulations addresses this issue by working with the logarithm of the conformation tensor, ensuring numerical stability.

1.3 The Log-Conformation Approach

In our simulations, we employ the log-conformation approach, which transforms the conformation tensor logarithmically:

The evolution equation is then solved for rather than directly for , which ensures that the reconstructed conformation tensor remains positive-definite. This approach significantly improves numerical stability, especially for high Weissenberg number flows.

The total stress in the fluid is the sum of the Newtonian solvent contribution and the polymer stress:

where is the solvent viscosity, is the elastic modulus, and is the identity tensor.

Part 2: Bursting Bubbles in Viscoelastic Media

Our first simulation explores the dynamics of bubble bursting in viscoelastic fluids, a phenomenon that differs significantly from the Newtonian case due to elastic memory effects.

2.1 Physical Description

When a bubble reaches a free surface, the thin film separating it from the surrounding atmosphere eventually ruptures. This triggers a series of events:

  1. Capillary waves travel along the interface
  2. The bubble cavity collapses
  3. A central jet may form and rise upward (Worthington jet)
  4. Depending on conditions, the jet may break into droplets

In viscoelastic fluids, elastic stresses dramatically alter these dynamics, often leading to enhanced jet formation and reduced droplet breakup due to elastic resistance to extensional deformation.

Problem Statement Simulate the bursting of a bubble at the interface between a viscoelastic fluid and air, examining how elasticity affects the jet dynamics and potential droplet formation. Key dimensionless parameters include:

  • Deborah number (De): Ratio of relaxation time to flow time scale
  • Elasto-capillary number (Ec): Ratio of elastic forces to surface tension
  • Ohnesorge number (Oh): Ratio of viscous forces to inertial and capillary forces
  • Bond number (Bond): Ratio of gravitational forces to surface tension

2.2 Implementation Details

The simulation uses Basilisk’s axisymmetric framework with a log-conformation approach for viscoelasticity:

#include "axi.h"                  // Axisymmetric coordinates
#include "navier-stokes/centered.h"  // Main Navier-Stokes solver
 
// Viscoelastic model implementation
#if !_SCALAR
#include "log-conform-viscoelastic.h"  // Tensor implementation
#else 
#include "log-conform-viscoelastic-scalar-2D.h"  // Scalar implementation
#endif
 
#define FILTERED                           // Smear density and viscosity jumps
#include "two-phaseVE.h"                   // Two-phase flow with viscoelasticity
#include "navier-stokes/conserving.h"      // Conservative momentum advection
#include "tension.h"                       // Surface tension forces

The bubble is initialized using a pre-calculated equilibrium shape that balances surface tension and gravitational forces:

event init (t = 0) {
  if (!restore (file = dumpFile)){
    // Read initial shape from data file
    char filename[60];
    sprintf(filename,"Bo%5.4f.dat",Bond);
    FILE * fp = fopen(filename,"rb");
    
    // Generate distance field from shape coordinates
    coord* InitialShape;
    InitialShape = input_xy(fp);
    fclose (fp);
    
    // Compute signed distance function
    scalar d[];
    distance (d, InitialShape);
    
    // Refine mesh based on distance function
    while (adapt_wavelet ((scalar *){f, d}, (double[]){1e-8, 1e-8}, MAXlevel).nf);
    
    // Convert distance function to volume fraction
    vertex scalar phi[];
    foreach_vertex(){
      phi[] = -(d[] + d[-1] + d[0,-1] + d[-1,-1])/4.;
    }
    fractions (phi, f);
  }
}

Physical properties are specified based on the dimensionless parameters:

// Density ratio (liquid : gas = 1000:1)
rho1 = 1., rho2 = 1e-3;
  
// Gas phase Ohnesorge number
Oha = 2e-2 * Oh;
  
// Viscosity from Ohnesorge numbers
mu1 = Oh, mu2 = Oha;
  
// Polymer relaxation time (liquid has elasticity, gas does not)
lambda1 = De; lambda2 = 0.;
  
// Elastic modulus (liquid is viscoelastic, gas is not)
G1 = Ec; G2 = 0.;
 
// Surface tension coefficient
f.sigma = 1.0;

Adaptive Mesh Refinement The simulation uses multiple criteria for mesh refinement to accurately resolve the interface dynamics: 

event adapt(i++){
  // Calculate interface curvature for refinement criterion
  scalar KAPPA[];
  curvature(f, KAPPA);
  
  // Adapt mesh based on multiple criteria
  adapt_wavelet ((scalar *){f, u.x, u.y, conform_p.x.x, conform_p.y.y, 
     conform_p.y.x, conform_qq, KAPPA},
     (double[]){fErr, VelErr, VelErr, AErr, AErr, AErr, AErr, KErr},
     MAXlevel, MAXlevel-6);
}

This adaptive approach concentrates computational resources where they’re needed most: at the interface, in regions of high curvature, areas with significant velocity gradients, and regions with rapidly varying conformation tensor components.

2.3 Expected Results

For Newtonian fluids (De = 0), bubble bursting typically results in:

  1. Rapid retraction of the film
  2. Formation of capillary waves
  3. Collapse of the cavity
  4. Potential formation of a Worthington jet
  5. Possible droplet pinch-off from the jet

For viscoelastic fluids (De > 0), we expect to observe:

  1. Enhanced jet formation due to elastic stresses
  2. Delayed or suppressed droplet pinch-off
  3. Longer and thinner jets at high Deborah numbers
  4. Beads-on-a-string morphology for moderate elasticity

Exercise: Investigating Deborah Number Effects

Modify the bubble bursting simulation to explore different Deborah numbers:

  1. Run simulations with De = 0 (Newtonian), De = 1 (moderate elasticity), and De = 10 (high elasticity), keeping other parameters constant
  2. Track and compare:
    • Maximum jet height
    • Time to first droplet pinch-off (if it occurs)
    • Jet morphology (thickness, length, shape)
  3. Explain the physical mechanisms behind the observed differences

This exercise will help you understand how elastic memory affects interfacial dynamics.

Part 3: Viscoplastic Drop Impact

Our second simulation investigates the dynamics of a viscoplastic (yield-stress fluid) drop impacting a solid surface.

3.1 Physical Background

When a viscoplastic drop impacts a surface, its behavior differs significantly from Newtonian drops due to the yield stress, which creates unyielded (solid-like) regions within the fluid. Key features include:

  1. Initial deformation driven by inertia
  2. Formation of unyielded regions where stress falls below the yield stress
  3. Limited spreading compared to Newtonian drops
  4. Potential formation of a solid-like cap at the drop apex
  5. Asymmetric retractions and frozen interfaces after impact

Problem Statement Simulate the axisymmetric impact of a viscoplastic drop onto a solid surface. Key dimensionless parameters include:

  • Weber number (We): Ratio of inertial to surface tension forces
  • Ohnesorge number (Oh): Ratio of viscous to inertial and surface tension forces
  • Plasto-capillary number (J): Ratio of yield stress to capillary pressure

3.2 Mathematical Formulation

For the viscoplastic fluid, we employ the Herschel-Bulkley model with regularization:

where:

  • is the deformation rate tensor
  • is its second invariant
  • is a regularization parameter
  • is the power law index (n = 1 for Bingham model)
  • is the dimensionless yield stress
  • is the Ohnesorge number

The regularization ensures numerical stability by avoiding the mathematical singularity at zero deformation rate.

3.3 Implementation Details

The simulation leverages a modified two-phase flow solver for viscoplastic fluids:

#include "axi.h"
#include "navier-stokes/centered.h"
#define FILTERED // Smear density and viscosity jumps
#include "two-phaseVP-HB.h"
#include "navier-stokes/conserving.h"
#include "tension.h"
#include "reduced.h"

Boundary conditions specify a solid wall at the left boundary:

// Left boundary: solid wall (no-slip and no flux)
u.t[left] = dirichlet(0.0);  // Tangential velocity = 0 (no-slip)
f[left] = dirichlet(0.0);    // Volume fraction = 0 (solid wall)

The drop is initialized with a downward velocity and positioned slightly above the wall:

event init (t = 0) {
  if (!restore (file = dumpFile)){
    refine((R2Drop(x, y) < 1.05) && (level < MAXlevel));
    fraction(f, 1. - R2Drop(x, y));
    foreach() {
      u.x[] = -1.0 * f[];  // Initial velocity in x-direction (toward wall)
      u.y[] = 0.0;         // No initial velocity in y-direction
    }
  }
}

Adaptive Refinement Strategy The simulation employs a sophisticated adaptive refinement strategy that changes based on the simulation stage: 

event adapt(i++){
  if (t < 1e-2){
    // Early impact phase: focus on interface and velocity fields
    adapt_wavelet ((scalar *){f, u.x, u.y},
    (double[]){fErr, VelErr, VelErr},
    MAXlevel);
  } else {
    // Later stages: add refinement based on curvature and deformation
    scalar KAPPA[], D2c[];
    curvature(f, KAPPA);
    foreach() {
      // Calculate deformation tensor components
      double D11 = (u.y[0,1] - u.y[0,-1])/(2*Delta);
      double D22 = (u.y[]/y);
      double D33 = (u.x[1,0] - u.x[-1,0])/(2*Delta);
      double D13 = 0.5*((u.y[1,0] - u.y[-1,0] + u.x[0,1] - u.x[0,-1])/(2*Delta));
      double D2 = (sq(D11)+sq(D22)+sq(D33)+2.0*sq(D13));
      D2c[] = f[]*(D2);  // Deformation rate in the drop phase
    }
    adapt_wavelet ((scalar *){f, u.x, u.y, KAPPA, D2c},
    (double[]){fErr, VelErr, VelErr, KAPPAErr, D2Err},
    MAXlevel);
  }
}

This strategy ensures fine resolution near the interface and in regions of high deformation, which is critical for accurately capturing the yield surfaces where the stress transitions across the yield threshold.

3.4 Expected Results

For viscoplastic drops, the impact dynamics depend strongly on the plasto-capillary number J:

  1. Low J (J << 1): Behavior approaches Newtonian drops with nearly complete spreading
  2. Moderate J (J ~ 0.1-1): Limited spreading with formation of unyielded regions
  3. High J (J > 1): Minimal deformation, approaching solid-like behavior

A key feature to observe is the formation of unyielded regions where the stress falls below the yield stress. These regions behave as solid-like domains embedded within the flowing fluid.

Exercise: Non-Newtonian Impact Dynamics

Compare viscoplastic drop impact with Newtonian drop impact:

  1. Run simulations with J = 0 (Newtonian), J = 0.1 (moderate yield stress), and J = 0.5 (high yield stress)
  2. Compare:
    • Maximum spread factor (ratio of maximum diameter to initial diameter)
    • Final equilibrium shape
    • Locations of unyielded regions (where applicable)
  3. Plot the spread factor vs. time for each case on a single graph
  4. Explain how yield stress affects the impact dynamics

This exercise will help you understand how yield stress modifies impact behavior compared to Newtonian fluids.

Part 4: Soft Impact Dynamics with Viscoelastic Fluids

Our third simulation explores the dynamics of viscoelastic drop impact, focusing on the effects of elasticity on the impact behavior.

4.1 Physical Background

When a viscoelastic drop impacts a surface, elastic stresses can lead to behaviors not seen in Newtonian impacts:

  1. Enhanced rebound due to elastic energy storage and recovery
  2. Suppression of satellite droplet formation
  3. Extended filament formation during rebound
  4. Complex fingering patterns during spreading
  5. Beads-on-a-string morphology in retracting filaments

Problem Statement Simulate the impact of a viscoelastic drop on a solid surface, examining how elasticity affects the impact dynamics. Key dimensionless parameters include:

  • Weber number (We): Ratio of inertial to surface tension forces
  • Ohnesorge number (Oh): Ratio of viscous to inertial and surface tension forces
  • Weissenberg number (Wi): Ratio of elastic to viscous forces
  • Elasticity number (El): Ratio of elastic to inertial effects
  • Bond number (Bo): Ratio of gravitational to surface tension forces

4.2 Implementation Details

The simulation uses a similar framework to the bubble bursting case, with axisymmetric formulation and log-conformation approach:

#include "axi.h"                  // Axisymmetric coordinates
#include "navier-stokes/centered.h"  // Main Navier-Stokes solver
 
// Viscoelastic model implementation
#if !_SCALAR
#include "log-conform-viscoelastic.h"  // Tensor implementation
#else 
#include "log-conform-viscoelastic-scalar-2D.h"  // Scalar implementation
#endif
 
#define FILTERED                           // Smear density and viscosity jumps
#include "two-phaseVE.h"                   // Two-phase flow with viscoelasticity
#include "navier-stokes/conserving.h"      // Conservative momentum advection
#include "tension.h"                       // Surface tension implementation

The drop is initialized as a sphere approaching the axis of symmetry (which represents the impact surface in the axisymmetric formulation):

event init (t = 0) {
  if (!restore (file = dumpFile)){
    // Refine mesh around interface
    refine((R2Drop(x, y) < 1.05) && (level < MAXlevel));
    
    // Initialize volume fraction
    fraction(f, 1. - R2Drop(x, y));
    
    // Set initial velocity field
    foreach() {
      u.x[] = -1.0 * f[];  // Droplet moves left with velocity -1.0
      u.y[] = 0.0;         // No vertical motion initially
    }
  }
}

Physical properties are specified based on the dimensionless parameters:

// Density ratio (liquid : gas = 1000:1)
rho1 = 1., rho2 = 1e-3;
  
// Gas phase Ohnesorge number
Oha = 2e-2 * Oh;
  
// Viscosity from Ohnesorge numbers
mu1 = Oh/sqrt(We), mu2 = Oha/sqrt(We);
  
// Polymer relaxation time (liquid has elasticity, gas does not)
lambda1 = Wi; lambda2 = 0.;
  
// Elastic modulus (liquid is viscoelastic, gas is not)
G1 = El; G2 = 0.;
 
// Surface tension coefficient
f.sigma = 1.0/We;

Scalar vs. Tensor Implementation The code provides two options for implementing viscoelasticity:

  1. Tensor Implementation: Full anisotropic stresses, physically accurate but computationally intensive
  2. Scalar Implementation: Simplified model, faster but less accurate for complex flows

The tensor implementation tracks the full conformation tensor, while the scalar implementation uses a simplified approach that captures the main features while being computationally more efficient.

4.3 Expected Results

Viscoelastic drop impact exhibits several distinctive behaviors compared to Newtonian impact:

  1. Weissenberg = 0 (Newtonian): Standard impact dynamics with spreading, potential rebound depending on Oh and We
  2. Low Weissenberg: Enhanced rebound due to elastic recovery
  3. Moderate Weissenberg: Formation of extended filaments during rebound
  4. High Weissenberg: Dramatic enhancement of rebound with possible beads-on-a-string formation in filaments

Exercise: Elasticity Effects on Impact

Investigate how elasticity affects drop impact behavior:

  1. Run simulations with Wi = 0 (Newtonian), Wi = 1 (moderate elasticity), and Wi = 10 (high elasticity)
  2. For each case, track:
    • Maximum spread factor
    • Rebound height (if applicable)
    • Time evolution of kinetic and elastic energy
  3. Visualize the stress distribution during impact for different Wi values
  4. Explain the physical mechanisms behind the observed differences

This exercise will help you understand how elasticity modifies impact dynamics through storage and recovery of elastic energy.

Part 5: Taylor-Culick Viscoelastic Retraction

Our final simulation explores the capillary-driven retraction of a viscoelastic fluid filament, known as Taylor-Culick flow.

5.1 Physical Background

When a cylindrical fluid filament with a free surface is suddenly released, surface tension causes it to retract. For Newtonian fluids, the retraction velocity reaches a constant value known as the Taylor-Culick velocity:

where is the surface tension, is the density, and is the filament thickness.

For viscoelastic fluids, the retraction dynamics are more complex due to elastic effects:

  1. Initial acceleration phase similar to Newtonian fluids
  2. Development of a bulbous end that grows as it collects fluid
  3. Formation of a neck behind the bulbous end
  4. Potential development of beads-on-a-string morphology
  5. Modified retraction velocity due to elastic stresses

Problem Statement Simulate the retraction of a viscoelastic fluid filament under capillary action, examining how elasticity affects the retraction dynamics. Key dimensionless parameters include:

  • Deborah number (De): Ratio of relaxation time to capillary time
  • Ohnesorge number (Oh): Ratio of viscous to inertial and capillary forces
  • Elasticity number (El): Ratio of elastic to inertial effects

5.2 Implementation Details

The simulation uses the axisymmetric framework with log-conformation approach for viscoelasticity:

#include "axi.h"                  // Axisymmetric coordinates
#include "navier-stokes/centered.h"  // Main Navier-Stokes solver
 
// Viscoelastic model implementation
#if !_SCALAR
#include "log-conform-viscoelastic.h"  // Tensor implementation
#else 
#include "log-conform-viscoelastic-scalar-2D.h"  // Scalar implementation
#endif
 
#define FILTERED                           // Smear density and viscosity jumps
#include "two-phaseVE.h"                   // Two-phase flow with viscoelasticity
#include "navier-stokes/conserving.h"      // Conservative momentum advection
#include "tension.h"                       // Surface tension implementation

The initial configuration consists of a cylindrical filament with a hemispherical cap:

event init(t = 0){
  if(!restore (file = dumpFile)){
    // Refine grid near the interface
    refine(x<(h0/2.0)+0.025 && y < hole0+(h0/2.0)+0.025 && level<MAXlevel);
    
    // Initialize the interface geometry (cylindrical filament with cap)
    fraction(f, y < hole0+(h0/2.0) ? 
        sq(h0/2.0)-(sq(x)+sq(y-(h0/2.0)-hole0)) : (h0/2.0)-x);
  }
}

Length Scales In this simulation, there are two important length scales:

  1. h0: The filament thickness, set to 1.0 (the characteristic length)
  2. hole0: The initial distance of the hemispherical cap from the origin, also set to 1.0

These dimensions determine the initial geometry of the filament, which affects the subsequent retraction dynamics.

5.3 Expected Results

The retraction dynamics of viscoelastic filaments show distinctive behaviors:

  1. Newtonian (De = 0): Formation of a bulbous end with diameter approximately twice the filament thickness, constant retraction velocity

  2. Low Deborah Number: Similar to Newtonian but with slightly modified bulbous end shape and retraction velocity

  3. Moderate Deborah Number: Formation of a neck behind the bulbous end, potential development of beads-on-a-string morphology

  4. High Deborah Number: Dramatic modification of retraction dynamics, with strong elastic resistance to extensional deformation in the neck region

Exercise: Viscoelastic Filament Retraction

Investigate how viscoelasticity affects filament retraction:

  1. Run simulations with De = 0 (Newtonian), De = 1 (moderate elasticity), and De = 10 (high elasticity)
  2. For each case, track:
    • Retraction velocity over time
    • Bulbous end diameter
    • Neck formation and evolution (where applicable)
  3. Compare the observed retraction velocity with the theoretical Taylor-Culick velocity
  4. Explain how elasticity modifies the retraction dynamics through extensional rheology

This exercise will help you understand how elasticity affects capillary-driven flows, particularly through resistance to extensional deformation.

Part 6: Numerical Considerations for Viscoelastic Simulations

Simulating viscoelastic flows presents unique numerical challenges that require special attention:

6.1 The High Weissenberg Number Problem

As mentioned earlier, traditional numerical methods often fail at high Weissenberg numbers due to the loss of positive-definiteness of the conformation tensor. This manifests as:

  1. Unbounded growth of stress values
  2. Numerical instabilities
  3. Simulation divergence

The log-conformation approach addresses this by:

  1. Working with the logarithm of the conformation tensor
  2. Ensuring the reconstructed tensor remains positive-definite
  3. Allowing stable simulations at much higher Weissenberg numbers

6.2 Adaptive Mesh Refinement Strategies

Effective adaptive mesh refinement is critical for viscoelastic simulations. Key refinement criteria include:

  1. Interface location: Captured by volume fraction gradients
  2. Velocity gradients: Important for capturing flow features
  3. Interface curvature: Critical for surface tension effects
  4. Conformation tensor gradients: Essential for capturing elastic stress variations

Balancing these criteria is important for efficiency and accuracy.

6.3 Multiphase Interface Treatment

For multiphase viscoelastic simulations, proper treatment of interfaces is crucial:

  1. Filtered approach: Smooths property jumps across interfaces for numerical stability
  2. Conservative formulation: Ensures momentum conservation during rapid interface deformation
  3. Tension model: Accurately captures surface tension forces

6.4 Stability Tips

Based on experience with these simulations, here are some tips for ensuring stability:

  1. Use a sufficiently small CFL number (typically 0.1 or less)
  2. Start with a moderate Weissenberg/Deborah number and gradually increase
  3. Ensure sufficient mesh resolution near interfaces and regions of high stress gradients
  4. Monitor energy conservation and stress growth for early signs of instability
  5. Use the tensor implementation for complex flows, switching to scalar only when necessary for performance

Debugging Strategy If your viscoelastic simulation becomes unstable:

  1. Check the kinetic energy for exponential growth
  2. Examine the conformation tensor components for unbounded growth
  3. Reduce the time step (lower CFL number)
  4. Increase mesh resolution near problematic regions
  5. Consider reducing the Weissenberg/Deborah number temporarily
  6. Ensure boundary conditions are properly implemented

Conclusion

In this assignment, we’ve explored the fascinating world of computational viscoelastic fluid dynamics through four complex multiphase flow simulations. These simulations demonstrate the rich and often counterintuitive behaviors that emerge when elastic and viscoplastic effects are introduced into fluid flows.

Key insights from our exploration include:

  1. The critical importance of specialized numerical techniques (like the log-conformation approach) for stable viscoelastic simulations
  2. The profound impact of elasticity on interfacial dynamics, from enhanced jet formation in bubble bursting to modified retraction behavior in filaments
  3. The formation of unyielded regions in viscoplastic flows, creating complex solid-fluid transitions
  4. The challenges and strategies for adaptive mesh refinement in flows with multiple time and length scales

These complex fluid simulations represent the culmination of our progression through computational fluid dynamics, building on earlier foundations in heat conduction, single-phase flows, and simple multiphase systems. The techniques and insights gained are applicable to a wide range of real-world phenomena, from polymer processing and food science to biological flows and advanced manufacturing.

Further Exploration

Consider these questions for deeper investigation:

  1. How would viscoelasticity and viscoplasticity interact in a combined model (e.g., an elasto-viscoplastic fluid)?
  2. What role does fluid microstructure play in determining macroscopic rheological behavior?
  3. How could these models be extended to incorporate temperature-dependent properties?
  4. What experimental techniques could validate the computational predictions made by these simulations?
  5. How might machine learning approaches be used to develop more accurate constitutive models for complex fluids?

Back to main website Home, Team, Research, Github

Graph View