A polymer nanocomposite reinforced with single walled carbon nanotube is simulated using the LAMMPS molecular dynamics code in order to obtain elastic properties. The simulation volume contains 10 polymer chains that are 260 monomers long surrounding one central carbon nanotube containing 930 atoms. Due to simulation restrictions, the atoms in the carbon nanotube are treated as a rigid body. The volume is modeled as periodic in all directions.
The rigid body restriction on the nanotube prevents calculation of longitudinal properties, therefore displacements are only applied normal to the nanotube axis in order to obtain transverse properties. The transverse properties are not directly affected by the nanotube, so a rigid nanotube is a reasonable approximation. After each displacement is applied a time-averaged value of stress is calculated. LAMMPS provides the stress contributions σ_ij^α from each atom as where superscript α corresponds to the atom, V is the atomic volume, v_i and v_j are the velocity components in their respective directions, and F_i^αβ and r_j^αβ are the components of force and separation between atoms α and β respectively. The stresses used to cerate the stress-strain curves are the average atomic stresses across the volume of the MD model, Here V is the simulation volume. These values of stress and strain are used to calculate Young's modulus using Hooke's law.

A uniform tensile strain is applied to the simulation volume normal to the axis of the nanotube. Strains are applied by resizing the simulation volume in the direction of the strain and rescaling the coordinates of the atoms to fit the new dimensions. Each iteration increases the length of the volume by 0.05 angstroms. The simulation is allowed to run until reaching 100% strain, meaning a displacement of 40 angstroms (4 nm) normal to the nanotube. After each iteration the system is allowed to equilibrate for 2,000 time steps, then the reported stresses for the next 1,000 time steps are averaged to produce a single value for stress at the current strain. Initially the response is linear-elastic and the stress-strain plots are used to calculate Young's modulus, according to Hooke's law. As the strain increases the polymer begins to fail and the stress-strain curves break down.

The simulation is repeated three times with different interface strengths between the polymer and the nanotube. The interface strength is based on the standard 12/6 Lennard-Jones potential, where ε is the depth of the potential well, σ is the distance at which the interparticle potential is zero, and r is the distance between the particles. For this paper the interface strength is adjusted using three different values of ε; ε_std=0.00286 eV, ε_weak =0.25ε_std=0.000715 eV, and ε_str=4ε_std=0.01144 eV.

The stress-strain curves for the polymer-nanotube composite are included in Figure 2. The weakest interfaces experiences several partial failures under relatively low strain as the polymer separates and regroups along the nanotube. The polymer-nanotube interface is not strong enough to constrain the polymers, which mass together at one end of the nanotube and begin to separate from it as the strain increases. As the strain continues to increase the polymer separates completely from the nanotube and clusters in the gap between. With the standard strength interface the polymer continues to be evenly distributed along the nanotube, but it separates from the polymers in the periodic volumes on either side. Eventually only one polymer chain remains to hold the system together. The strong interface is very similar, except that the polymer chains remain even closer to the nanotube.
For relatively low strain the response of each nanocomposite is linear-elastic, therefore Young's modulus can be calculated via Hooke's law. For the system with the weak interface three values of Young's modulus were calculated, corresponding to the initial slope of the stress-strain curve and the slopes after the first two partial failures. Averaging these three values results in an approximate Young's modulus of 3.67 GPa. The standard strength interface and the strong interface were calculated using only the linear portion of the graph, resulting in Young's moduli of 4.36 GPa and 4.22 GPa, respectively. The strong interface may have a lower elastic modulus than the medium interface because the polymer chains are constrained too closely against the nanotube. This prevents them from occupying the empty space and tying the simulation volumes together.

References

1. Frankland, SJV and Harik, VM and Odegard, GM and Brenner, DW and Gates, TS, The Stress-Strain Behavior of Polymer-Nanotube Composites from Molecular Dynamics Simulation. Composites Science and Technology, volume 63, pages 1655-1661, 2003.

Go here to request the source code.

A level-set based image segmentation method based on a paper by Li et al. [1] has been implemented. The method eliminates the need for reinitialization by adding energy terms which force the level set to remain close to a signed distance function. The evolution equation is given by: where μ, λ, and ν are user-defined parameters. The μ term is an internal energy term that depends only on φ and penalizes the level set when it deviates from a signed distance function. The λ and ν terms are external energy terms that drive the motion of the level set towards the desired image features. Here g is an edge indicator function given by where G_σ is the Gaussian kernel and I is the image data. The δ_ε is a smoothed dirac delta function

Two solution methods for this system of equations are compared. In the explicit solution time is discretized using the forward Euler method and solving for the next time step φⁿ⁺¹, The semi-implicit solution method is more complicated. The discretized equation is which can be rewritten as This is the standard form of a Helmholtz equation. Rather than write a custom solver to solve this system of equations, we use Intel's Math Kernel library Helmholtz solver in two dimensions. Spatial derivatives for both the explicit and the implicit method are solved using central differences.

