iMechanica - large deformation - Comments

Re:NLGEOM in ANSYS

shrimad — Fri, 01 Aug 2008 17:08:00 -0400

Biswajit,

Thanks for the reply.

If I consider a small problem; a square as an axisymmetric model (giving rise to a cylinder in 3D). The material model is viscoelastic. If I have a line on the top surface and bottom surface and press the top surface to about half the height with the lower surface fixed (contact analysis, no friction) I get a deformed rectangle. I do the same with ANSYS and with my FEA code in FORTRAN. What happens is:

Without NLGEOM ON, I get an exact match with my code and ANSYS, but the volume is not conserved.

With NLGEON ON, I get a slightly larger rectangle and the volume is conserved in ANSYS in this case.

I want to conserve the volume in my FEA code also, and hence I am trying to understand what modification do I need to do in my FEA code? My FEA code follows Newton Raphson Method and the same algorithm as ANSYS for viscoelastic stress calculations. I am stuck in my research because of this. Please guide me with whatever you know.

Thanks,

Shriram

Research Assistant
Structural and Multidisciplinary Optimization Lab.

University of Florida

Re: NLGEOM in ANSYS

Biswajit Banerjee — Thu, 31 Jul 2008 23:24:37 -0400

The NLGEOM command activates corrections for large rigid body rotations and translations when a small strain constitutive relation is used. It is also activated for large strain material models. For details on the implementation you can go to Chapter 3. Structures with Geometric Nonlinearities of the ANSYS Theory Manual.

Using NLGEOM should have no impact on whether the incompressibility condition is satisfied. Incompressibility is usually satisfied using a Lagrange multiplier approach which can lead only to approximate satisfaction of incompressibility in a numerical code. However, there are are stress update algorithms that satisfy incompressibility exactly, often by using an exponential map - see Weber and Anand, 1988 (or so).

A cut/paste job from the ANSYS manual follows.

-- Biswajit

3.3. Large Rotation

If the rotations are large but the mechanical strains (those that cause stresses) are small, then a large rotation procedure can be used. A large rotation analysis is performed in a static (ANTYPE,STATIC) or transient (ANTYPE,TRANS) analysis while flagging large deformations (NLGEOM,ON) when the appropriate element type is used. Note that all large strain elements also support this capability, since both options account for the large rotations and for small strains, the logarithmic strain measure and the engineering strain measure coincide.

3.3.1. Theory

Large Strain presented the theory for general motion of a material point. Large rotation theory follows a similar development, except that the logarithmic strain measure ((Equation 3–6)) is replaced by the Biot, or small (engineering) strain measure: (3–37)where: [U] = stretch matrix [I] = 3 x 3 identity matrix 3.3.2.

Implementation

A corotational (or convected coordinate) approach is used in solving large rotation/small strain problems (Rankin and Brogan(66)). "Corotational" may be thought of as "rotated with". The nonlinearities are contained in the strain-displacement relationship which for this algorithm takes on the special form:(3–38)where: [Bv] = usual small strain-displacement relationship in the original (virgin) element coordinate system [Tn] = orthogonal transformation relating the original element coordinates to the convected (or rotated) element coordinates

The convected element coordinate frame differs from the original element coordinate frame by the amount of rigid body rotation. Hence [Tn] is computed by separating the rigid body rotation from the total deformation {un} using the polar decomposition theorem, (Equation 3–5). From (Equation 3–38), the element tangent stiffness matrix has the form: (3–39)and the element restoring force is: (3–40)where the elastic strain is computed from: (3–41) is the element deformation which causes straining as described in a subsequent subsection.

The large rotation process can be summarized as a three step process for each element:

Determine the updated transformation matrix [Tn] for the element.

Extract the deformational displacement from the total element displacement {un} for computing the stresses as well as the restoring force .

After the rotational increments in {Δu} are computed, update the node rotations appropriately.

All three steps require the concept of a rotational pseudovector in order to be efficiently implemented (Rankin and Brogan(66), Argyris(67)).

Chandra, Thanks for the

shrimad — Thu, 31 Jul 2008 11:25:47 -0400

Chandra,

Thanks for the reply.

