Diffusion Calculations

NairnMPM can do diffusion calculations coupled with stresses and strains through concentration-induced expansion.

Diffusion Modeling

To be able to model diffusion in composite materials where different phases may absorb different amounts of solvent, all diffusion calculations are done in terms of a concentration potential for the solvent in the material. Diffusion, as derived by mass flow, can be written as:

      [math]\displaystyle{ {\partial c\over \partial t} = - \nabla \cdot J = \nabla \cdot \left({Dc\over RT} \nabla \mu\right) }[/math]

where c is concentration (e.g., moles per unit volume) and J is flux given by minus the bracketed term where D is the diffusion coefficient (in area/sec), μ is the chemical potential (in J/mol), R gas constant (in J/(mol-K)), and T is temperature. This modeling writes the chemical potential is terms of a concentration potential

      [math]\displaystyle{ \mu = \mu^{(ref)} + RT \ln \bar{c} \qquad {\rm where} \qquad \bar{c} = \frac{c}{c_{ref}} }[/math]

is a concentration potential. Here [math]\displaystyle{ \mu^{(ref)} }[/math] is reference chemical potential when [math]\displaystyle{ c=c_{ref} }[/math]. Taking gradient of the chemical potential and substituting into the diffusion equation leads to

      [math]\displaystyle{ {\partial \bar{c}\over \partial t} = \nabla \cdot D \nabla \bar{c} }[/math]

Because [math]\displaystyle{ \bar{c} }[/math] is dimensionless, it can be defined using any consistent units for [math]\displaystyle{ c }[/math] and [math]\displaystyle{ c_{ref} }[/math]. When modeling composite materials, [math]\displaystyle{ c_{ref} }[/math] is a material property. Thermodynamic equilibrium in composites is when all phases reach the same chemical potential. The goal of diffusion modeling therefore is to pick [math]\displaystyle{ c_{ref} }[/math] for each material such that all materials have the same [math]\displaystyle{ \mu^{(ref)} }[/math].

Once [math]\displaystyle{ c=c_{ref} }[/math] is picked for each material, all internal calculations are done using [math]\displaystyle{ \bar{c} }[/math]. Furthermore, concentration boundary conditions are expected to be in terms of dimensionless [math]\displaystyle{ \bar{c} }[/math] instead of in terms raw concentration.

Reference and Saturation Concentration

The initial implementation of diffusion in NairnMPM was to model diffusion of solvents into solids including expansion due to solvent (e.g., swelling of wood due to moisture content). Most results for such solvent effects are based on mass fraction solvent content and expansion coefficients are measured as swelling strain per weight fracture solvent. As a result, NairnMPM assumes [math]\displaystyle{ c }[/math] and [math]\displaystyle{ c_{ref} }[/math] are both mass fraction solvent and the reference state is set to the saturation limit for solvent content or [math]\displaystyle{ c_{sat} }[/math] (which is specified for each material type in basic material properties). By this model, concentration can never exceed saturation concentration meaning that the [math]\displaystyle{ \bar{c} }[/math] is restricted to the interval [0,1] and equilibrium conditions correspond to all particles being at the same [math]\displaystyle{ \bar{c} }[/math] (i.e., at the same fracture of their saturation concentration).

The saturation model is the default method for diffusion calculations but it can be changed by redefining [math]\displaystyle{ c_{ref} }[/math] with any other units. For example, one could model concentration in moles per unit volume (molar concentration). This model still uses [math]\displaystyle{ \bar{c} }[/math] but that potential can now exceed 1 or can have any non-negative value. For this approach to work in composites, [math]\displaystyle{ c_{ref} }[/math] can be any state but each material must define [math]\displaystyle{ c_{ref} }[/math] such that all materials have the same chemical potential at their reference concentrations. Note that despite this change in meaning of [math]\displaystyle{ c_{ref} }[/math], it is still entered using the csat material property.

Stress-Free Concentration

Diffusion effects are coupled to mechanics by changes in concentration causing shrinking or swelling resulting in changes in strain. Linear swelling is modeled as

      [math]\displaystyle{ \Delta\varepsilon = \beta(c-c_0) = \beta c_{ref}(\bar{c}-\bar{c}_0) }[/math]

where [math]\displaystyle{ \bar{c}_0 }[/math] is the stress-free concentration potential. It is the potential at which a material point would have zero strain. Note that aniostropic materials may shrink and swell differently in different directions.

Each simulation can pick any non-negative [math]\displaystyle{ \bar{c}_0 }[/math]. Any material points that do not specify a different concentration will be set to potential [math]\displaystyle{ \bar{c}_0 }[/math]. But note that if a material point starts with a different potential, the first time step will see an incremental strain of

      [math]\displaystyle{ d\varepsilon = \beta c_{ref}(\bar{c}-\bar{c}_0) }[/math]