The explicit and semi-implicit methods produces very similar results, but the semi-implicit scheme has a significantly broader range of stable input parameters. The semi-implicit method was able to run stably with a Δτ of 60.0 for some problems, while the explicit method never ran stably with Δτ above 10.0. Larger time steps allowed the semi-implicit simulations to complete in a fraction of the time.

Two sample solutions are shown below, where blue lines are partially resolved zero-contours of the level set, and the final zero-contour is shown in red. The initial level set for both problems was defined as a rectangle near the edge of the image, where φ(x)=-6.0 for all x inside the level set and φ(x)=6.0 for all x outside the level set. Clearly the method is capable of evolving the level set around complex image contours whether the edges are sharp or ill-defined, even without reinitialization. In the second example, varying the parameters allows us to capture different contours of the image.

References

1. Li, C. and Xu, C. and Gui, C. and Fox, M., Level Set Evolution Without Re-initialization: A New Variational Formulation. IEEE Computer Society Conference on Computer Vision and Pattern Recognition, volume 1, pages 430-436, 2005.

See the full pdf here, or go here to request the source code.

Pressure-dependent materials such as soil or concrete require advanced material models to predict yielding. A commonly used yield criterion is the two-invariant Mohr-Coulomb model, as described in (1.1), where the σ_i's correspond to principal stresses, c is the cohesion, and φ is the friction angle. On the π-plane, the Mohr-Coulomb yield surface is similar to a distorted Tresca yield surface with lower tensile strengths and higher compressive strengths. It is well suited to any porous or granular material in compression loading. The corners of the Mohr-Coulomb yield surface visible on the π-plane, as well as the pressure vertex in the meridian plane, can cause problems for numerical integration techniques. Several different smoothed approximations have been developed in order to simplify calculations, and it has been shown that these methods may actually be more accurate than the unmodified Mohr-Coulomb model.

This paper will derive the Drucker-Prager smoothed approximation to the Mohr-Coulomb plasticity model, including a detailed account of the integration algorithm and numerical implementation. The Drucker-Prager yield surface and plastic potential are defined by with deviatoric stress and mean normal stress (pressure) , with and δ as the Kronecher delta. The hardening laws are given by Here λ is the cumulative plastic multiplier and α₀, β₀, a, and k are material parameters.

The Drucker-Prager material model is non-linear, so it is integrated with an iterative numerical scheme. In this paper a simple Newton-Raphson scheme is used in the material subroutine. Although the Drucker-Prager method only depends on two invariants, in this paper the more general three-invariant return mapping scheme is used. This is more difficult to implement, but can be applied to a broader range of material models. The method is a spectral method and therefore frame-invariant, and so can be easily applied to a broader range of material loading conditions. It requires the repeated calculation of eigenvalues and eigenvectors which can become computationally expensive in large models, but for the small numerical models presented here this is irrelevant.

Numeric Examples

Two numerical experiments are discussed. These are designed to test both the global loading iteration loop and the local material subroutine iteration loop. The common initial stress state is p₀=-50 kPa and q=0 kPa, ie a hydrostatic state. Both experiments use vertical strain control and horizontal stress control. During loading, the vertical strain increment is Δε₁ = -0.2% and the horizontal stress is kept constant. Both experiments stop at ε₁ = -3%.

Axisymmetric Loading

For the case of axisymmetric loading, a strain Δε is applied along the axis while the radial stresses σ_r are kept constant. In order to apply these complicated loading conditions, a global Newton-Raphson loop is implemented. The Newton-Raphson iteration scheme solves for the appropriate stresses and strains to meet the boundary conditions. We can write this as a residual vector R and a vector of unknowns X, This problem is difficult to write a residual for, so we rewrite it Now the Newton iterations can be defined by where C is the consistent tangent operator from the material subroutine. Because the problem is strain-controlled, Δε₁ must equal the -0.2% loading. In order to ensure this, before solving an iteration step we back-substitute and solve for a σ₁ that satisfies the Δε₁ boundary condition. Now we can solve the iteration step for Δε₂ and Δε₃.

Under this loading, the Drucker-Prager is elastic for one loading step before going plastic. Once the material response is plastic, the material begins to soften. The convergence for this experiment is quadratic, as expected from the use of the consistent tangent operator. In the final step the iterations converge to 1e-16, which is machine precision, not a reflection of the method or the implementation.

Plane Strain Loading

For the plane strain loading, a strain Δε is applied in the vertical (ε₁) direction. The strain ε₂ is maintained at zero, and only the stress is controlled in the third direction. Again, a global Newton-Raphson iteration scheme is used to apply these conditions. The same Newton-Raphson loop as in axisymmetric loading is used here, with a minor modification to control two strains instead of just one. Here we solve for the σ₁ and σ₂ values that correspond to the desired strain state before inverting the matrix.