I couldn't find the technical details in the ANSYS help menu. By theory manual do you mean the 'help' section or it is something different? Please let me know.

Thanks once again.

Shriram

It uses a Hencky or

Chandra Veer Singh — Wed, 30 Jul 2008 22:22:52 -0400

It uses a Hencky or logarithmic strain measure, which is ln(U), where U is the right stretch matrix obtained by polar decomposition of the deformation gradient. For calculating Hencky strain, it uses the eigen values and eigen vectors of U. You can read more in ANSYS theory manual, there is a chapter on structural nonlinearity.

Chandra Veer Singh

New items on deformable dielectrics

Zhigang Suo — Sun, 02 Sep 2007 07:35:20 -0400

I posted a theory of polyelectrolyte gels.
Xuanhe Zhao has just posted a new paper from my group on dielectric gels.
A thread of discusion on Maxwell stress and electrostriction has just started.

Modélisation numérique du procédé de forgeage radial

Olivier Pantalé — Wed, 29 Aug 2007 09:49:08 -0400

Dr Olivier PANTALE

Assistant Professor
G2TR - LGP - ENI Tarbes

47, avenue Azereix - BP 1629
F-65016 Tarbes Cedex

Papers on nonlinear deformation of dielectric elastomers

Zhigang Suo — Wed, 30 May 2007 10:25:51 -0400

Dear Aman: In addition to this paper, members of my group have posted the following papers on iMechanica:

These papers have shown how the basic formulation may be used to analyze some of the intriguing experimental observations reported recently. Please do let us know if you have comments.

Accepted: A nonlinear field theory of deformable dielectrics

Aman Haque — Wed, 30 May 2007 09:54:12 -0400

Dr. Suo,

Reading the manuscript, I couldnt curb my curiosity..can I have copies of the papers Zhao et al. and Zhou et al.?

Many thanks

Aman

charge control

Bill Greene — Sun, 27 May 2007 16:33:33 -0400

Hi Jinxiong,

I actually think that the charge control strategy described by Wei is quite a good approach to this problem. In fact, after reflecting on this more, I think my suggestions about the possible need for an arc-length algorithm may have been a bit hasty.

I think it is more likely that the convergence difficulties are simply due to numerical ill-conditioning near the bifurcation point. Have you done any numerical experiments increasing the size of the imperfection to see to what extent that improves convergence?

Bill

Back to voltage control if we use arc length method?

Jinxiong Zhou — Sat, 26 May 2007 21:11:52 -0400

Hi Bill,

Another interesting question. Previously Wei mentioned that we must introduce some triky treatment of potential on the boundary where charge was applied and controled as load increments. I am still thinking that the voltage-charge curve should be continous even in the jump area. I mean following the jump line the voltage decreases continuously with charge almost fixed. If the jump is big, small variation of charge may correspond to big change of voltage. This is the reason why New-Raphson method fails to converge.

However, the votage can increase or decrease automatically if we use arc length method, because we actually use the length of arc as parameters. If that is the case, maybe we can back to potential loading and we can avoid any additional treatment of potential!

Can I have your comments on this point?

Jinxiong

Arc-length method for mixed DOFs?

Jinxiong Zhou — Sat, 26 May 2007 20:52:03 -0400

Hi Bill,

Thanks for your valuable comments! You suggested a feasible way to obtain the critical load value, and I will try it later.

I know ABAQUS uses the arc-length method or you previously mentioned the Riks mehtod for postbuckling analysis. And the reported researches show that this method is effective for static postbuckling analysis. Of course, most of these reported researches involve only only type of DOFs, the displacements. What interests me is you mentioned that the method is disable for mixed DOFs. I have no experience of using ABAQUS to analyze postbuckling analysis and I am not clear what behind this invalidation of arc-length method.

Is there anybody who use Riks method in ABAQUS to analyze stability problems with different types of degree-of-freedoms (displacements plus whatever else)? Can anybody give us some insightful suggestions?

Jinxiong

Equilibrium paths and arc length methods

Bill Greene — Sat, 26 May 2007 10:56:36 -0400

Jinxiong,

Thanks very much for your detailed reply.

