Some Q&As on FE Theory/Modeling and Simulation of Systems and Devices

[Back to SKs Home Page]

Question:

While using the Finite Element Methodology to model the response of MEMS devices under multi-physics energy domains, the software used may mesh the model into thousands of elements, some of which might be in the nano scale.

However, as explained in your lecture "Simulation of MEMS and Nano Devices" here at AAU, you had indicated that there is a physical limit to the size of problems that can accurately be modeled by FE methodology. In particular you had said at the nano-scale level, strong inter-molecular forces exist which require the use of methodologies based on molecular dynamics and quantum mechanics instead of micro/macro mechanics.

Now, when some of your elements in FE mesh have a size of a few hundred nanometers, wouldn't these intermolecular effects make your FE mesh invalid.

Drs. Enyew at the Addis Ababa University - Dec 2000.

Answer:

This is actually an excellent question.

Now, what you need to remember is the fact that meshing/discretization  is done in an FE model,

1) to represent the geometry accurately and also handle geometrical,  material and loading irregularities in the model

2) to relax the domain so that sufficient degrees-of-freedom (and hence equations) are generated for convergence (more technically speaking, it is to approximate the solution over each element so that the solution is very well represented over the whole regime/domain) .

These are strictly mathematical and modeling requirements. A finer mesh will result in a converged solution. From a global perspective, the fact that we have refined the mesh into smaller elements (some of which can be in nano-scale) better approximates the nodal solution and guarantees convergence. From the element point of view, the displacements (or any of the primary unknowns such as electrostatic potential difference-voltage, stream function, temperature) may represent reality very well. But element-level secondary variables (stresses, electric flux, magnetic flux density, heat flow) may not be realistic due to the molecular forces we mentioned. However, globally, the behavior determined through FE simulations is acceptable.

Do not mistake this for other mathematical programs that perform tasks on loans and financial records.

Question:

I would like to carry out a multi-physics analysis of an electrostatically actuated micro device. I have an in-house FE simulation code for a boundary value problem defined by partial differential equations. What are my options?

Student in MEMS Class - UCSD - Oct 2000.

Answer:

Having FEA code based on general PDEs does of course provide you an advantage. Most physics problems in different regimes are defined by the same differential equation (Ordinary or partial). For example, Poisson's Equation . (ku) = f mathematically defines (models) heat transfer, irrotational flow of an ideal fluid or electrostatics problems. For different energy regimes, k, u and f represent, of course, different physical quantities. (k->conductivity in heat transfer, density in ideal fluid flow, and dielectric constant in electrostatic problems)

For more on physics problems defined by similar mathematical relationships, you may follow the links (courtesy of Reddy's FE Book). [Poisson's Equation] ---- [Differential Equation in 1D].

Therefore by specifying the appropriate constants, you will be able to solve the relevant/corresponding physics problem. More importantly, if your FE code models only one of the energy regims, you can then use analogy to solve for the energy regime you are interested in.

In the case of physics problems where there is coupling between the energy regimes (at best mildly), then you will have to use what is called sequential coupling. In cases like thermal-structural coupling, a one pass coupling will be sufficient where a one-pass of temperature from thermal analysis to a structural analysis will give acceptable results. In other highly coupled problems like fluid-structure interaction, a cyclical/recursive coupling is required since fluid-flow is affected by structural displacements which will inturn will affect the structure displacements. The analysis needs to be done cyclically until equilibrium is achieved between the two regimes. In electrostatically (also in electromagnetically) acuated micro devices, the coupling between the electrostatic/electromagnetic and structural regimes are considerable, therefore you will need to carry out a recursive analysis. Potential changes in the mesh in both regimes have to be addresses by modifying the mesh. If this process could be automated, it will save you a lot of headache.   

Question:

Our research group doesn't have much experience with the fatigue behavior of composite materials. Our aim is to acquire some knowledge about this subject within our department. Till now we have developed a rather simple experimental setup to follow the fatigue deterioration of plain woven E-glass/epoxy composites. Since we hope to describe the fatigue behavior of the material with some kind of law later on, it is important to accurately determine the stresses in the composite material.

Now there are some problems : - an accurate elastic theory to describe the behavior of plain woven composites is not developed yet. Most of the theories are micro mechanical in nature and are only useful for very simple loading configurations. - when talking about deterioration, phenomenon like delaminations and interlaminar stresses are very important in determining the deteriorated state of the material. The classical plate theory does not take into account these effects. - with our experimental setup composite specimens are loaded as a cantilever, so the more common boundary conditions for plates on 4 simply supported edges don't apply to this case.  As a first approach we used the following approximations :

