User login


You are here

Time-History Analysis of Simple Structures under Shock

Dear all mechanicians,

I modeled a disk with outer diameter of 48 mm and inner diameter of 24 mm. The material is defined as isotropic; in this case, it's glass-ceramics (E = 90GPa, v = 0.24). The inner diameter is constrained to transverse directions, while it is subjected to half-sine wave impulse of 500G (G = 9.81 m/s^2) within 2 ms. The disk moved up under this loading, and I extracted the vertical displacement at the tip with respect to time (time-history analysis). The time step was 1e-6 second, and the duration was 4 ms. The disk was made of 8-node hexahedron element, with one element in the thickness. Three formulations were tested: (1) reduced integration [RI], (2) selectively-reduced integration [SRI] and (3) full integration [FI].

After conducting convergence study (starting from a very coarse mesh model to a very fine model), the conclusion is this: the model with FI is converged faster. It means that when I was comparing maximum displacement, coarse-mesh model will give similar value to the fine-mesh model. Whilst for SRI and RI, coarse mesh will give lower value, and after mesh refinement, it will give similar value to the FI results.

I did simulation with explicit solver (LSDYNA). The conclusion is the same.

What could be the cause of the slow convergency of RI and SRI? Can anyone refer me to the literatures on this dynamic problem?

Thank you!


I can't understand how your problem involves a shock. Also, I don't see why you want to use reduced integration. Is it because you want to reduce the compute time? What sort of stabilization do you use for your reduced integration calculations?

It would help if you could upload a figure showing your problem set up and your results.


Dear Biswajit,

This is to simulate the sudden stop of the hard disk drive.

 Dear Yudhanto,

 I conducted some work on Hard disk drives with seagate. I guess one layer of solid element may be insufficient to model the vibration behavior of the disk. The results converge faster with shell elements.


Luo Jun


I think I have a better idea about the problem though I still don't understand some of the details of the setup.

There are few details that were left out in the original question:

1) The thickness of the disk.  This is crucial if you wish to determine why you get the rate of convergence that you do.

2) The density of the material.  This is needed if you want to find the rate at which the stress wave propagates and whether the shock speed is faster than the speed of sound in the material.

3) The mass of the disc which is needed if the 500g applied acceleration is to be converted into a force. 

4) The direction of the applied force - is it normal to the disk or radial.  Since you specify it in units of g-force it appears to be of the centrifugal type. 

5) The actual type of element that you have used.  The 8-noded brick that you mention (with selective reduced integration) is not a solid element but actually a "continuum-based" shell element in Belytschko's language.  See for example
Nonlinear Finite Elements for Continua and Structures by T. Belytschko, W. K. Liu, and B. Moran, John Wiley and Sons, 2000. 

Let us assume that the thickness of the disk is 1 mm and that the density of the material is 5180 kg/m^3 (this is true for certain glass-ceramics).  Then the shear wave speed in the material is 2650 m/s and the bulk wave speed is 3340 m/s.  If the impluse is a true shock then it should travel across the thickness and along the plane of the disk at speeds greater than (say) 3340 m/s.

If you want to resolve the shock, the time increment should be appropriately determined.   Your convergence results will strongly depend on what time step you choose.  Also, depending on the direction of the impulse, you might have to use standard hex elements.  In that case, you will have to use a number of elements across the thickness to prevent locking.

From what I've read so far in this thread, it appears that you are indeed using shell elements (even though they have eight nodes).   That is why your results converge  (do they really converge to the correct solution?) to a solution faster with the fully integrated element.  However, I can't recall what selective reduced integration means in the context of continuum-based elements.  The book by Belytschko explains the details of some of the formulations that you have looked at in a concise manner.  You'll probably find the exact answer you are looking for in that book.


I appreciate the "refinement" on the questions, Biswajit; since it'll lead to better understanding about the problem.

(1) The thickness of the disk = 0.5 mm

(2) Density of glass-ceramics = 2530 kg/m^3

(3) Applied acceleration was defined as displacement function in ANSYS

(4) The direction of applied acceleration is normal to the disk; and it's applied to the inner diameter of the disk.

