A small version of the system, called FEAPpv, is available free of any charge. See projects. Zienkiewicz and R. Taylor and published by Elsevier, Oxford, PDF files require Acrobat Reader to print. This is a site to post questions, get answers, and interact with other users. We encourage users to help other users.
It is just in its beginning stage but already contains some useful information for new and potential users. Please have a look and provide feedback and content.
The file is in 'zip' form. The stress is calculated using the strain rate and the material law provided by the user. The volumetric strain using the small strain formulation is independent of density.
In the crushing direction, honeycomb has no Poisson's effect and stress integration over the initial surface is acceptable. The effect on strain is small during elastic deformation and can be corrected in the plastic phase by using a modified engineering stress-engineering strain material curve. For materials like crushable foam, with a small Poisson's ratio, this formulation can be applied successfully in certain situations. However, this formulation has to be used very carefully for other materials.
Shell elements have fewer limitations than solid elements. For crash applications, the main shell deformation is bending. The small strain formulation has no effect on the bending description if membrane deformation is small. The small strain formulation can be applied to some elements for which the time step is reaching a user specified value. If the critical time step is small, compared to the initial one, this formulation gives acceptable results and is more accurate than removing the deformed elements.
By computing the derivative of shape functions at each cycle, large displacement formulation is obtained. The large strain formulation results from incremental strain computation. Stresses and strains are therefore true stresses and true strains.
Stress rate is a function of element average rigid body rotation and of strain rate. The time step is computed at each cycle. Large element deformation can give a large time step decrease. For overly large deformations a negative volume can be reached and it then becomes impossible to invert the Jacobian matrix and to integrate the stresses over the volume. Summation over EQ. In the case of lower indices, summation is over the number of space dimensions. For upper case indices, summation is over the number of nodes.
The nodes in the sum depend on the type of entity considered. When the volume is considered, the summation is over all the nodes in the domain. When an element is considered, the sum is over the nodes of the element. All the dependency in the finite element approximation of the motion is taken into account in the values of the nodal variables.
From EQ. The test functions are next substituted into the principle of virtual power EQ. The expression applies to both the complete mesh or to any element. It is pointed out that derivatives are taken with respect to spatial coordinates and that integration is taken over the current deformed configuration. The equations of motion for rotational degrees of freedom are complicated if written in the global reference frame.
They are much simpler if written for each node in the principal reference frame attached to the node. The resulting equations are the standard Euler equations. The equation of motion for rotational degrees of freedom is thus very similar to that for translational degrees of freedom.
Equation 3. It is shown hereafter that expressing the shape functions in terms of intrinsic coordinates is equivalent Finite elements are usually developed with shape functions expressed in terms of an intrinsic coordinates system to using material coordinates. They simply relate coordinates in the physical world to the EQ. Writing EQ. External forces and the mass matrix can similarly be integrated over the domain in the intrinsic coordinate system.
These spatial derivatives are obtained by implicit differentiation. Considering the velocity gradient e. Usually, it is not possible to have closed form expression of the Jacobian matrix. As a result the inversion will be performed numerically and numerical quadrature will be necessary for the evaluation of nodal forces. For full integration, the number of integration points is sufficient for the exact integration of the virtual work expression.
The full integration scheme is often used in programs for static or dynamic problems with implicit time integration. It presents no problem for stability, but sometimes involves "locking" and the computation is often expensive. Reduced integration can also be used. In this case, the number of integration points is sufficient for the exact integration of the contributions of the strain field that are one order less than the order of the shape functions.
The incomplete higher order contributions to the strain field present in these elements are not integrated. The reduced integration scheme, especially with one-point quadrature is widely used in programs with explicit time integration to compute the force vectors.
In two dimensions, a one point integration scheme will be almost four times less expensive than a four point integration scheme. The savings are even greater in three dimensions. The use of one integration point is recommended to save CPU time, but also to avoid "locking" problems, e.
Shear locking is related to bending behavior. In the stress analysis of relatively thin members subjected to bending, the strain variation through the thickness must be at least linear, so constant strain first order elements are not well suited to represent this variation, leading to shear locking. Fully integrated first-order isoparametric elements tetrahedron also suffer from shear locking in the geometries where they cannot provide the pure bending solution because they must shear at the numerical integration points to represent the bending kinematic behavior.
This shearing then locks the element, i. On the other hand, most fully integrated solid elements are unsuitable for the analysis of approximately incompressible material behavior volume locking. The reason for this is that the material behavior forces the material to deform approximately without volume changes.
Fully integrated solid elements, and in particular low-order elements do not allow such deformations. This is another reason for using selectively reduced integration. Reduced integration is used for volume strain and full integration is used for the deviatoric strains. However, as mentioned above, the disadvantage of reduced integration is that the element can admit deformation modes that are not causing stresses at the integration points.
These zero-energy modes make the element rank- deficient and cause a phenomenon called hour-glassing: the zero-energy modes start propagating through the mesh, leading to inaccurate solutions. This problem is particularly severe in first-order quadrilaterals and hexahedra. To prevent these excessive deformations, a small artificial stiffness or viscosity associated with the zero-energy deformation modes is added, leading in EQS 3.
For each time step in a particular analysis, the algorithm used to compute results is: 1. For the displacement, velocity and acceleration at a particular time step, the external force vector is constructed and applied. A loop over elements is performed, in which the internal and hourglass forces are computed, along with the size of the next time step. The procedure for this loop is: 2a. The stress rate is calculated: EQ. Cauchy stresses are computed using explicit time integration: EQ.
The internal and hourglass force vectors are computed. The next time step size is computed, using element or nodal time step methods see chapter 4. After the internal and hourglass forces are calculated for each element, the algorithm proceeds by computing the contact forces between any interfaces. Finally, time integration of velocity and displacement is performed using the new value.
However, in practical finite element analysis, a few effective methods are used. The procedures are generally divided into two methods of solution: direct integration method and mode superposition.
Although the two techniques may at first look like to be quite different, in fact they are closely related, and the choice for one method or the other is determined only by their numerical effectiveness.
In direct integration the equations of motion are directly integrated using a numerical step-by-step procedure. In this method no transformation of the equations into another basis is carried out. As the solution is obtained by equation written at discrete time points includes the effect of inertia and damping forces. The variation of a step-by-step procedure, the diverse system nonlinearities as geometric, material, contact and large deformation nonlinearity are taken into account in a natural way even if the resolution in each step remains linear.
The mode superposition method generally consists of transforming the equilibrium equation into the generalized displacement modes. An eigen value problem is resolved. The eigen vectors are the free vibration mode shapes of the finite element assemblage. The superposition of the response of each eigen vector leads to the global response. As the method is based on the superposition rule, the linear response of dynamically loaded of the structure is generally developed.
In the following, we describe first the resolution procedure in direct integration method when using an explicit time discretization scheme.
Then, the procedures of modal analysis are briefly presented. The implicit method will be detailed in chapter A few commonly used integration methods exist in the literature [54].
This section deals with time integration of accelerations, velocities and displacements. The general algorithm for computing accelerations, velocities and displacements is given.
Stability and time step aspects are then discussed. It can be shown that it is always unstable. An integration scheme is stable if a critical time step exists so that, for a value of the time step lower or equal to this critical value, a finite perturbation at a given time does not lead to a growing modification at future time steps.
It can be shown that it is conditionally stable. This integration scheme is the unconditionally stable algorithm of maximum accuracy. The time step h may be variable from one cycle to another. It is recalculated after internal forces have been computed. It is pointed out that the solution of the linear system to compute accelerations is immediate if the mass matrix is diagonal.
Figure 4. The problem is similar in landing of aircrafts. It can be studied by an analytical approach where the dropping body is modeled by a simple mass-spring system Figure 4.
If h is the dropping height, m and k the mass of the body and the stiffness representing the contact between the body and the ground, the equation of the motion can be represented by a simple one d. It is compared to the analytical solution given by EQ. The difference between the two results shows the time discretization error. A numerical procedure is stable if small perturbations of initial data result in small changes in the numerical solution.
It is worthwhile to comment the difference between physical stability and numerical stability. Numerical instabilities arise from the discretization of the governing equations of the system, whereas physical instabilities are instabilities in the solutions of the governing equations independent of the numerical discretization.
Usually numerical stability is only examined for physically stable cases. For this reason in the simulation of the physically unstable processes, it is not guaranteed to track accurately the numerical instabilities. Numerical stability of a physically unstable process cannot be examined by the definition given in above.
We establish the numerical stability criteria on the physically stable system and suppose that any stable algorithm for a stable system remains stable on an unstable system [36]. On the other hand, the numerical stability of time integrators discussed in the literature concerns generally linear systems and extrapolated to nonlinear cases by examining linearized models of nonlinear systems.
The philosophy is the following: if a numerical method is unstable for a linear system, it will be certainly unstable for nonlinear systems as linear cases are subsets of the nonlinear cases. Therefore, the stability of numerical procedures for linear systems provides a useful guide to explore their behavior in a general nonlinear case. To study the stability of the central difference time integration scheme, we establish the necessary conditions to ensure that the solution of equations are not amplified artificially during the step-by-step procedure.
Stability also means that the errors due to round-off in the computer, do not grow in the integration. It is assured if the time step is small enough to integrate accurately the response in the highest frequency component. The equations to relate displacements, velocities and accelerations in a discrete In direct integration method, at time tn the solutions for the prior steps are known and the solution for the time time scale using the central difference time integration algorithm are given in section 4.
A spectral analysis of this matrix highlights the stability of the [ ] integration scheme: [ ] The numerical algorithm is stable if and only if the radius spectral of A is less than unity. In the other words when the module of all eigen values of A are smaller than unity the numerical stability is ensured.
Three cases are to be considered: 1. The corresponding part of the boundary of the stability domain is the segment analytically defined by 2. Their modulus is equal to 1. To precisely define the stability domain, we must also and the boundary of the domain.
Rewriting the central difference time integration equations from EQ. Stability is given by the set of conditions from EQ. It is the same for the right inequality of the first expression.
We may yet remark that damping has changed the strict inequality into a large inequality, preventing from weak instability due to a double eigen value of modulus unity. It is important to note that the relation EQ.
This will be studied in the next section. Except for the case of full diagonal, contrary to [M ]. We therefore often compute the viscous forces using the velocities at the preceding mid-step, which are explicit. Stability is again given by the set of 2 2m m conditions EQ.
As described in the preceding sections, the computation of the viscous forces by using velocities at time steps [ ] leads to obtain a non-diagonal matrix [C] which should be inverted in the resolution procedure.
To avoid the Rayleigh equation into EQ. Stability is obtained as before by means 2 2 of the set of conditions from EQ. It is advantageous to separate the two contributions, restrictions of the time step then becoming lighter. This value is less than the critical time step obtained by: 4. Using now the values higher than or equal to the critical time step, the solution diverges towards the infinity as shown in Figure 4.
In case of m k k m dropping body example, the mass-spring system can be compared by analogy to a two-node mass-spring system where the global stiffness is twice softer. The critical time step is then computed using the nodal time step of the entire element refer to the following sections for more details on the computation of nodal time step.
This reduces computational time considerably as no matrix inversion is necessary to compute accelerations. The stability condition is given in the last sections. The time step restriction given by EQ. These changes in the material and the geometry influence the value of and in this way the critical value of the time step. The above point can be easily pointed out by using a nonlinear spring with increasing stiffness in example 4. It can be shown that the critical time step decreases when the spring becomes stiffer.
Therefore, if a constant time step close to the initial critical value is considered, a significant solution error is accumulated over steps when the explicit central difference method is used.
Another consideration in the time integration stability concerns the type of problem which is analyzed. For example in the analysis of wave propagation, a large number of frequencies are excited in the system. That is not always the case of structural dynamic problems. In wave shock propagation problem, the time step must be small enough in order to excite all frequencies in the finite element mesh.
This requires so short time step that the shock wave does not miss any node when traveling through the mesh. The condition EQ. If a uniform linear-displacement bar element is considered, see Figure 4. Combining EQ. This is shown for the first time by Courant et al.
In spite of their works limited to simple cases, the same procedure can be used for different kinds of finite elements. The characteristic lengths of the elements are found and EQ. Regarding to the type shape of element, the expression of characteristic length is different. The step used for time integration or moving forward in time can be calculated using two different methods.
The method used depends on the type of simulation being performed. The time step used by the solver is the largest possible time step, as determined by the Courant condition that will maintain stability.
If the default large strain formulation is used, the time step is computed at each cycle. If the deformation is too large, negative volumes can result, which make it impossible to invert the Jacobian matrix and to integrate the stress over the volume. This is either the beginning of the analysis or the time at which the small strain formulation is initiated. If the sound speed is constant, the time step thus becomes constant. Using this formulation, the time step has no effect on the computation since the initial volume is used.
This is the default setting. The element time step is computed at the same time as the internal forces. The characteristic length and the sound speed are computed for each element in every cycle. The computed time step is compared to a minimum time step value and a scale factor is applied to insure a conservative bound. This is the default for brick and quadrilateral elements. This is the default for shell elements.
This only works for shell and brick elements. The nodal stiffness is one half of eigen value from element stiffness matrix; for a truss element this value is equal to the diagonal term of the stiffness matrix. It is computed from the accumulation of element and interface stiffnesses.
These stiffnesses are obtained during internal force computation. For a regular mesh, the element time step and nodal time step conditions are identical. For example, take a truss element, Figure 4. As for the element time step, minimum time step and scale factors are required. The default value is for the scale factor is 0.
If the minimum time step is reached, the analysis can either be stopped or a mass scaling formulation can be applied. In this latter case, mass is added to the affected nodes so that the time step remains constant at the minimum value. This option can be enabled using the same third keyword as used in the element time step control.
It must be checked that added masses do not affect the accuracy of results. If one uses the nodal time step, the element time step is ignored. The interface time step control depends on the type of interface being used. For the interfaces in which the contact conditions are defined by applying kinematic conditions no time step restriction is required. Therefore, the interfaces are stable if a time step scale factor less than or equal to 0.
BeamDyn Driver Input File 4. BeamDyn Primary Input File 4. Blade Input File 4. Echo File 4. Summary File 4. Results File 4. Coordinate Systems 4. Geometrically Exact Beam Theory 4. Linearization Process 4. Damping Forces and Linearization 4. Local Articulation Axes Axis Transformation and Euler Rotations Direction and Phase Angle Conventions Ocean Environmental Conditions.
Ocean Waves Regular Wave Linear Regular Wave Second Order Stokes Wave Irregular Waves Formulated Wave Spectra Pierson-Moskowitz Spectrum Gaussian Spectrum User Defined Wave-Spectrum
0コメント