- the plain woven material can be modelled as two cross-ply layers (0/90).
- to acquire the interlaminar stresses, I did some reading in the libraries and found about the "Generalized Laminate Plate Theory" in some articles published in "Computers & Structures" and "The International Journal for Numerical Methods in Engineering". I think this theory would be very appropriate for determining the interlaminar stresses in the plain woven composite. This is important because a damage model should need these stresses to find out whether or not deterioration will continue.

Therefore I would like to ask you :

- are there any commercially available implementations of this theory that run on a plain PC computer with Windows 95/98/NT or do any guidelines exist how you can implement it yourself in C or Fortran ? Besides the composite specimens are simple rectangular plates, so there is no need for a preprocessor or a modeller that can handle various shapes of the specimen. I just need a FEM code that implements the theory and solves the equations. And more important : is it possible to retrieve the source as well, because the degradation law for the fatigue modelling should be implemented later on in this FEM code ?

WVP - University of Ghent, Belgium.

Answer:

A number of years ago, I had written a FORTRAN program for an extension of the GLPT to composite shells.  My doctoral thesis was on modeling Discretely Stiffened Cylindrical Shells using the layerwise Theory and my work was supervised by Dr. Reddy who fist formulated the theory. 

If you need it, I am willing to give you a version of my program (source listing). You can try to understand the theory and usage of the program it (with user manual) and apply it for your purpose. Let me know if you are interested.

Question:

1) which version of FORTRAN did you use and in which environment (Unix, GNU, Fortran compiler, ... )?
2) is there any user manual or a few guidelines relating to the format of the input file ( if you have one, I would be pleased to receive one of them ) ?
3) did you publish some articles that are treating the plate shell theory you are implementing in this code ? Then I could search them in the library over here.

Answer:

* I used FORTRAN 90 and the program is compilable with FORTRAN PowerStation. I had used dynamic array dimensioning which is supported by this compiler. It is easy to remove the DAD and compile it by any compiler.
* I am glad that you received the source code. I hope it will benefit you.
* I had a paper published at Composite structures: Vol 41, 1998, pp 13-26
* I don't have an on-line version of the user's manual. I will mail it to you later. Regards and good luck...

Question:

I have been working for a while with your program "COSHELL". However I have difficulties with the finite element modelling of two problems. I have read your publication "Local behavior of discretely stiffened composite plates and cylindrical shells" in Composite Structures, but I cannot detect the error in my input files. I hope you could give some advice. Don't be discouraged by the number of text lines. I have written
out all my assumptions to clear out the problem for you.  And to further simplify the case for you, I have reduced the problems to a very simple finite element mesh with four four-noded elements. It is a square plate with 200 mm side and 3 mm total thickness (two layers of 1.5 mm and zero angle).

The first problem is the axial tension of an unstiffened orthotropic plate. The attached figure "figures.jpg" shows a schematic drawing of the problem. A uniform tensile stress is applied to one side of the orthotropic plate. The first input file "axial.txt" that I created, uses the following settings:


1) PX (uniform axial load) = -10 Newton/millimeter. I assume that PX is applied to the plane x = 0, like in the figure in your publication; that's why PX is a negative value. However I am not pretty sure what the impact is of setting the hypothesis to "uniform axial load". It is not clear to me in which plane the uniform axial load is defined.
2) NBSF = 0 (array for natural (force) boundary conditions). I have set this value to zero because I thought that the value of PX imposed already the force boundary conditions.
3) the essential boundary conditions I used are the boundary conditions that are shown in your publication for the buckling of the orthotropic stiffened plate.
The calculation runs without any errors, but the strains and displacements are far from uniform through the height of the orthotropic plate. Moreover they are extremely large. I have three questions:

1) if I consequently use "Newton" and "millimeters", is PX then the value of the applied axial stress in N/mm (along the midplane of the orthotropic plate) ? This would be the total force of 2 kN (see figure), divided by 200 mm. I ask this because the displacements are extremely large for such a loading case. The elastic properties should then be expressed in N/mm^2.


2) in your publication you mention that the buckling plate is simply supported and as a consequence you impose w=0 on the boundaries. Don't you constrain the Poisson-effect in the thickness direction then, because w=0 is set for all interface nodes and not only for the midplane?


3) I have set the RADIUS value extremely large for the three interfaces to simulate a flat orthotropic plate, but is there perhaps any parameter to indicate that it is a flat plate ?