which could potentially be a large strain that acts like an impact to the material point. It is usually better to start all material points at [math]\displaystyle{ \bar{c}_0 }[/math] and let them naturally develop swelling strains by diffusion from boundary conditions are ramping them up with a property ramp.

Activating Diffusion

In scripted files, diffusion is activated with the command

Diffusion (YesOrNo),<(conc0)>,<(nolimit)>

In XML input files, diffusion is activated with the <Diffusion> command, which must be within the <MPMHeader> element:

<Diffusion reference = '(conc0)' style='1' nolimit='(nolimit)'/>

where

• (YesOrNo) must be "Yes" or "No" to activate or not activate diffusion calculations. It can also be "solvent" to indicate diffusion of a liquid into another material. In XML input files, the presence of a <Diffusion> command activates diffusion (it can add a style='1' attribute).
• (conc0) is used to set a concentration potential for the stress-free concentration potential or to enter [math]\displaystyle{ \bar{c}_0 }[/math]. The default (conc0) is 0. All material points will be initialized to [math]\displaystyle{ \bar{c} }[/math] of (conc0) unless explicitly set to a different potential. Note that for historic reasons, XML input files sets [math]\displaystyle{ \bar{c}_0 }[/math] with a "reference" attribute. The use of that name should not be confused with [math]\displaystyle{ c_{ref} }[/math] for each material, which is set with the csat material property.
• (nolimit) is used to select the diffusion model. If entered as "limit" (or "0"), the calculations interpret [math]\displaystyle{ c_{ref} }[/math] in each material to be defining a saturation constant and to limit [math]\displaystyle{ \bar{c} }[/math] to the interval [0,1]. To change to a model where [math]\displaystyle{ c_{ref} }[/math] has another meaning and [math]\displaystyle{ \bar{c} }[/math] can be any non-negative value, enter the text "nolimit" (or "1"). In XML input files, the "nolimit" attribute must be an "0" for "limit" mode or and any integer greater than zero for "nolimit" mode. The default (when omitted) is the "limit" mode. The code output will say "0 <= c/csat <= 1" for "limit" mode and "c/csat >= 0" for "nolimit" mode.

By default, diffusion uses update methods analogous to FLIP methods used in mechanics. This update, however, sometimes results in concentration oscillations on particles within one cell. Diffusion simulations with oscillations can be improved by using periodic FMPM(k) for diffusion updates.

Note that diffusion models fluid transport through materials by transport methods nearly identical to those used to model poroelasticity. Because they share same methods, a simulation can activate diffusion (with above commands) or poroelasticity (with comparable Poroelasticity commands), but cannot activate them both. Any simulation, however, can combine diffusion or poroelasticity with thermal calculations and conduction. Note that when choosing archiving options, the terms "concentration" and "porepressure" are synonyms or either can be used and the archiving will store concentration terms for diffusion calculations or pore pressure terms for poroelasticity calculations.

Note that diffusion calculations are often improved by using FMPM(k) methods.

Diffusion Material Properties

Concencentration changes are coupled to stress and strains through concentration expansion coefficients defined for the materials. By default, all moisture expansion coefficients are zero which decouples diffusion and strains. By entering non-zero values, the coupling will occur. Isotropic materials have a single solvent expansion coefficient (beta) while anisotropic materials will have two or three solvent expansion coefficients for swelling differently in different directions (e.g., see orthotropic materials).

The rate of diffusion is controlled by the solvent diffusion constant defined for each material. Isotropic materials have a single solvent diffusion constant, D, as defined above (and entered as property D). For anisotropic materials, D is replaced by a diffusion tensor and material properties set two or three diffusion properties for diffusion in various material directions (e.g., see orthotropic materials).

Archived Concentrations

All internal calculations use [math]\displaystyle{ \bar{c} }[/math] or concentration potential, but when concentrations and concentration gradients are archived, they are converted to dimensioned values using:

      [math]\displaystyle{ c = c_{ref}\bar{c} = {\tt csat}*\bar{c} }[/math]

This conversion applies both to particle archives and to global archiving. The archived results will therefore have same units the simulation assumed when entering each material's csat property.

Diffusion Boundary Conditions

When diffusion is activated, you can set, the possible concentration boundary conditions are:

• You can set concentration on the grid.
• You can set a concentration flux on particle surfaces.
• Gradually apply initial concentration field to particles using a PropertyRamp Custom Task.
• You can set initial particle concentration when creating material points. Note, however, that initial concentrations different than the [math]\displaystyle{ \bar{c}_0 }[/math] will cause strains to immediately evolve toward the changed state. The net effect will be an instantaneous "impact" that might cause undesirable dynamic effects. Using a PropertyRamp Custom Task is preferred method to set initial concentrations.
• Rigid particles can provide moving concentration boundary conditions.