I’m working on the simulation of human tissues, specifically the simulation of the aortic valve. These types of tissues have highly anisotropic behavior that I want to consider.
In Calculix, the *ELASTIC card with the TYPE=ENGINEERING CONSTANTS allows for performing anisotropic cases with 3 Young’s modulus and 3 Poisson’s coefficients (one for each direction: x, y, z).
The application of this model on a plate in the x-y plane (with one Young’s modulus in the x-direction and the other along the y-direction) is straightforward and works without problems.
The same plate, but in a new plane different from x-y, still works by using the ORIENTATION card.
Now, the application of this process on a curved surface seems more difficult to define the ORIENTATION card.
I have attached a picture of the valve, with the two directions of fibers marked in red. On the right, there is an illustration of a single valve element with fibers viewed from the bottom .
I believe I need to define a local axis(x’,y’,z’) for each element. I can calculate the normal z’ by taking the cross product of the element nodes. My question is, how do I choose the correct rotation to align each of my x’ and y’ local axes along the fibers? Or perhaps there is a better approach.
For information, I’m using shell elements, and the mesh shown in the picture is not the one I’m using for calculations but just a software representation.
Yes. That is basically the same approach I took working with fiber reinforced sandwich structures.
I wrote a script called auto-orient that reads a mesh and creates the needed orientations for elements in specified sets. You can find it here.
The repository contains an example that shows how it works.
This script contains some assumptions (see the README) that might or might not work for you. Feel free to fork it and adapt it to your situation.
Two points w.r.t. those materials. First, where are you going to get the required parameters for that model? Second, it is a nonlinear material. This will result in an iterative calculation that might or might not converge. Personally, I would at least run an analysis with a linear material first, as the manual suggests.
I would suggest explicitly moving to he20r / C3D20R elements. Calculix expands shells to 3D elements internally anyway.
Thank you very much for your quick response and your script! I took the time to read your code and also tried things on my side, so I am only replying now. I apologize for the delay.
I cannot directly use your script for 3D elements because my Calculix is coupled to a fluid solver, and the coupling only works with shell elements for now.
For the fiber-reinforced material, I indeed need to develop an optimization process, but I have already done this before, so it should be fine.
As the manual and you suggested, I will run an analysis with a linear material first. This is what I’m currently doing. However, my main problem is not the material properties themselves but the orientation of each local axis of my element.
I don’t understand how you choose the x’ and y’ orientations in your script. For example, I want to fix the x’ axis along the circumferential fiber of the valve. I can calculate z’ by taking the cross product of the element nodes to define it along the normal.
I’m still examining your code, but could you please summarize for me how you define each axis from a mathematical point of view?
First, in the function set_normals, I calculate the normal vector for the element. For a C3D20 element I calculate three edge vectors (nodes[4]-nodes[0]), (nodes[1] - nodes[0]) and (nodes[3]-nodes[0]). I sort these vectors by length and then create the normal vector by calculating the cross product of the two longest edge vectors. Note that these normal vectors are “normalized” to a length of 1.
[updated 2023-07-14]
In the function write_orientation, the other axes, locx and locy are calculated. I do this by first calculating the dot product of the normal vector with the global x axis. If that dot products is close to 1, the normal lies in the global x direction and the local x and y are set as shown in the code.
Otherwise, the local y direction is the cross product of the normal and the x-axis and the local x-axis is the cross product of the local y direction and the normal.
It is done in this way because I am assuming that (a) the global x axis is the long axis of the shape under investigation and (b) this is 0 degree direction of the laminates.
In your case you will probably need to develop functions to generate the local x and y directions based on the position. From the images it looks like the fiber directions are basically an intersection of the shell with a plane. And furthermore these planes look perpendicular.
So for example you could intersect the plane formed by the normal with one of the fiber planes. That would give you the local x or y direction. Taking the cross product of the local x or y with the normal would give you the other normal axis.