From your paper and this discussion, it seems clear that there are two equilibrium paths. Both of these paths should be reasonably continuous. For the perfect structure, I'm assuming that the secondary path intersects the primary path at least once at a bifurcation point-- possibly on section B of the curve in your figure 6. When you introduce an imperfection, you convert this bifurcation problem into one that includes an equlibrium path with a limit point. This path is also continuous (i.e. no "jumps"). However, in the vicinity of the bifurcation point (singular point), this imperfect structure is so ill-conditioned numerically that it is extremely difficult to obtain a converged solution. I think that is why you are seeing what appears to be a jump; its just very difficult to obtain converged solutions for the intermediate points on the equilibrium path.

You may be able to obtain some additional insights by studying the behavior of the perfect structure. Specifically, if you trace the primary equilibrium path and observe sign changes in the determinant of the jacobian, it should be possible to estimate the location(s) of the bifurcation points.

I agree with your comments about arc-length methods. You'll notice I made a couple of comments about these in another reply on this thread. As I mentioned below, I have never tried these when the arc length contains two different types of degrees-of-freedom. I suspect there may be some sticky issues of how to scale the different DOF types in computing the arc length. One reason why I say this is that ABAQUS contains an arc length algorithm that is very effective in problems with only displacement DOFs. However, this algorithm is disabled when the problem contains different DOF types (e.g. displacement and voltage).

Bill

Thank you for your insightful questions

Jinxiong Zhou — Fri, 25 May 2007 21:02:21 -0400

Hi, Bill and Wei,

Thank you for your insightful questions and comments. Today I am busying myself writing and debugging another code so I missed your helpful discussions. So I will read carefully your questions and comments and try my best to answer your questions.

First, I will talk about something about tracking the unstable path. From figure 6 we can see that if we use voltage as loading increments, when the voltage reaches the peak value, the voltage must decrease. At this point, even you decrease the increments of voltage factitiously, the solution will go back and follow the previous solutions if no special treatment is introduced. Alternatively, we use charge as loading increments because charge is monotonically increasing. This technique works only for very small jump when the solution at the peak is near to the solution on the horizontal line. If the jump is big, i.e., the difference between the latest known homogeneous solution and the first inhomogeneous solution is big, the last homogeneous solution is not a good guess of the inhomogeous solution and in this case the New-Raphson iteration method can not converge. To be honest, this remains an issue for our programm and so far our programm can not fallow the unstable path when big jump occurs.

Actually, our problem is the same as the post-buckling problem in strucutral mechanics or the computation of birfurcation in nonlinear dynamics. As far as I know, the problem is solved and there are some approaches are available. One of the method I used many years ago is the parametric continuation method, e.g. the arc-length continuation method. The key point of the continuation method is to convert increments of a physical parameter, e.g., charge or potential in our method, to a geometric parameter-the arc length. The increments of physical parameters may increase or decrease during the whole solution diagram, while the arc length of the solution path is definitely increasing monotically. This point is manifested if you just use a one degree freedom system. So you have an additional variable and you need an additional equation added to the New-Raphson method. This additional equation is called arc-length equation. Therefore, even the original Jacobin matrix is singular, the matrix of the augmented system is no longer singular. In this way the method can follow any unstable path and can obtain the whole solution diagram. This method can be implemented at the cost of more complicated algorithm and programming. We have not implemented continuation method yet, but hopefully this can be done in the follow up works.

I also agree with you that for engineering applications only the first stable solution is important and in most cases the actuators will breakdown before the second stable state is observed. I am not clear whether this big jump is useful for some potential applications.

Hope that my reply can answer your first question. I will read your other questions and Wei's brilliant answers carefully and hopefully can involve into this hot discussion.

Jinxiong

numerical and analytical

Wei Hong — Fri, 25 May 2007 15:11:33 -0400

Yes, you are right, Bill.

To obtain the unstable result of homogeneous deformation, one also needs to contrain the deformation so that it is uniformly compressed. We have given it a little bit thought before making the plot, such as combine the DOFs of the vertical displacement of the top-nodes. However, different from the charge control, the deformation constraint is artificial, and the only result it can get is to have the calculation follow the analytical curve all the way, which is sort of borring anyway.