Because this model was not good, I tried a second one where I applied the external force myself instead of using PX. I also reduced the force the make the displacements smaller. The settings are:

1) NBSF = 9, because the force will be applied in the 3 interface nodes (= 2 layers (N=2) + 1) for each of the three nodes at the tensile side of the plate.
2) the {q_u} force vector is assigned a value for each of the 3 interface nodes. The problem is that I don't know how I have to distribute the force between each of the nine interface nodes. I think this depends on the degree of the interpolation function.   However, when the calculation is done, the displacements are not uniform again.

The second problem is an unstiffened cylindrical panel that is subjected to internal pressure. It is a part of a composite pressure vessel that is shown on the left in the figure.

The settings are:

1) PW (uniform pressure load) = 10.0
2) ILOD = 1 (internal pressure)
3) the boundary conditions are the ones that you used in your publication for a stiffened cylindrical panel under a point load.
4) NBSF = 0, because I assumed that the PW value already imposes the force boundary conditions.
Again I hesitate what the units of PW are ? Is this N/mm^2, when you use the units "Newton" and "millimeter" ? Because the w-displacements of the node-number 5 are 600 millimeters!

I would be very grateful if you could take a look at these problems, because I did a lot of work trying to solve the difficulties, but it did not succeed.

Answer:

If you put the integration order parameter, NRED (4th parameter at data line 3) to either

0 - full integration
1 - selective reduced for transverse shear terms
2 - selective reduced for transverse normal terms
3 - selected reduced for transverse shear and normal stresses

then you will get reasonable displacements and stresses. I had analyzed it myself successfully. I will attach the result files in the next e-mail.

i.e.,
change
0.0 10.0E0 1.0E00 1 4 1 0 1 1.0 0   to the following line
0.0 10.0E0 1.0E00 1 0 1 0 1 1.0 0


The complete reduced integration could give erroneous results for cases like this depending on the loading, geometry and lamination scheme. However, using full integration could also give errors on a very thin plate like yours (t/b = 100) when subjected to transverse loads.  In general, you need to keep in mind that plate and shell elements are sensitive to t/b ratio, loading, lamination schemes and boundary conditions. "Shear" and "membrane" locking are major problems and you need to keep your eyes open on that. This is common in most plate/shell elements; but the inclusion of the transverse normal terms in layerwise theory introduces one more complicating factor. Hope this helps..... If u still have questions and the numbers still look funny, let me know. I will take care of it.

 

1) PX (uniform axial load) = -10 Newton/millimeter. I assume that PX is applied to the plane x = 0, like in the figure in your publication; that's why PX is a negative value. However I am not pretty sure what the impact is of setting the hypothesis to "uniform axial load". It is not clear to me in which plane the uniform axial load is defined.


Px is applied to plane x=0;

2) NBSF = 0 (array for natural (force) boundary conditions). I have set this value to zero because I thought that the value of PX imposed already the force boundary conditions.

good...

3) the essential boundary conditions I used are the boundary conditions that are shown in your publication for the buckling of the orthotropic stiffened plate. The calculation runs without any errors, but the strains and displacements are far from uniform through the height of the orthotropic plate. Moreover they are extremely large.

Question:

I have three questions:

1) if I consequently use "Newton" and "millimeters", is PX then the value of the applied axial stress in N/mm (along the midplane of the orthotropic plate) ? This would be the total force of 2 kN (see figure), divided by 200 mm. I ask this because the displacements are extremely large for such a loading case. The elastic properties should then be expressed in N/mm^2.

Answer:

Yes

Question:

2) in your publication you mention that the buckling plate is simply supported and as a consequence you impose w=0 on the boundaries. Don't you constrain the Poisson-effect in the thickness direction then, because w=0 is set for all interface nodes and not only for the midplane?

Answer:

All interfaces at that node will be constrained. w=0 at node 1 means that w-1, w-2, and w-3 at node 1 are all 0.0. Hence Poisson effect is also eliminated.

Question:

3) I have set the RADIUS value extremely large for the three interfaces to simulate a flat orthotropic plate, but is there perhaps any parameter to indicate that it is a flat plate ?

Answer:

That is perfectly fine...

Question:

Because this model was not good, I tried a second one where I applied the external force myself instead of using PX. I also reduced the force the make the displacements smaller. The settings are:

1) NBSF = 9, because the force will be applied in the 3 interface nodes (= 2 layers (N=2) + 1) for each of the three nodes at the tensile side of the plate.
2) the {q_u} force vector is assigned a value for each of the 3 interface nodes. The problem is that I don't know how I have to distribute the force between each of the nine interface nodes. I think this depends on the degree of the interpolation function. However, when the calculation is done, the displacements are not uniform again. Thickness-wise, I use linear interpolation....for 2 layers, the distribution will be therefore 25%, 50% and 25%..etc, etc...

The second problem is an un-stiffened cylindrical panel that is subjected to internal pressure. It is a part of a composite pressure vessel that is shown on the left in the figure. The settings are:
1) PW (uniform pressure load) = 10.0
2) ILOD = 1 (internal pressure)
3) the boundary conditions are the ones that you used in your publication for a stiffened cylindrical panel under a point load.
4) NBSF = 0, because I assumed that the PW value already imposes the force boundary conditions.
Again I hesitate what the units of PW are ? Is this N/mm^2, when you use the units "Newton" and "millimeter" ? Because the w-displacements of the node-number 5 are 600 millimeters !

Answer:

The integration parameter should take care of that.....

Question:

Thank you very much for your help. I was not aware of the fact that such an integration scheme could make such trouble. I still use version 1.0 of your program COSHELL, and I saw in your output result files, that there is mentioned something about the Qi3, Qi4 and Qi5 terms concerning the type of integration scheme used.

What does this "reduced integration" refer to ? To the integration of the number of Gauss points ?

Answer:

Yes, "reduced integration" makes huge huge difference in the numerical behavior of models. Most plate and shell elements suffer from this numerical behavior, so called-membrane or shear locking. The reason is that interpolation functions used for these elements end up making the shear or membrane components stiff for thin plates and shells. If interested, I can give you references on this. It is an interesting area; but causes a headache for the unsuspecting user.

Question:

I have one more little question for you. I have been reading some publications of Reddy and you and others about transverse shear and normal stresses.  It is not always clear to me what the implication is of the assumed displacement fields on the transverse continuity and traction-free boundary conditions.
I have tried to summarize it for myself:

1) if w(x,y,z,t) = w_0(x,y,z,t), then the transverse normal strain is zero. The transverse shear strains are constant in each layer. When the transverse shear stresses are derived from the constitutive equations, the transverse shear stress continuity at layer interfaces and the traction-free boundary conditions are NOT satisfied for the transverse shear stresses. Robbins and Reddy reported that layerwise models that neglect transverse normal strain, satisfy the traction-free boundary conditions for transverse shear stresses at the laminate edge only in the integral sense and not in the local sense, despite the level of refinement through the thickness. A more accurate prediction of the transverse shear stresses can be obtained by using the equilibrium method. However Reddy mentions in a study of the ARALL-laminates "Interlaminar stresses (sigma_xz, sigma_yz, sigma_zz) are easily computed from the equilibrium equations of 3D-elasticity when exact analytical solutions are available. An approximate technique is used here to integrate the equilibrium equations and compute the derivatives sigma_xz,z and sigma_yz,z directly from the finite element computation. The equations (with derivatives) are:


sigma_xz,z = - (sigma_xx,x + sigma_xy,y)
sigma_yz,z = - (sigma_xy,x + sigma_yy,y)
(end of citation)


However I don't understand why the condition of the existence of an exact analytical solution has to be fulfilled.  And if you would calculate sigma_zz from the values of sigma_xz and sigma_yz, will the estimate be reasonable ?  Finally Chaudhuri and Seide have demonstrated that the traction-free boundary conditions are not satisfied for the equilibrium method in the case of asymmetrically laminated plates.

Answer:

OK....to start with a good reference is a work done by my friend Don Robbins who was also Dr. Reddy's student. He had published a paper in 92 or 93. You may look for it. If you want, I can make a copy of a report he wrote with Dr. Reddy and mail it to you.

Now, the main thing you need to know when using the equilibrium method to evaluate transverse stresses in FE solutions is that since they are functions of the in-plane stresses themselves, their accuracy depends very much on how accurately the in-plane stresses were calculated to start with. My experience which is supported by the works of others is that at stress-concentration points like supports, the in-plane stresses are not evaluated very accurately. This translates to bad results in your transverse shear/normal stresses at these points. If you, on the other hand, have an exact analytical solution for the in-plane stresses, then your transverse stresses can reasonably be calculated anywhere in your element. Remember stresses from FE analysis are reasonably accurate at Gauss points only. As you move away from the Gauss points, your stresses are interpolated, thereby introducing approximations.

Question:

2) if w(x,y,z,t) = w_0(x,y,z,t) + W(x,y,z,t), like in your model, the transverse normal strain is taken into account as well. When the transverse stresses are derived from the constitutive equations, will transverse stress continuity (both for shear stress and normal stress) and traction free boundary conditions automatically be satisfied ?

Answer:

Depends on the refinement through the thickness. Remember in 3D Layerwise theory (which is essentially descretized 3D), your elements are interpolated both surface and depth-wise. Now if you use enough elements through the thickness, then the conditions will be satisfied very well. If the descretization in the thickness direction is inadequate, then the continuity conditions will not be met. This is very much expected.

Actually, that is why Don Robbins is saying below....

According to Reddy: "The transverse shear stresses, as computed from the constitutive relations, represent the average transverse shear stress within a particular thickness subdivision. However continuous transverse shear stress variation through the thickness can be computed via the 3-D equilibrium equations after having obtained the in-plane stresses accurately." According to Robbins and Reddy: "Note that the form of the natural boundary conditions guarantees that the transverse shear stresses will satisfy the traction-free edge conditions on a local basis as the number of subdivisions through the thickness is increased."

Question:

So, when deriving the transverse shear stresses from the constitutive equations, the traction-free edge condition is not satisfied, because the shear stresses represent the average transverse shear stress within a particular thickness subdivision, but if you increase the number of elements through the thickness, the average transverse shear stress in the boundary layers will tend to zero.  And what when you use the equilibrium method in this case, are then all conditions satisfied ?

Answer:

I will qualify your above statement by saying:

"when deriving the transverse shear stresses from the constitutive equations, the traction-free edge condition CAN BE satisfied" instead of saying "condition is not satisfied".  And I have explained the reason above, i.e., more sub-divisions depth-wise satisfies the reqt.

Question:

In your program COSHELL the in-plane stresses and the transverse normal stress are  written out. Could you use here the constitutive equations or better the equilibrium method to calculate the transverse shear stresses ?

Answer:

Yes, indeed...that is one of the reasons I like the layerwise formulations. It is essentially a 3D solution; hence at Gauss points, you could get accurate in-plane and transverse stresses because these stresses are directly calculated from the constitutive equations.  Both Robbins and mine humble work support that.  I consider this as one of the powers of the layerwise theory....

Question:

The reason why I wrote you about transverse shear stress continuity, is that I started thinking of the way that you apply loads. For example consider a cantilever beam, which is clamped at one side and with three free edges. Suppose that the cantilever beam is loaded with a transverse line loading F [N/m] at its free edge, opposed to the clamped edge. I have attached a drawing "loadcases.jpg" to make it clear. When you want to apply this load in your COSHELL program, there are basically two ways:


(1) you define this line loading as a transverse load, associated with the transverse displacement w_0(x,y,t) of the midplane of the plate.
(2) you define the line loading as a transverse shear stress and you distribute the parabolic shear stress distribution over the various integration points. Then transverse shear stress continuity becomes an important matter.

The same applies to a bending moment M [Nm/m] (see also attached drawing):


(1) you can define the bending moment as two opposite forces, applied to the upper and lower surface of the plate.
(2) you can define the bending moment as a triangular normal stress distribution, which is contributed to the various integration points.

In fact, I am not sure which way to choose in both cases. I have the feeling that the virtual work, done by both definitions, will not be the same, since the interpolation functions of the finite elements play their role. Due to the separate surface- and depth-wise integration, the weight factors of the several integration points through the thickness do not seem to be easily defined in case of the parabolic shear stress distribution or the triangular normal stress distribution.

What do you think ?

Answer:

 

Actually, you have asked a very good question.

In the virtual work and FE formulations, I assume that the load could be applied at any of the interfaces. Which means, as long as the loads are applied with a magnitude that reflects their weight factors, then you could apply them at the right interface.  However, I haven't given the user an option just to specify the load and the program to internally distribute them to the right interfaces (hence degrees-of-freedom) using interpolation functions. This is just like the way you specify uniform load over the surface and the program automatically integrates them and distributes them to the nodes using the shape functions.

I should mention that you seem to have understood this point very well.

Now, given this, what are our options?

Case I - I like your idea of putting the load in the middle; primarily because the parabolic distribution is not done by the program and it is a lot of work for you to calculate the magnitudes at each interface and put them manually. Moreover, to get good results, you may need to use a lot of layers.

One loading pattern I am inclined to consider is applying it at the top; if it is an external load. Think about this....In reality, where is the load applied? At the outer skin or in the middle of the plate?

Case II) Applying the equal and opposite couple at top and bottom fibers will work perfectly fine.

Hope this answers your question....