The results are similar to the case of axisymmetric loading. The additional strain constraint in the horizontal direction from plane strain causes higher stresses under loading. Convergence for the plane strain loading is still quadratic, due to the consistent tangent operator.

References

1. See the full paper here.
2. R.I. Borja and J.E. Andrade. Critical state plasticity. Part VI: Meso-scale finite element simulation of strain localization in discrete granular materials. Computer Methods in Applied Mechanics and Engineering, 195(37-40):5115-5140, 2006.
3. R.I. Borja, K.M. Sama, and P.F. Sanz. On the numerical integration of three-invariant elastoplastic constitutive models. Computer Methods in Applied Mechanics and Engineering, 192(910):1227-1258, 2003.

See the full pdf here.

The advection-diffusion equation models situations where a substance is both transported and diffused through a medium. The one dimensional advection-diffusion equation for a particle concentration φ can be written Where u is the transportation velocity and v is the diffusion coefficient, both of which are assumed to be constant throughout the domain.

Solving the advection-diffusion equation using the standard Galerkin method can lead to oscillatory solutions caused by spatial numerical instabilities. The Petrov-Galerkin method is designed to alleviate these problems by essentially adding a viscous damping force. The implementation here is the Streamlined Upwind Petrov-Galerkin (SUPG) method in one dimension.

A specific model is implemented to show the differences between the Galerkin and Petrov-Galerkin methods. The simulation is a one meter segment of a long tube filled with solution. There is a steady state concentration of 5% at x=0 and 20% at x=1. The velocity u = 2.0 m/s and the diffusion coefficient v = 0.025 m²/s.

Finite Element Solutions

The Galerkin finite element discretization follows from the general form of the advection-diffusion equation. Using the test function w(x) and integrating over the domain, From this, the weak form is For a uniform finite element mesh, the discrete equations become The Petrov-Galerkin formulation modifies the Galerkin test function by adding an additional term where γ=αΔx/2. In one dimension it is possible to select α such that the solution at the nodes is exact, which is based on the Peclet number, which will be shown to have a large effect on the stability of the solution. The modified test function leads to the same terms as the regular Galerkin formulation, plus additional Petrov-Galerkin stabilization terms, The fourth term depends on the second derivative of φ, but for linear shape functions this term drops to zero. So for linear shape functions in a uniform grid, the discrete form of the Petrov-Galerkin is If non-linear shape functions are used, the discarded fourth term must additionally be implemented.

It can be shown that the stability of the Galerkin method depends on the Peclet number. If the Peclet number is less than one, then the Galerkin discrete solution approaches the exact solution. For Peclet numbers greater than one, the Galerkin discrete solution becomes oscillatory and the Petrov-Galerkin formulation is required. To illustrate this, the example problem in the introduction is run with several different Peclet numbers. To achieve the different Peclet numbers the resolution of the mesh is varied, using 10, 20, 50, and 100 elements.

The 10 and 20 element meshes result in Peclet numbers of 4.00 and 2.00, respectively. As shown below, the Galerkin method is unstable for these meshes and oscillates badly as x increases. Due to the proper choice of α, the Petrov-Galerkin method achieves the exact solution at the nodes, although the choice of linear shape functions limits the accuracy between nodes.

For the finer resolution meshes (50 and 100 elements), the Peclet numbers are 0.80 and 0.40, respectively. The Galerkin method is stable for these values, and compares favorably to the Petrov-Galerkin and analytic solutions.

For all cases the Petrov-Galerkin calculates exact solutions at the nodes, due to the choice of α in one dimension. However, the finer meshes give a better resolution of the solution. For the Galerkin method, increasing mesh resolution results in more accurate values at the nodes.

SUPG Parameters

The choice of γ in the Petrov-Galerkin formulation plays a major role on the stability and accuracy of the solution. The γ term essentially adds an artificial viscosity to the solution to dampen oscillations. Incorrect choices of γ can lead to the solution being over-damped and inaccurate, or under-damped and oscillatory. In one dimension, γ can be calculated such that the solution is exact at the nodes.

The solution with several values of γ has been plotted below. With low damping values the solution remains unstable, and with high values the solution is damped so much that the solution is no longer accurate.

Non-uniform Mesh

The stability problems of the Galerkin formulation are spatial instabilities, caused by problems with the mesh. As shown in the previous Section, refining the uniform mesh eliminates the oscillations, at the cost of increased computational effort. Selectively refining the mesh around key areas allows the Galerkin solution to remain stable without unnecessary increased computation.

In the example problem, the mesh is refined near the spike in concentration. The mesh on x=[0.9,1.0] is ten times finer than the mesh from x=[0.0,0.9]. The Peclet number is calculated based on the element length in the refined area. For this simple problem and simple non-uniform mesh, this Peclet number still results in exact solutions for the Petrov-Galerkin formulation.

For a more adaptively refined mesh, the calculation of the γ becomes more difficult, and it may not be possible to calculate a γ resulting in the exact solutions.