Victor, I've been thinking (slowly, as that is my way!) about how to incorporate simple Weibull statistics to my FEA results for brittle components so the survival/failure probability is more statistical and reflecting the volume-dependence of the component. My Weibull understanding is not deep, but put simply, wherever there is tensile stress in a brittle material there is a finite probability of failure. This probability is a function of the value of the stress and the volume over which it acts (because if the volume is large there is a higher probability of a fatal flaw that initiates failure).
I think to do this I need to incorporate two material properties (constants) into the integral, namely the Weibull parameter m and the characteristic stress s0. What I am after is integral of (s/s0)^m dV, where s is the stress in the element and dV is the volume of the element. My rudimentary algebra tells me it is a fairly simple matter to get s0 outside of the integral, so at a minimum integral of s^m would be OK. My rudimentary algebra does not see a way to get the ^m bit out of the integral (maybe there is a way), so I would need you to incorporate this into the integral function. Is this a possibility?
I will usually not even have a value for s0, so my analyses would be comparative rather than quantitative, ie., design A is 80% less likely to fail than design B.
I see Mecway will produce tabulated stress values for elements, but I don't see a way to produce tabulated element volume values. Am I wrong? If I could do this, I could do the required calculations in a spreadsheet. I guess deformed volumes would be better, but I guess it will not make much difference. I looked at the liml file but I don't think the stress results are accessible.
There may be a complication: I am not sure whether to include compressive stresses into the integral, my instinct is to exclude them as brittle materials are generally failing in tension. Maybe someone out there knows enough to help me with this, but I think I have someone else who can advise, so I will update if I get better info.
Thanks for reading.
Comments
I don't know about what stress to use or whether to exclude compressive stress or anything. But once you know that, you should be able to incorporate it into the formula, possibly using the heaviside() function to exclude negative values of some variable.
I’m also interested in your problem, and I have tried to implement this statistical approach. I think I have got some progress but there is a point I’m not really sure if it is properly working. Heaviside function is giving me values bigger than 1 inside the elements when I’m working with quadratic elements.
If the Mean_Volume_Integral is computed at the nodes that’s not an issue but if it is done over the Gauss points It may affect the result cause I’m not sure if it is properly eliminating elements under compression. Meanwhile I had to go to linear elements.
Could you make some progress?. I would really like to know how are you approaching the problem.
I share you my temporary file. I’m still working.
Usually, this would just be a discretization error that reduces with mesh refinement as the proportion of affected elements reduces. Do you think that's the case here or it somehow remains significant?
I'm not sure linear elements are better at this, they just keep their value everywhere inside the element bounded by the maximum and minimum node values because they're approximating the step function with a linear one.
I have checked with a materials guy at University of Surrey and he says to exclude the compressive stresses. I haven't worked out how to relate the integrated numbers to failure probability (it seems complicated) but I think even without the detailed understanding we can say that a bigger number is more likely to fail than a smaller.
I’m trying to get it indirectly from available variables but I’m not finding the way.
I think there is no way to sort the need of the volume of each individual element for computing its individual contribution to the overall probability of survival / probability of failure.
Another option would be to mesh the model with all the elements with the same size but I don’t think this is always possible.
However, the Python API gives access to individual element volumes as well as the integral tool, all field variables and all numeric material properties so that could be a better way for producing a single summary result or a list of values per element. But it can't yet generate a new field variable to display graphically.
I’m struggling with the computation of the POF by means of a unique Script.
I already know exactly what to do but I do not have skills with the API environment to build the script yet so I’m firstly solving it in Excel, postprocessing the available information, but some things are missing.
I know and I have computed the volumes of each element, the Gauss point coordinates and their weights, but the Stresses at the Gauss points are missing, and I can’t get to extract them from Mecway.
¿How do I extract the Principal Stresses at the Gauss Points in Mecway ?
My goal is to solve it in Excel , validate the procedure with a known solution and find some angel who help me to build the script and share it for other users.
I suppose another way would be to interpolate the node values back to the Gauss points using element shape functions.
Though I don't know anything about this Weibull statistics, I'm surprised you need to use Gauss points and element volumes. Is it a method which fundamentally requires a discretization of the stress field and it's easiest to use the existing mesh for that? Or does it require some operations which aren't available in Mecway so you have to apply them to the mesh externally?
If you're doing element integration yourself to get around the issue of heaviside() producing values outside the range of [0, 1], I don't think that will really solve the problem since the field can still be outside that range elsewhere in an element even if it's only 0 or 1 at every integration point.
As I have understood, first I need to compute the failure probability of each element by means of a volumetric integral on each finite element.
The field to be integrated is a complex combination of the principal stresses. I don’t see the way to introduce it as a custom solution in MECWAY, and if I could, I would have to integrate, element per element, to get each element probability of survival (Sk). Once each finite element probability of survival is computed, then it needs to be multiplied by its volume (Vk) before adding all together to get the POF (probability of Failure).
Stot= SUMPRODUCT (S1:Sk ; V1:Vk) k= Number of elements
POF=1-Stot
Vk is the way to introduce the volumetric dependance into the POF and how each finite element probability is "transferred" to the overall body.
I will update if I get any progress with ccx and “.dat file (*EL PRINT) ”
Thanks. I really appreciate your help.
NOTE: I have done some testing with the API scripts and I think that or I’m really bad (most probably) or “mw.message( str(model_elements_ids[1] ) )” doesn’t work when there is only one element.
-----------------
import math
model_elements_ids= mw.selected_elements()
mw.message("Volumen del elemento" + str(3) + " is: " + str(mw.volume(1)) + " m3")
mw.message( str(model_elements_ids ) )
mw.message( str(model_elements_ids[1] ) )
-----------------
My understanding of Weibull is based on a single very good talk plus a couple of emails (and I still keep getting it wrong). It's a model that describes a probability of survival for a component under stress. Let's imagine a ceramic component under uniform uniaxial stress. At low stress survival probability might be close to 100%, at high stress survival probability might be close to 0%. The Weibull modulus is a number in the model that controls the sharpness of the transition from high survival probability to low survival probability. A high Weibull parameter gives a sharp transition (highly predictable and repeatable failure stress); a low Weibull parameter indicates unpredictable failure due to material flaw inconsistency.
Lets say the probability of survival of the component at a given stress is P. If we then consider a similar component with 4x the volume, the probability of survival is the conditional probability that all four volumes simultaneously survive, or P^4.
With non-uniform stress in a component you need to integrate the volume-weighted stress for each element rather than integrating the probability of survival. From the integrated stress (tensile only) a survival probability for the component can be inferred. I am still working out how exactly this last idea is applied and I'll post a more mathematical treatment when I have made some progress.
It will involve integrating the tensile stresses over the volume, using the Volume integral function in Mecway on a field formula something like:
heaviside(s.xx) * s.xx + ...
depending which stress components you wish to include.
The integral term in the Weibull formula is actually integral (stress^m) dV, with m being the Weibull modulus for the material, which is something I need to think more about. Straightforward enough to include in the formula, I guess.There are a number of assumptions implicit in the Weibull model and the applicability may be limited because multiaxial stresses are many times more likely to see a failure than their sum would suggest, but you guys know way more about this than I do.
I will post any progress that I make, this is a bit further up the shopping list again now.
PROBABILITY OF FAILURE MODELS IN FINITE ELEMENT ANALYSIS OF BRITTLE MATERIALS
CONSTANTINOS GEORGIADIS
Computers & Structures Vol. 18, No. 3, pp. 537-549, 1984. (Computer subroutines outlined)
It describes the background and numerical method. Let me know if you find the way to implement it in MECWAY. I only saw this possible if all the elements had the same size (volume) so the weighting individual volumes become a constant that could be moved out of the integral. I made some progress with Excel, but it rapidly becomes painful.
NOTE: "A Weibull Brittle Material Failure Model for the ABAQUS Computer Program (Joel Bennett)" talks about a UMAT but I didn't pay too much attention as I don't have programing skills. ¿Maybe someone has alrready developed it in the Calculix forum?
The sum is over the elements in V.
The outer integral is over the volume of a single element V_k.
The inner integral is over what volume ∆V? It refers to Appendix B which describes integrating over a sphere. Is ∆V a sphere within the element?? Are those Gauss points it refers to the element's Gauss points or those of some other smaller shape within it?
I'm thinking you might be able to create a field variable for the inner integral - somehow - and then use the API for the outer integral and sum.
The overall Risk of Rupture in a model is a summatory of individual contributions of the Ne elements.
Each of this individual contribution is a common function evaluated at each element (blue) and weighted by the element volume Vk (green) .
The blue term could be evaluated in Mecway through Solution/ New Table/ Element values ((Sn/S0)^m) , but the element volume Vk to weight each term is not a variable inside the custom formula menu.
I see it like a weighted average using the volume as weighting parameter Vk/AV to take into account that elements with more volume have more relevance in the risk of rupture.
Note: I’m not sure if I can post that image here. ¿?¿?
I still don't understand what the domain of the inner integral is. But if you can evaluate that using the table data, then great.
If the only obstacle is obtaining a list of element volumes, here's a script to generate that:
I guess Mr. GEORGIADIS will not mind if I post some extracts of his work. It is >20 years old paper.
There are three integrals.
-Over the whole volume when we add element by element . (see formula 4.8)
-Over each element, when we integrate through all the Gauss Points (see formula 4.9)
-Over each Gauss Point when we integrate for all directions around each Gauss Point (Spherical formula) .
I guess the sequence would be :
Solve the system.
Obtain the 3 Principal Stresses at the Gauss Quadrature Points. (First problem as Mecway provide the Principal Stresses at the nodes)
If one PS is negative then is set to 0 (heaviside).
Now we are ready to proceed.
1-Element Risk of Rupture S is the result of multiplying the element volume (green) by the integral of the Weibull function inside the element (blue).
2- This Weilbull integral is evaluated in two steps. First Red and Second Blue.
Red integral “T” assigns a value to each Gauss Point. It could be evaluated by means of a custom formula as far as the principal stresses were provided at the Gauss points. The spherical integral is substituted by a numerical scheme. Do not mix T’s Gauss Points and Weights with the Element Gauss Points and Weights. They are different. Use linear squeme to evaluate T as a first try.
3- Once this red custom formula is evaluated at each Gauss Point, MECWAY could do itself the blue integral if the user requests a Table / Element value of this custom formula. Here we have a second potential problem similar to the previous one. I’m not sure if Mecway uses Gauss points or Nodal values when evaluating the “Element value”.
4-Export the table and manually multiply each Element value obtained in 3 by its element volume (formula 4.9).
5-Sum each of the terms obtained in point 4 to get the Risk of Rupture. (see formula 4.8).
I stumbled across this topic because I'm trying to understand how to design parts made out of brittle materials. In my research I've found quite a few papers and having found this topic I'll post some links to see if they help anyone here.
For the time being I'll post those that specifically talk about using computers and particularly FEA.
Here are a couple on some NASA software for analysing brittle materials:
CARES/LIFE Ceramics Analysis and Reliability Evaluation of Structures Life Prediction Program
"The primary objective of this report is to describe a public domain computer program that is coupled with the general purpose finite element codes MSC/NASTRAN(1), ANSYS(2), and ABAQUS(3) for predicting the fast-fracture and/or service lifetime failure probability of monolithic ceramic components...
1 MacNeal-Schwendler Corporation, Los Angeles, CA.
2 Swanson Analysis Systems, Inc.,Houston, PA.
3 Hibbitt,Karlsson and Sorenson, Inc.,Providence, RI."
https://ntrs.nasa.gov/api/citations/20030014949/downloads/20030014949.pdf
Surface Flaw Reliability Analysis of Ceramic Components With the SCARE Finite Element Postprocessor Program
"The SCARE (Structural Ceramics Analysis and Reliability Evaluation) computer proqram on statistical fast fracture reliability analysis with quadratic elements for volume distributed imperfections is enhanced to include the use of linear finite elements and thecapability of designing against concurrent surface flaw induced ceramic component failure. presently coupled as a postprocessor to the MSC/NASTRAN general purpose, finite element analysis program. The improved version now includes the Weibull and Batdorf statistical failure theories for both surface and volume flaw based reliability analysis."
https://ntrs.nasa.gov/api/citations/19870007654/downloads/19870007654.pdf
One from Los Alamos
A Weibull Brittle Material Failure Model for the ABAQUS Computer Program
"A statistical failure theory for brittle materials that traces its origins to the Weibull distribution function is developed for use in the general purpose ABAQUS finite element computer program."http://www-eng.lbl.gov/~shuman/NEXT/MATERIALS&COMPONENTS/Quartz/weibull_ABAQUS.pdf
I've not read them fully because they need a better understanding of the workings of FEA than I have at the moment but I am still trying to get a better understanding of how to analyse stresses in brittle materials so if anyone has any tips I would appreciate any advice I can get.