FEAP User Forum
FEAP => Input File Issues => Topic started by: bharathn on April 12, 2019, 05:20:55 AM
-
Dear all, this is a trivial question but nevertheless important.
In the Holzapfel model, the orientations specified are in the global coordinate system? Or in the local? I have a cylindrical ring like structure. There are fibres in the wall. I would like to know whether, when I specifiy 1 0 0 as the fibre direction, they point in the radial direction (denoted by the knot number) or whether they point the general 1 0 0 direction. Note that it is not a perfect cylinder but a wavy one.
My apologies if the question is not worded properly.
Regards,
Bharath
-
In ver8.4 the direction is a global Cartesian direction.
For your model you will have to modify the 'pfiber.f' so the fibers in the reference direction for each element are correct. If the problem is almost cylindrical you may be able to use the coordinates passed in 'elcoor.h' (reference ones) to define the reference orientation. Be sure to check that the element you are using computes the xref(:) values before calling constitution.
You will need something like:
theta = atan2(xref(2) - xv(2),xref(1) - xv(1))
! Calculate new a0 for each point
a(1) = cos(theta)*av(3) - sin(theta)*av(1)
a(2) = sin(theta)*av(3) + cos(theta)*av(1)
a(3) = av(2)
! Normalization
a0(:) = a(:)/sqrt(a(1)**2+a(2)**2+a(3)**2)
-
Dear professor,
What if it is not cylindrical and I want to generalize it? Can I use the same piece of code you provided?
-
You will need to use element to project cylinder basis into element ones, or if not close to a cylinder whatever medical image you might be using.
-
Dear Professor,
I have edited the first few lines of code in pfiber.f to the following:
subroutine pfiber(d, a0, f, j, sig, dd, ftype, xref, xcur) ! xref and xcur have been added as new inputs!!!
c-----[--+---------+---------+---------+---------+---------+---------+-]
c Purpose: Fiber models for finite deformation constitution
c N.B. - Use I_4 not \bar{I_4}
c 1.) Holzapfel-Gasser model
c Psi(I_4) = 0.5*h1/h2*[exp(h2*(I_4 - 1)**2) - 1)]
c I_4 = A_I*C_IJ*A_J = (F_iI*A_I)*(F_iJ*A_J) = a_i*a_i
c a_i = F_iI*A_I
c Inputs:
c d(*) - Material parameters
c a0(3) - Structure vector
c f(3,3)- Deformation and displacement gradient
c j - Jacobian determinant
c ftype = 1 uses F; else uses G = displacement gradeint
c Outputs:
c sig(6) - Cauchy stress
c dd(6,6) - Moduli (current configuration)
c-----[--+---------+---------+---------+---------+---------+---------+-]
implicit none
real*8 d(*), a0(3), f(3,3), j, sig(6), dd(6,6), xref(3)
real*8 xcur(3)
integer ftype, i
real*8 a(3), tau(6), i4, ee, h1,h2, dpsi1, dpsi2, aLocal(3),
& a0Local(3), theta
c Compute current structure tensor
! *************** Am I correct in using xcur and xref? ***************
theta = atan2(xref(2) - xcur(2),xref(1) - xcur(1))
! Calculate new a0 for each point
aLocal(1) = cos(theta)*a0(3) - sin(theta)*a0(1)
aLocal(2) = sin(theta)*a0(3) + cos(theta)*a0(1)
aLocal(3) = a0(2)
! Normalization
a0Local(:) = aLocal(:)/sqrt(aLocal(1)**2+aLocal(2)**2+ &
aLocal(3)**2)
if(ftype.eq.1) then ! Use deformation gradient
a(:) = 0.0d0
do i = 1,3
a(:) = a(:) + f(:,i)*a0Local(i)
end do ! i
else ! Use displacement gradient
a(:) = a0Local(:)
do i = 1,3
a(:) = a(:) + f(:,i)*a0Local(i)
end do ! i
endif
There are a few syntax errors of course but I would like to know if this is what I am supposed to end up doing.
Regards,
Bharath
-
Sorry, the 'xv' was the origin for the cylinder (assumed in x-y plane with cylinder along z-axis). You do not use any xcur, the deformation gradient maps the vector to the current orientation. You just need the alocal to be in the plane of the cylinder surface.
-
Dear professor, I made the changes as you have recommended. Is there any problem I could validate this on?
c~ *************** ENSURE CORRECT LOCAL FIBER ORIENTATION ***************
theta = atan2(xref(2),xref(1))
aLoc(1) = cos(theta)*a0(1) - sin(theta)*a0(2)
aLoc(2) = sin(theta)*a0(1) + cos(theta)*a0(2)
aLoc(3) = a0(3)
c~ Normalization
a0Loc(:) = aLoc(:)/sqrt(aLoc(1)**2+aLoc(2)**2+ aLoc(3)**2)
if(ftype.eq.1) then ! Use deformation gradient
a(:) = 0.0d0
do i = 1,3
a(:) = a(:) + f(:,i)*a0Loc(i)
end do ! i
else ! Use displacement gradient
a(:) = a0(:)
do i = 1,3
a(:) = a(:) + f(:,i)*a0Loc(i)
end do ! i
endif
c~ *************** END ***************
-
Try a cylinder with just circumferential ones and then with radial. Finally at an angle
-
Dear professor,
I have attached the images of fibre orientations in the 100, 010, and 110 directions. I have also attached the input file. Note that the first digit represents the radial direction and the second, the circumferential.
The 100 orientation seems to not have much of an impact on the expansion of the cylinder. The 010 orientation exhibits the greatest resistance to expansion. The 110 is in between. I believe this should be what we are looking for?
One observation is that the 010 orientation seems to not have a uniform displacement distribution as we move radially, unlike the 110 and the 100 orientations. Any idea whether this is realistic?
-
I would expect a circumferential one to still give cylindrical response. I would output the initial orientation of the fibers at the various gauss points to see what they are. You should be able to see if they are purely circumferential and of unit length.
-
I have outputted the fiber orientations along with the coordinates at which they are located. To me they seem alright.
An inspection of the output file for the circumferential system shows a rather poor convergence even though I have allowed for quite a number of tang/utang iterations. Could this be the issue?
-
There might be a problem with the plotting of displacement distribution for a cylinder in FEAP (even if you plot it in a cylindrical coordinate system). You should check your displacement with another software such as Paraview, Tecplot.
By the way, please have a look at the paper where the model was proposed,
https://www.biomech.tugraz.at/index.php/puplications#2000
The fibers are not supposed to resist the compression load meaning they do not work while under compression.
-
Thank you for your suggestion! I even checked the output file manually for the displacements and the circumferential orientation gives me non-cylindrical displacements. Let me check it again!
With regards to not working in compression, in the current case, it is only expansion for the circumferential fibres, isn't it? Or am I missing something? For the radial fibres, they are under compression and therefore do not contribute anything.
-
Dear professors,
I have plotted the topmost layers of the cylinder with the circumferential fibres. I used the demo in http://nurbscalculator.in/ . I have attached the NURBS control point data along with the images.
I have compared the output I got from the simulation with the output I would get if it were a 'perfect' circle. Do you think the difference between the two is large enough to see the non-cylindrical plots that we saw in the FEAP window earlier?
Regards,
Bharath
-
The radial plot seemed to have rather large changes in circumference when you added the fibers. However, it has a lot of symmetry for only a few elements. Something seems to be not working well. I will run your mesh to see if I can spot anything.
Is the correct mesh still the one you posted eariler?
-
If you mean the input file, then yes.
Thank you professor!
-
I note you have the fiber axes at 1 1 0 -- this would mean the specified a0(:) would be 0.701 0.701 0.0 -- to have unit length. Is this what you are using or did you modify what the 1 1 0 means. When I changed to 1 0 0 and ran with cylindrical orientation the cylinder was circular.
-
Dear professor, I had tried the 100 (fibres in the radial direction), 010 (fibres in the circumferential direction) and the 110 (at an angle) configurations. When I use 100, I got the cylindrical results that I had attached earlier. The 010 configuration gave me the non-cylindrical displacement. You had suggested that even for 010, it should provide us with a cylindrical displacement. I had probably attached the input file after having tried the 110 configuration.
With regards to the a0 values, I am using the same system (110 ~ a0 = 0.7071 0.7071 0).
Did you get a cylindrical result when you used the 010 fibre direction?
-
I do get a cylindrical shape as shown
-
Okay, then maybe there is something wrong with my formulation. Do you have any idea if I could use the file from FEAp Ver 8.5? Also, would it be possible to see the output file for the radial case? I would like to confirm where the problem lies.
From my outputs for the 100 and the 010 directions, I see the following discrepancies in the convergence of the residual:
######################### FOR 0 1 0 CONFIGURATION ############################
Computing solution at time 1.0000E-01: Total proportional load 1.0000E-01
Residual norm = 2.9199336E+00 1.0000000E+00 t= 0.87 0.01
End Triangular Decomposition t= 0.87 0.01
Residual norm = 4.8878692E-01 1.6739659E-01 t= 0.88 0.01
End Triangular Decomposition t= 0.88 0.01
Residual norm = 1.9583752E-01 6.7069169E-02 t= 0.89 0.01
End Triangular Decomposition t= 0.89 0.01
Residual norm = 7.7890117E-02 2.6675304E-02 t= 0.90 0.01
End Triangular Decomposition t= 0.90 0.01
Residual norm = 3.1134057E-02 1.0662591E-02 t= 0.91 0.01
End Triangular Decomposition t= 0.91 0.01
Residual norm = 1.2430845E-02 4.2572356E-03 t= 0.92 0.02
End Triangular Decomposition t= 0.92 0.02
Residual norm = 4.9646143E-03 1.7002490E-03 t= 0.93 0.02
End Triangular Decomposition t= 0.93 0.02
Residual norm = 1.9825770E-03 6.7898015E-04 t= 0.94 0.02
End Triangular Decomposition t= 0.94 0.02
Residual norm = 7.9172505E-04 2.7114488E-04 t= 0.95 0.02
End Triangular Decomposition t= 0.95 0.02
Residual norm = 3.1616084E-04 1.0827672E-04 t= 0.96 0.02
End Triangular Decomposition t= 0.96 0.02
Residual norm = 1.2625065E-04 4.3237507E-05 t= 0.97 0.02
End Triangular Decomposition t= 0.97 0.02
Residual norm = 5.0413894E-05 1.7265425E-05 t= 0.98 0.02
End Triangular Decomposition t= 0.98 0.02
Residual norm = 2.0130661E-05 6.8942186E-06 t= 0.99 0.02
End Triangular Decomposition t= 0.99 0.02
Residual norm = 8.0381633E-06 2.7528583E-06 t= 1.00 0.02
End Triangular Decomposition t= 1.00 0.02
Residual norm = 3.2095678E-06 1.0991921E-06 t= 1.00 0.02
End Triangular Decomposition t= 1.01 0.02
Residual norm = 1.2815252E-06 4.3888847E-07 t= 1.02 0.02
End Triangular Decomposition t= 1.02 0.02
Residual norm = 5.1168007E-07 1.7523688E-07 t= 1.02 0.03
End Triangular Decomposition t= 1.02 0.03
Residual norm = 2.0429628E-07 6.9966070E-08 t= 1.03 0.03
End Triangular Decomposition t= 1.03 0.03
Residual norm = 8.1566904E-08 2.7934507E-08 t= 1.04 0.03
End Triangular Decomposition t= 1.04 0.03
Residual norm = 3.2565282E-08 1.1152748E-08 t= 1.05 0.03
End Triangular Decomposition t= 1.05 0.03
The 010 configuration stops the iterations as soon as the residual norm hits 10e-08. However, the 100 configuration continues till it hits 10e-13. Shouldn't the residual be lower for the 010 case before the solver decides to move to the next load step?
######################### FOR 1 0 0 CONFIGURATION ############################
Computing solution at time 1.0000E-01: Total proportional load 1.0000E-01
Residual norm = 2.9199336E+00 1.0000000E+00 t= 0.89 0.02
End Triangular Decomposition t= 0.89 0.02
Residual norm = 1.5349910E+01 5.2569382E+00 t= 0.90 0.02
End Triangular Decomposition t= 0.90 0.02
Residual norm = 5.6304700E+00 1.9282870E+00 t= 0.91 0.02
End Triangular Decomposition t= 0.91 0.02
Residual norm = 4.6535173E+00 1.5937066E+00 t= 0.92 0.02
End Triangular Decomposition t= 0.92 0.02
Residual norm = 3.2524666E-02 1.1138838E-02 t= 0.93 0.02
End Triangular Decomposition t= 0.93 0.02
Residual norm = 5.3513725E-07 1.8327035E-07 t= 0.94 0.02
End Triangular Decomposition t= 0.94 0.02
Residual norm = 2.5086218E-14 8.5913660E-15 t= 0.95 0.02
End Triangular Decomposition t= 0.95 0.02
Thank you so much once again!
-
Who has the license for the copy you are using? It may be possible to update what you have.
-
Dear professor, I am verifying that now and I will get back to you as soon as I find out!
-
Dear professor,
I have sent you a message with further details!