Solid mechanics simulations are mainly based on the finite element method that adopts the Lagrangian description—in which computational mesh is carved in the material body and deforms with the body. In post-processing, analysts naturally form their interpretation from results at Lagrangian mesh, for example, displacements and reaction forces at nodes, stresses and strains at integration points, nodes, or centroids of elements. Those results are attached to specific material points through Lagrangian mesh and are called results in Lagrangian representation. Investigating Lagrangian results satisfies analysts’ needs in solid mechanics problems most of time.

But in some dynamic or quasi-static problems in which solid material behaves like a flow, we may occasionally be more curious about results at some specific spatial locations rather than material points. For example, in a roll-to-roll process where a web travels along a designed path line for some functional processing from a unwinding roll to a rewinding roll, engineers may specifically be interested in the tension, traveling speed, temperature, or other physical quantities of web at some spatial locations such as some rollers since they may serve as parameters of control system or relate to product quality. Results of material expressed at spatial points are called results in Eulerian representation.

Abaqus has the ‘probe’ function in the visualization module which enables users to approximately get results at spatial locations through Lagrangian mesh in individual solution frame. However, it does not provide a convenient way of extracting the time history of results at spatial locations throughout the whole simulation. Yet, the time history of Eulerian results in these simulation is more valuable as it reveals when and how the steady state is achieved, which some times exposes physical information about the characteristics of the problem.

This project develops two general-purpose Abaqus Python scripts based on two methods that can automatically extract the time history of field outputs at some spatial location of interest. The developer wishes the scripts can save some labor and willpower of analysts so their energy can be reserved for pure problem solving rather than tedious external processing.

Implementation details refer to this blog post.

If a web laying on a plane surface does not have a perfectly flat and straight shape, it has camber. The cambered web laying on the plane surface that flattens but forms a curved shape has a uniform camber; the cambered web that has does not flatten on the plane and forms a baggy lane or edges has local cambers [1]. Cambered webs have severe steering issues as handled on the web line.

We developed an implicit dynamic FE model and studied the deformation of cambered web with uniform curvature transporting on web line and the boundary conditions in the final steady state. The web camber was modeled by an initial linearly varied thermal expansion over width of the web, which was implemented by a spatial temperature field defined by the Abaqus user subroutine UTEMP. We found cambered web steers to its longer side and ultimately forms a S-shape in the span.

The reason of using implicit dynamic analysis type over explicit dynamic is, first, implicit dynamic analysis provides cleaner displacement results while the default explicit dynamic analysis has more noise; second, the subroutine which directly assigns a spatial temperature field is only available in Abaqus/Standard. There is no quick way to implement this in Abaqus/Explicit.

Lateral displacement and lateral stress of cambered web transporting on two aligned rollers

Sometimes in a modeling and analysis work, successfully simulating a problem is not easy, successful interpretation of results is harder. The latter requires a solid theoretical knowledge, a close experience of the physical problem,  and the intuition formed by both. As a thorough understanding of problem is gained, articulation and representation of the understanding is also hard. Most of the finite element codes are focusing on improving the algorithm part, which is tries to provide a better tool to simulate more physical problems, accurately and efficiently. The post-processing part of those codes, which helps to articulate and represent our understanding, lags far behind. In this problem I developed a little Python script that can extract relevant data from ODBs of Abaqus and visualize them. Most of the time for the cambered web problem, the middle line of the cambered web can represent the movement and deformation of the web. So I specifically visualized the lateral displacement, slope, and curvature of the middle line of the cambered web between the two rollers (blue and red in the animated figures below). I think this helps us to intuitively “see” the deformation and boundary conditions related with beam theory.

Visualization of the lateral displacement, slope, and curvature of the middle line of the cambered web

We also simulated the cambered web transporting on 3 rollers that steers on 2 spans and the cambered belt problem.

A cambered belt steering on 2 spans

Steering of a cambered belt

1. J. J. Shelton, “Effects of Web Camber on Handling”, Proceedings of the Fourth International Conference on Web Handling, Stillwater, OK, 1997

