Some Q&As on FE Theory/Modeling and Simulation of Systems and Devices
[Back to SK’s Home Page]
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. |
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. |
||
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. |
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 Ñ. (kÑu) = 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. |
||
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). - 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. |
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. |
||
1) which version of
FORTRAN did you use and in which environment (Unix, GNU, Fortran compiler,
... )? |
* 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 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
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.
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. The settings are: 1) PW (uniform pressure
load) = 10.0 |
If you put
the integration order parameter, NRED (4th parameter at data line 3) to
either
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. |
||
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. |
Yes |
||
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? |
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. |
||
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 ? |
That is perfectly fine... |
||
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. 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: |
The integration parameter should take care of that..... |
||
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 ? |
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. |
||
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.
|
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. |
||
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 ? |
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." |
||
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 ? |
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. |
||
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 ? |
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.... |
||
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:
What do you think ? |
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.... |