For the plots in Fig. 4, actually there is a little discrepency. The block does show some inhomogeneous deformation on the corner. But as the region is small, it contributes little to the overal measurements, say the voltage and total charge.

For the interesting transition points near C or D, although thinking as an ideal thermodynamic system, the preferable line would be horizontal line CD intersects with the primary pathes. However, to have the coexistent states, as the system does in between C and D, it has to pay extra energy to form an transition band between the two phases. Which can be seen on fig 7 C and D, the transition band is having a width of about the thickness of the film. The extra cost in forming the band prevent the film to jump from the primary curve to CD directly. Instead, it overshoots a little before making the transition. The amount of the overshoot depends on the formation energy of the transition band.

You can also think in this way: The CD segment corresponds to the coexistent states of the two thickness, so you do need both stable states. However, the transition band from thick to thin states has finite width. If the necking region is even smaller than that width (e.g. the state a little bit prior to Fig. 7C), the coexistent of states simply can not happen. Comparing to the state of having a partially necking in the middle, the system would rather prefer a uniform thick state. That's why we always have an overshoot before making the turn.

You guessing might also be right that where to make the turn also depends on the algorithm itself, as the discrepency between loading and unloading. There is no proof that N-R method would prefer an energetically more favorable state, and the convergence region of N-R method is initialvalue dependent and somewhat arbitrary or fractal.

Wei

Numerical solutions for alternate equilibrium paths

Bill Greene — Fri, 25 May 2007 13:15:56 -0400

Wei,

Thanks for the
detailed and thoughtful reply. Do you mind a couple more questions and
comments? One is based on a comment from your last reply that I thought I had
understood from the paper.

I've included parts of your reply below (italics font) so I
can refer to them directly.

For the neo-Hookean material, the curve does not go up
again, which means the coexistent states include one thin state of thickness 0.
Or, simply there are no coexistent states. Once you have a separation,
the system goes all the way to 100% thin state. Such a change is too
significant for Newton-Raphson method to handle. So the calculation fails
instead of the collapse of the film.

In my previous reply, my
first thought was that your charge-control algorithm would follow the
equilibrium path almost until the block reached zero-thickness. However, after
reflecting on this more, I think it may be necessary to constrain also the displacement
increments in addition to the charge. This would be a Riks (or Crisfield) type
algorithm that you may be familiar with. I must admit, though, that I have no
experience with these algorithms for mixed degree-of-freedom types (e.g.
displacement and potential) in the solution vector.

For Boyce-type of materials, if we don't introduce any imperfection,
and use a square
block and a uniform mesh, the calculation can go into part of the unstable
state. This does not mean our calculation preduces artificial results.
It's just the curve is in the unstable regime for an "homogeneous
deformation state". But if you look at the deformation pattern, it
is not uniform any more. The block is more or less tilted. This means the
cost of the transition region is so expensive that the system can not afford to
have coexistent states. Instead, it has mostly a thick region and some
little transition at the edge. Such a inhomogeneous deformation state is
actually stable.

However, if you keep on loading, the transition region
will grow and graduate eat up the whole block. That's when the curve will
deviate from the analytical solution.

I thought this is exactly the problem you show results for
in Figure 4 where the numerical solutions lie right on the analytical curves? I
had assumed these numerical solutions had uniform displacements because they
agree so well with the analytical solutions. What am I misunderstanding?

For the difference in loading and unloading, I
guess it is because the imperfection we introduced does not have a same
contribution to the transitions. (One being local thinning and the other being
local popup) It is not a big difference and the imperfection we introduced is
artificial, I don't think that's a physical phenomenon that worth to be
studied.

One of the main reasons for my questions is to try and
understand the little C-D segment of the equilibrium path in your Figure 6. I'm
guessing that it actually intersects the B and E sections of the primary
equilibrium path but that it is numerically very difficult to get converged
solutions closer than you have to the primary path? Does this speculation agree
with your numerical experiences? Here is another guess about the "loading" and
"unloading" behavior. I'm guessing that the Newton-Raphson algorithm is just
(somewhat) arbitrarily converging to one equilibrium path or the other based on
whether you are "loading" or "unloading" in your solution algorithm. If so, this
"choice" doesn't really imply anything about how the solid actually behaves.

Thanks again.

Bill