The Wound-On-Tension (WOT) is the tension in the outermost lap of a winding roll during winding. It is the most influential input that determines the wound roll stresses. WOT is dependent of winder types, material properties of web, winding parameters, and geometry of winders. We developed an explicit model to account for all those aspects. Winder types of pure center winder, center winder with nip, surface winder, and hybrid winder were studied. Web material properties are known to be state dependent. This is accommodated by implementing a user subroutine VUMAT in Abaqus. New test methods were developed for out-of-plane shear modulus. Experimental verifications of the WOT were conducted and good agreement was found. A parametric model was developed by Python scripts, which facilitates use of the model by other analysts.

Stress components of a roll in pure center winding

Stress components of a roll in center winding with a nip roller

Measurement of out-of-plane shear modulus of web stacks

One of the fundamental rules in web handling is a web always tends to enter a cylindric roller normally—the middle of web will be perpendicular to the axis of the roller. When a roller is misaligned in a web line, the web will steer laterally on the roller until a normal entry condition is achieved, thus, induce a bending deformation in the span of web prior to the misaligned roller. The shear forces developed between web and the misaligned roller have the potential to cause instability of web in spans or on rollers. We used explicit finite element analysis to investigate the instability of web due to a misaligned roller.

Simulations showed the evolution of instability of web due to a gradually misaligned roller. First troughs form in the span. When misalign angle achieve certain value, a micro wrinkle enters the misaligned roller and moves across it, forming a wrinkle on the roller. Last, a foldover formed on the roller. The defect on the web then transported through the whole virtual web line downstream. This evolution is exactly consistent with the lab observation.

A wrinkle and foldover formed and passed through the web line.

A closer look of at the simulation results reveals that the misaligned roller induces the development of a pocket of negative lateral stress prior to the web’s entry point to the roller, which causes the formation of a micro wrinkle. The wrinkle then climbs on and over the roller, a foldover ultimately forms.

Development of the negative lateral stress during the evolution of troughs – wrinkles – foldovers

The discontinuous lateral displacement of web due to wrinkles

This work was later extended to help the development of a failure criterion that predicts the critical buckling stress and critical misaligned angle for web instability in a web line. The criterion provides a guidance for web line design and assembly and defect tracing.

The key of successful modeling this problem is the connector constraint and boundary condition that depict the movement of the misaligned roller. If they are not set up right, we may see some really odd behaviors like this one, which may be viewed as the problem of a loose roller or a special process?

Misaligned roller is swaying about its center point.

A nip set consisting of two nip rollers is widely used on web lines to adjust web line tension, feed in inks for printing, conduct lamination, and etc. The nip rollers may have cover material such as rubber or foam. Here we are curious how the velocities and tension are changed by the different configuration of nip set. Specifically, we studied two types of the problem: the type-1 problem has the tension control on both upstream and downstream ends of the nip couple; the type-2 problem has the velocity controlled upstream and tension controlled downstream.

Two types of problems

We developed an implicit dynamic Abaqus model and studied the problem of a web transiting through a nip roller set. We studied different nip couple configurations: rubber-covered nip roller against rigid roller and rubber-covered nip roller against rubber-covered nip roller. We investigated how the nip loads and rubber material properties affect the tension and velocity upstream and downstream of a nip set for two types of problems with two nip configurations. This study provides additional detailed information for classic tension control theory with nip couples. The contact status (slipping / sticking) developed between nip rollers and the web has insights to the more accurate control in printing and lithography applications.

The maximum principle stress developed in nip-cover

The contact status: slipping / sticking between nip rollers and web

Previous methods used to measure surface profiles in web handling area are hard to obtain high accuracy and efficiency simultaneously. We explored the fringe projection technique and developed a system that used a commercial projector and a digital camera and produced non-contact, full-field measurement of surface profiles with high accuracy. We applied the system to measure and reconstruct the troughs of webs and the side surface of rolls. The system has a good potential to be used in more areas of surface profiles measurement.

We reported this work in the 10th International Conference on Web Handling and won the best paper
award. Details can be found in the presentation slides iweb09.

Measurement of web troughs due to fixed end tension

