Hi everyone,
To put this in context, , I am making a Matlab algorithm for topology optimization that consider fatigue breaking. I am basing on the Ole Sigmund algorithm, which discretices the 2D piece into a rectangular mesh.
For analysing the fatigue resistance of a piece candidate (discreticed on little squares), I need the stress tensor everywhere, so I decided to compute the stress on the center of each squared element, with FEA analysis, because displacements of the nodes are already computed .
Then, I apply the standard 4-node quadrilateral element stress analys: I know the 2 displacements (u,v) of each 4 node, and I apply interpolation to know the deformation on the center of the element, and so the 3 stress components. I do it with the B matrix and C matrix:
B=0.5/l*[ -1 0 1 0 1 0 -1 0; 0 -1 0 -1 0 1 0 1; -1 -1 -1 1 1 1 1 -1];
(l=side of the element)
C=(E/(1-nu^2))* [ 1 nu 0;
nu 1 0;
0 0 (1-nu)/2 ];
Stress=C*B*nodes_displacements
The problem is that once I compute those stresses, and so the Von Misses stress, the final result doesn't seem to be correct, mostly on parts of the piece being flected.
I also tried taking the 16 closer elements for interpolating the stress on the center of the element, but it didn't work either. I did that last thing with a much bigger B matrix, that has as "input" the 32 displacements of those 16 points around the element:
B=[ -0.0013 0 0.0352 0 -0.0352 0 0.0013 0 -0.0117 0 -0.0117 0 0.0013 0 -0.0352 0 0.0352 0 -0.0013 0 0.0117 0 0.0117 0 -0.3164 0 0.3164 0 0.3164 0 -0.3164 0;
0 -0.0013 0 0.0117 0 0.0117 0 -0.0013 0 0.0352 0 -0.0352 0 0.0013 0 -0.0117 0 -0.0117 0 0.0013 0 -0.0352 0 0.0352 0 -0.3164 0 -0.3164 0 0.3164 0 0.3164;
-0.0013 -0.0013 0.0117 0.0352 0.0117 -0.0352 -0.0013 0.0013 0.0352 -0.0117 -0.0352 -0.0117 0.0013 0.0013 -0.0117 -0.0352 -0.0117 0.0352 0.0013 -0.0013 -0.0352 0.0117 0.0352 0.0117 -0.3164 -0.3164 -0.3164 0.3164 0.3164 0.3164 0.3164 -0.3164 ]/(2*l);
So my question is, what would be the simpliest way to fix this problem, and so to get the correct stress values?
I would prefer a way to fix it using only the 4 nodes around the element, because using 16 makes everything more complicated.
Thanks a lot for your answers.
Your strain-displacement matrix [B] is not correct. The Q4 element has bilinear interpolation functions, so the [B(x,y)] cannot be constant but it is a function of x and y – see any classical finite element textbook (Bathe, Zienkiewicz, Cook). If you want to develop a FE code on your own, I suggest that you give a look at the basic theory on the finite element method ; you will find all the equations you need.
Furthermore, be careful because your model has a sharp intersection that gives rise to a stress singularity (infinite stress) no matter how you refine the mesh.
You can compare your B matrix with Bathe's one (see p.173 of K.-J. Bathe. Finite Element Procedures. PRENTICE HALL, Upper Saddle River, New Jersey 07458, 1996.)
Thanks a lot to both Denis and Victor, I appreciate the effort.
I also read in some books that the B matrix may be (x,y) dependant, as the stress values are different all over the element's surface. However, for my work (fatigue analysis), I need a concrete value of the stress, a number that represents the stress on the element in question. That is why I decided to particularise the matrix for (x,y)=(0,0) , which is the centrer of the element.
So I can't really understand where is the problem of doing this particularisation... Or should I do a different thing to get this stress value?
Let me clarify that I'm getting introduced into this topic, so my questions might seem illogical. Thanks again for your effort.
Hi Miguel. It is my understanding that your B matrix corresponds to the case node1=(-l/2,-l/2), node2=(l/2,-l/2), node3=(l/2,l/2), node4=(-l/2,l/2), nodes_displacements=[u1,v1,u2,v2,u3,v3,u4,v4]. Is it true? I think that you have to verify your algorithm on simplest test problem, say, quadrilateral plate under uniform pressure.
Hello Miguel Diaz ,
Your B matrix must contain the elements dimensions a and b if it is in cartesian coordinates for a (2a*2b) element. If you are using a=b=1, then the coefficient is (1/(4ab))=0.25 and not the cited one (0.5/I). If you are using a (l*l) centered element, then N1(x,y)=(1/(l*l))*(0.5*l-x)(0.5*l-y), N2, N3, N4 : in this case your B matrix seems to be true when verifying all the shape element derivatives in the matrix. You can verify this on the finite element books like my colleagues have mentioned.
The easiest way is to compare your computations with theory is to use one or two elements. In this case you can analytically calculate the stresses components and after that the equivalent stress.
- Badges
- Science method