(5) It is 8-noded hexahedron element. In ANSYS, it's called "SOLID45". Shear locking can be "cured" by selecting extra displacement shape (formulation based on Wilson, Taylor, Doherty and Ghaboussi, 1973). The time-step was 1e-6 seconds, which I believe, it's small enough for this particular problem. When I reduced the time-step, e.g. 1e-7 seconds, the result was the same as the one with 1e-6 seconds.    

In ANSYS/LSDYNA, when we select 8-noded hexahedral element (it's called SOLID164), two element formulations are available; they are constant stress (reduced integration) which is one-point integration with viscous hourglass control, and selectively-reduced integration with enhanced strain. The latter is formulated based on Puso's paper "A highly efficient enhanced assumed strain physically stabilized hexahedral element", IJ Num. Meth. Engg., 2000.

Btw, here, I am concerned with the FE software users who are not really aware with the convergence problem that they may encounter during the shock simulation of HDD when they use lower-order hex element. 

Thank you again for the references.  




Thanks for the further clarifications in your comment.

In my previous response to your post I had written "The 8-noded brick that you mention (with selective reduced integration) is not a solid element but actually a `continuum-based' shell element". That statement is not correct. Continuum-based shell elements are two-dimensional elements and not "bricks".

Your problem brings up a few interesting issues.

  • How appropriate it is to use an elastic material to model these disks? When you apply a pulse to the disk does it vibrate for ever in vacuum? I doubt that. In that case, how long does it take for the elastic waves to be damped out? If the damping time is small (say a second or so), some dissipation has to be incorporated into the material model. Could an expert on dissipative elasticity point us to some literature on that subject?
    Since you don't include dissipation in your simulations, how do you determine when the simulation has converged? Do you wait for the motion to reach a steady state? Or do you include some artificial viscosity and damp out the motion?
  • Clearly the air around the disk also damps the motion. If believe this is what you mean by the ``air-bearing''. How strong is this effect compared to the dissipation within the material itself? How would you include this damping effect in a numerical simulation?
  • The rate at which the elastic wave propagates will be affected by the discretization of the mesh. See for example the schematic shown in the figure below.

    The location of the stress front as seen by the mesh at the first timestep is shown by the green dot, the blue dots show the second timestep, and the red dots show the third timestep. Clearly the rate at which the wave propagates depends on the mesh discretization and the chosen timestep. Also reflections from the bottom of the disk will affect the nodes at the top of the disk much earlier if you use only one element across the thickness.

    Hence the issue here is not really whether the element can support bending (though that is part of the issue) but how well you want to represent the stress waves. The accuracy of you results will depend on how finely you discretize the disk in the thickness direction.

Dear Biswajit,

In terms of material, there's a simplification in solving HDD shock problem. The material is defined as elastic only. Hence, the elastic wave dan damping may not be captured during simulation; the oscillation occurs after the shock has been terminated.

Air bearing is built below the slider when the slider flies over the spinning disk. To represent the air bearing people usually use spring-damper element or apply force - moment to several points of the slider.   

Since I have not published my work yet, you may refer to following paper to gain physical problem:

QH Zeng and DB Bogy "Numerical simulation of shock response of disk-suspension-slider air systems in hard disk drive", Mycrosystem Technologies (8), 2002, pp. 289 - 296, Springer Verlag.  

Similar research were also done by Eric Jayson and Frank Talke of Univ California San Diego, David Shu Dongwei - Luo Jun et al of Nanyang Technological University (Singapore), etc.


Following your advice I went and did some more reading of the literature on disk drives. My curiosity on this subject is due to the combination of shocks and fluid structure interaction that I deal with in my research (in an entirely different area). Please e-mail me the paper by Zeng and Bogy if you can.

I assume that you are primarily interested in small disks meant for portable machinery. In that case, accelerations of the order that you mention should lead to significant deformations - particularly because the disks are thin. Do you have some numbers on what deformations you get?

When you say "... the oscillation occurs after the shock has been terminated" what exactly do you mean? Clearly the shock pulse will continue to the outer radius, reflect back, interfere with other pulses, and so on. But will never damp out if you use a perfectly elastic material. However, the response of the disk will be significantly different if the material is dissipative. Are you saying that these effects are not of interest as far as your particular problem is concerned?

I don't really understand what it is that you are trying to find out. Is it whether the disk will break under the impulse? If that is the case, what failure criteria do you use. Or are you trying to find out whether the disk will impact other components of the drive system?

I went and looked up some material on air bearings on the web. From what I gather, these bearings are in the spindle subassembly at the center of the disk. However, the air bearing you are talking about is below the slider. I am confused here and an explanation of the situation will be really helpful.


The paper has been sent to your email.

Allow me try to answer you the best I can. 

What I'm trying to simulate here is very simple actually:

Disk with outer diameter of 48 mm, inner diameter of 12 mm, thickness of 0.5 mm is modeled in ANSYS using 8-node hex elements. The material property defined is glass-epoxy with E = 90 GPa, ν = 0.24, ρ = 2.53 x 10-6 kg/mm3. The nodes at the inner diameter of the disk is subjected to half-sine acceleration impulse with the peak of 500G (G = 9.81 m/s2) and the duration is 2 milliseconds. The direction of applied shock is vertical or upward. The inner diameter is constrained to move in translation directions, but it can move vertically. The applied shock is terminated at 2 ms, while the disk response (deformation) is of interest. The displacement of disk outer diameter (disk-tip) is measured with respect to time (time-history analysis). From here we get the shock response of the disk in terms of displacement. For different element formulation, ANSYS will give different value of displacement. As a benchmark, maximum displacement is 0.1 mm. Damping ratio is not yet included in this simulation. The objective is simply to find tip displacement of the disk with respect to time. The failure of the disk is not yet a consideration, and it is expected that impulse 500G will result in disk deformation only, not failure (break) of the disk.

Second problem is disk with suspension. Below the suspension, between the disk and the suspension, there's a tiny little block called "slider". The slider is flying on the top of the disk when the HDD is in operating condition (meaning: read/write operation is ongoing). The reason why it's flying is because of the air-bearing under the slider. It's built up by the unique contour of air-bearing surface of slider when the air (as a result of spinning disk) is passing it. People may model this air-bearing with spring-damper element. The stiffness (force and moment as well) can be obtained by solving Reynold equations.

The "bearing" that you mentioned could be spindle bearing which holds the disk. Another "bearing" is pivot bearing, where the suspension is attached to.

Zeng and Bogy paper might give you some idea how people use finite element method to solve this problem.    

Dear Biswajit,

The half-sine wave impulse was to replace impact response of disk drive. It's usually defined as half-sine wave, otherwise described differently by company's specification.

People who resort this problem with LS-DYNA or ANSYS/LSDYNA will be normally directed to the default element technology, which is reduced integration (constant stress) or RI. Another option is to select selectively-reduced integration or SRI. Some strategies to reduce hourglassing effect are introduced for both options. For RI, LSDYNA provides their own formulation plus four variants of Flanagan-Belytschko. For SRI, Puso hourglass control is introduced in LSDYNA. I used all of the formulations, and the results using Flanagan-Belytschko (stiffness form) converged faster than Puso or other variants of Flanagan-Belytschko.

While in ANSYS, we can have three options which include RI, SRI and FI. For FI, we have freedom to choose "include extra displacement shape" or "exclude extra displacement shape". Excluding extra displacement shape will generally result in more reasonable outcome.

Dear Luo,

Yes, I've read several papers published by you and coworkers (at NTU) on shock of HDD. My question is what was the element formulation did you use when you model the disk? Since you're using ANSYS-LSDYNA, you can have two options: RI or SRI. (Referring to a paper "Study of shock response of the HDD with ANSYS-LSDYNA", J. Magn. & Magnetic Mat'ls, 2006). 

Indeed, shell element will give faster convergency when we compare the displacement of disk-tip. However, solid element is used whenever two contact surfaces are to be defined for lower and upper suspensions.  

Thanks, Biswajit and Luo.





Dear Yudhanto,

Thanks for your interest in the paper. In that study, I adopted the default element type in ANSYS-DYNA: The Belytschko-Tsay (default) element formulation with hourglass control. This element uses one point quadrature. It's always a problem as we lack comparable shock test data. Fully integration formulation may also introduce some numerical problems such as stiffening. BTW, you can still use shell element to simulate the disk when there are two suspensions. This is what I did in that paper. The contact algorithm in DYNA considers the shell thickness.



Department of Mechanics

Huazhong University of Science and Technology, Wuhan, China.

Thanks, Jun. I may try to use shell element for two-suspension case.

As for fully-integrated 8-node hexahedron, ANSYS provides extra-displacement shape as I have explained to Biswajit. Shear locking can be avoided using this technique.

What contact property did you use in LSDYNA? I understand that you defined "automatic surface-to-surface". It's based on penalty method, isn't it? Since you were simulating non-operating shock, dynamics of head slap may be governed by predefined contact property. In other words, different contact property may give different head-slap dynamics.   

Dear Yudhanto,

Correct! Great to learn DSI is doing the same thing. Have you gone to op-shock?  Prof. Bogy did very good job in  this area, but I am still not so sure some of the details in their op-shock modelling. The FH of the slider is on the order of several nanometers and the air bearing dynamics are coupled with a relative larger components like suspensions. How can the solution of the suspension part be so accurate? Do you have any ideas on this point?

I left NTU in Dec of 2005. Quite hard for me to continue this work without industry support.



 Department of Mechanics

Huazhong University of Science and Technology, Wuhan, China.

Dear Jun

I have some preliminary simulation results on op-shock, but I need to "cook" them up further to obtain more reasonable data.

You're right, Jun: Bogy group (UC Berkeley) are leading researchers in this matter. Others are Talke group (UC San Diego), Shu Dongwei (I believe you remember this name - NTU), Yap Fook Fah group (NTU), DSI (S'pore), etc. Each of them comes up with their own focuses. So it's quite exciting to see the growing enrichment in this field.

Answering your question: as you may know, flying height modulation during- and post-shock is generally governed by air-bearing stiffness (vertical, roll and pitch), dimple-slider contact, stiffness of suspension, shock amplitude and duration, and other factors. Nonlinearity of the air-bearing is believed to give significant effect when HDD is subjected to high-G shock. In numerical viewpoint, I'd say that preprocessing part has very significant role; including the selection of appropriate element formulation for disk-suspension system. Accuracy must come first; then, we need to make adjustment on solving computational time.

Currently, the gap is in the experimental validation of numerical results. Some people may think that it's not really useful to compare numerical solution and experimental result since both conditions are separated by initial assumptions and unknown factors. Nevertheless, I believe it's still worthwhile to validate numerical results with experiments since we want (at least) to see the trend of the result; not the exact value.

What do you think, Jun? 

Indeed, you are again correct. Without any industrial collaboration, the process of getting it done is somehow unsmooth. The research direction may be deflected somewhere else as well.  

Dear Yudhanto,

         Good!  I am looking forward to the appearance of your work. We also did some work on op-shock modelling in Prof. Shu group. But the suspension dynamics and air bering dynamics were not fully coupled then.  I guess there should be some more progress on this topic in Prof. Shu's group till now. Prof. Shu (my postdoctoral advisor) group is also doing some shock tests. You know, quite hard topic. Not so easy to capture the shock behavior at the head-disk interface. I guess that's why we see so many numerical simulations but seldomly find corresponding experimental work. 

You said it. The experimental part is essential for several reasons. Firstly, to play with real things help us to understand their function mechanism and thus help we set up meaningful numerical models. Secondly, numerical simulations need experimental verification. Otherwise we may get plenty of beautiful useless results. And vice versa.

Department of Mechanics

Huazhong University of Science and Technology, Wuhan, China.

Henry Tan's picture

Dear Arief and Biswajit,

Your discussions on shock is very interesting to me.

Therefore, I linked this blog to my topic on explosion science and engineering (

Welcome to post something on my blog. It is too silent there.

Subscribe to Comments for "Time-History Analysis of Simple Structures under Shock"

Recent comments

More comments


Subscribe to Syndicate