Measurement of roll side surface

Date: 11.2008
My role: modeling analyst and experimenter

The in-plane orthotropic material model works quite well in capturing the elastic behavior of spunbond nonwoven under small deformation¹. Although the model is a basic one in elasticity, some material constants such as the Poisson’s ratio and shear modulus are difficult to measure accurately for nonwoven. We proposed a practical material calibration way that is based on fitting the deformed shape of samples in uniaxial tensile tests with an inverse FEA procedure, to avoid the direct measurements of those parameters. This method transfers some of the experimental effort into analysis and provides an accurate material calibration for modeling of nonwovens.

The in-plane orthotropic material model has 4 independent constants, which can take the form of two Young’s moduli \(E_1\) and \(E_2\) in two respective principle directions—machine direction (MD) and cross machine direction (CMD), one Poisson’s ratio \(\nu_{12}\) (called major Poisson’s ratio) relating the induced deformation in CMD to the loaded deformation in MD, and one in-plane shear modulus \(G_{12}\). The strain-stress relation connected by the material compliance matrix is

$$\begin{bmatrix}\epsilon_{11} \\ \epsilon_{22} \\ \gamma_{12} \end{bmatrix}=\begin{bmatrix} \frac{1}{E_1} & -\frac{v_{12}}{E_1} & 0 \\ -\frac{v_{12}}{E_1} & \frac{1}{E_2} & 0 \\ 0 & 0 & \frac{1}{G_{12}}\end{bmatrix} \begin{bmatrix}\sigma_{11}\\\sigma_{22}\\\tau_{12}\end{bmatrix}.$$

Among the 4 material constants in the model, the two Young’s moduli can be measured by uniaxial tensile tests and have the best accuracy and precision. The in-plane major Poisson’s ratio of a typical nonwoven web is usually grater than 1 due to their composed microstructure. The severe necking deformation in a traditional uniaxial Poisson’s ratio test makes the calculation of transverse strain confusing. Thus the value of major Poisson’s ratio obtained may not be reliable. Furthermore, an accurate shear modulus test for nonwoven is even harder because of its extremely low resistance to shear buckling; thus, it’s difficult to gather enough meaningful data in a pure shear test of nonwoven before it buckles. A shear test using the picture frame may be viable but it causes more effort in machining clamps and preparation of samples.

To avoid the uncertainties of measurements on Poisson’s ratio and shear modulus, one more uniaxial test in some other angle—\(\theta\) besides the two principle directions is conducted. The orthotropic theory gives the expression of modulus in the \(\theta\) direction

$$\frac{1}{E_\theta}=\frac{\cos^4(\theta)}{E_1}+\frac{\sin^4(\theta)}{E_2}+\frac{1}{4}\bigg[ \frac{1}{G_{12}}-\frac{2\nu_{12}}{E_1} \bigg]\sin^2(2\theta).$$

Therefore, knowing the measured values of \(E_1\), \(E_2\), and \(E_\theta\) in \(\theta\), the shear modulus \(G_{12}\) and major Poisson’s ratio \(\nu_{12}\) are related by the above equation. There is only one variable left need to be determined. This allows us to transform the direct measurement problem into a single variable optimization problem—with the aid of finite element analysis used in an inverse and iterated fashion. The optimization objective is set to get an agreed deformed shape of test samples and FEA. The specific procedure is as follows:

1. Three uniaxial tensile tests up to 1% strain are conducted in MD, CMD, and \(\theta\) direction. The three Young’s moduli were measured. And the deformed shape of sample in the MD test were recorded.
2. A finite element model which simulates the MD tensile test of the same dimension sample is created. The material model in simulation is set as orthotropic elastic model with inputs of measured Young’s moduli \(E_1\) and \(E_2\), an assumed major Poisson’s ratio \(\nu_{12}\), and the calculated shear modulus \(G_{12}\) according to the above equation.
3. A python script was developed to run multiple simulations, iterating Poisson’s ratio and corresponded shear modulus until the deformed shape of FEA matches that of the test. Then the complete material model is calibrated.

In seeking an agreed deformed shape between FEA and experiments, we’ve practiced using the lateral contraction only at the middle of the sample as the criterion. A more accurate calibration would use digital image correlation (DIC) to have a full field deformation measurement of nonwoven sample and minimize an average measure of displacement difference between FEA and test over multiple selected control points.

In summary we proposed and verified a practical material calibration way to fit parameters of orthotropic model for nonwoven. The method requires only three uniaxial tensile tests and avoids the uncertainty in measuring Poisson’s ratio and shear modulus. Since the method is based on an inverse FEA procedure of which the objective is to match the test deformation, the calibrated material model is especially accurate for future finite element analyses of nonwoven under similar tensile loadings. Nevertheless, it should be noted the proposed method is based on the validity of orthotropic model on depicting the elastic behavior of nonwoven, which is only validated in small deformation. If one aims to model behavior of large deformation where the material nonlinearity exhibits, more advanced constitutive models need to be developed. The orthotropic framework may still be valid, but material parameters must be made state dependent.

1. H.S. Kim, “Orthotropic theory for the prediction of mechanical performance in thermally point-bonded nonwovens”, Fibers and Polymers,(2004) 5: 139. https://doi.org/10.1007/BF02902928 ↩

Date: 06.2007–06.2008
My role: main modeling analyst

Rotary draw bending is the metal forming process widely used in aerospace and automotive industries to form straight metal tubes into designed spatial geometric shapes. Viewed from the perspective of analysis, the process is severely nonlinear from multiple sources—the geometric nonlinearity due to the large rotation, the material nonlinearity due to metal plasticity, and the state nonlinearity due to complex contact interactions between dies, mandrels and tubes. An accurate solution that can predict the defects of material during forming and capture the geometric features of tubes after bending as well as after springback is essential. This solution is only viable through finite element analysis.

We used the finite element code LS-DYNA as the core solver. In the bending process the tube’s kinetic energy produced by motion is negligible compared to its internal energy by deformation, the forming process is quasi-static, thus was modeled by the explicit solver so the convergence issue caused by nonlinearities could be avoided. The springback of bent tubes was computed by the implicit solver, which was realized by LS-DYNA’s feature of seamless transfer from explicit to implicit.

Movements of parts

Belytschko-Tsay reduced integration shell elements were used in the forming simulation to discretize the whole model. Stiffness-based hourglass control was used to inhibit hourglass mode. The reduced integration elements were switched to full integration after the model was transferred into the implicit solver so that stresses during the unloading is more accurately computed.

Mesh

The first pack of tubes were made of a kind of DQAK steel, whose complete elasto-plastic behavior was modeled. The plasticity was modeled as the power law isotropic strain hardening with Von Mises yield criterion. All the dies, clamps, and mandrels are modeled as rigid bodies to save computational cost; yet their geometric features are important to the forming shape of tubes, so their geometric features were preserved as much as possible in the model and a high density mesh were used for small features like fillets. Mechanisms between the last mandrel ball and mandrel axel, and between mandrel balls themselves were models by spherical joints type of kinematic constraints. Both time scaling and mass scaling were used to accelerate computation, the sensitivities to results were carefully tested,.

The analysis were not only developed to investigate the bending process for the DQAK steel, they were also intended to serve as a predictive tool that would be ultimately integrated in the design and optimization process for other material and geometry demands. An automatic means of pre-processing would save significant modeling effort. Thus, a parametric model was created to cover a wide range of scenarios and was implemented by ANSYS APDL. The abstracted model has 66 parameters to completely describe the geometry, material, contact, and processing aspects of the problem. Hence, the pre-processing labor left was to input values of those parameters. Then the FEA model would be automatically created and the corresponding keyword file of LS-DYNA would be submitted for computation.

In post-processing, we analyzed tube’s stress/strain, thickness variation, section distortion, and wrinkling under certain working conditions. With the aid of the parametric modeling, we conducted stacks of computation to optimize processing parameters for specific targets, for example, the rotary angle of bend die according to the final formed angle of tube after springback.

Thickness variation during the process

Springback after clamps and dies are released