View Full Version : Anybody familiar with Boyce mode
Hi!
I am coding Boyce model for polymers and had a question very specific to the model. Is there anybody who is familiar with Boyce model and has implemented it sometime?
Waiting for response
Ashu
I'm not sure which Boyce model you mean.
Hi!
This is the model developed by Mary C. Boyce at MIT for modeling PC large strain behavior.
LARGE INELASTIC DEFORMATION OF GLASSY POLYMERS.
PART I: RATE DEPENDENT CONSTITUTIVE MODEL
Mary C. BOYCE, David M. PARKS, Ali S. ARGON
Mechanics of Materials 7 (1988) 15-33 15
North-Holland
Ashu
Ooo- sorry, I'm not really familiar with this model. I think Dave Parks is a member here, maybe he can help.
I am quite familiar with that model, and I have implemented it as a UMAT and VUMAT for ABAQUS.
- Jorgen
Hi Jorgen!
As you know the model is framed in terms of deformation gradient, can I use DFGRD1 as F at any instant. Also, incase you are familiar with the UMAT implementation of Mary Boyce, she doesnt seem to use DFGRD0 or DFGRD1; instead she uses TF and TDTF which she obtains by modifying some of the internal ABAQUS subroutines. Is that really compulsory to do so or one can use DFGRD0 and DFGRD1 instead?
Ashu
Hi Jorgen!
Is it significant to carry out Hughes and Winget rotation while simulating Boyce model in ABAQUS or one can do without it as well?
Actually at first I am trying to get the code running at an approximate level of accuracy and thats exactly why I am a bit reluctant to go for Hughes and Winget rotation.
Ashu
Ashu,
A few comments:
:arrow: I am not sure what Mary Boyce did, but you can use the deformation gradients provided by ABAQUS: DFGRD0, DFGRD1. That should work fine.
:arrow: You don't need the Hughes rotations. There are different ways to integrate the equations. I suggest that you start with a simple forward Euler approach and Explict integration. Once that works, you can switch to something else.
- Jorgen
Hello Jorgen,
The other problem is that stress is a function of log(Felastic). As I obtain Fe=F*inv(Fp), it turns out that sometimes Fe does not have all positive eigenvalues (and thus log operation gives -inf) even if I symmetrize it, and thus the program gives an error.
I start with an identity matrix for Fe, calculate Fp, Dp, and gammap and Fpdot from there. Then I update Fp and recalculate Fe (this all happens at same F) using F=Fe*Fp. This is what sometimes gives me a matrix which has non positive eigenvalues.
Is my approach not correct? What is the way to ensure positivities of eigenvalues?
Ashu
Your approach sounds right. You should not need to make the eigenvalues positive, they should simply come out right automatically.
It sounds like you are using a simple forward Euler integration to get Fp at t+dt. For that to work, you need to make sure that the time increments are small enough.
- Jorgen
Thank you Jorgen,
The back stress in my model is a function of Vp where Fp= Vp * R.
The question is whether there exists any subroutine in ABAQUS which can do the polar decomposition or I have to write such a subroutine by myself?
Also, I am doing implicit integration rather than explicit. Perhaps I could not convey my approach clearly in my last posting. Will write about it in detail if any difficulty arises.
One more thing I wanted to ask is, whether you are familiar with Boyce and Mulliken's model? (the one that deals with high strain rate effects)
Cheers,
Ashu
I'm not sure if there is a built-in model for calculating the polar decomposition. I have developed my own approach that is based on the matrix square root.
And yes, I am familiar with the Boyce-Mulliken model.
- Jorgen
Hello Jorgen!
Mulliken's (or Boyce's as well) model assumes F= FeFp. Now when I start the iterations, I have only F. But inorder to enter the iteration looop, I must have atleast Fe which would start the calculation of Fp, gammap,Dp and so on until and unless the terms satisfy the equations. The question is, what should be taken as the initial value of Fe or Fp to start with.
The second question is, what is the yeild surface or the elastic plastic transition criterion? Should it be tau=epsilon or gammap=epsilon?
The other question is again about the kinematics,
Mulliken's paper says: Fp dot = Inv(Fe)*Dp*Fe*Fp while I think it should be Fp dot = Dp*Fp
I would really appreciate if you can clarify these doubts.
Regards,
Ashu
Here are my answers:
:arrow: At time t=0, Fe=Fp=1. At all other times, you need to keep Fp (or Fe) as a state variable.
:arrow: There is not a strict yield surface. The plastic flow rate is a continuous function of the applied stress.
:arrow: The equation in the Mulliken paper is correct. The key to the equation is that there is a tilde over the Dp. I recommend that you work through the kinematics equations.
- Jorgen
Hello Jorgen!
Thanks for the reply. I had one more confusion. Mulliken's paper says
N=T_prime/norm(T_prime) (Eq 17)
What is meant by norm(T_prime) ?
Does it mean T_prime/(tau*sqrt(2))
where tau= sqrt(0.5*trace(T_prime*T_prime))
Also, if I keep my time increment very small, can I get away implementing only elastic Jacobian, as I think it does not affect the solution but only the time of convergence?
Awaiting your advice
Ashu
Hello Ashu,
T_prim is the deviatoric part of T.
norm(T_prim) is some norm of T_prim. I don't think that their paper specifies what norm they used, but it does not matter much since most norms typically differ only by a constant. You can use your definition of tau as the norm.
You may be able to get a way with an elastic Jacobian. Give it a try and see how it works...
- Jorgen
Hello Jorgen!
I was trying to use SPRIND function to get eigenvalues and eigenvectors of a matrix but apparently it doesnt seem to work. I have 3 types of matrices:
1--> [a b 0;b c 0; 0 0 1]
2--> [a b 0;c d 0;0 0 1]
3--> [a b c;d e f;g h i]
As far as I can understand, the problem is with specifying NDI and NSHR in the arguments but am not able to figure out what to specify for each of the cases. Even SPRIC doesnt seem to give correct principal values.
Can you please tell me how to use this function for these cases ?
Ashu
*SPRIC was a typo for SPRINC..
sorry for that
The variables NDI and NSHR are "Number of DIrect" components and "Number of SHeaR" components. These values are input arguments to the UMAT and VUMAT.
I typically do not use the SPRINC function, for the simple reason that it is not available for other FE programs, and I want my code to work for other FE programs as well. Instead I have written my own functions to extract the eigenvalues and eigenvectors.
I would be very suprised, however, it the SPRINC or SPRIND functions don't work as described in the manuals.
- Jorgen
Hello Jorgen,
I have been trying to get the routine for finding U=sqrt(C) but am facing many problems. I am able to find the eigenvalues of C, and hence the eigenvalues of U (=sqrt (eigenvalues of C)). But I also need the roatation matrix, R, which I could use to rotate the diagonal U matrix to the configuration in which C is prescribed. I could do it for a 2 D case, but have no idea how to obtain it in 3 D.
Can you please tell me how should I go about obtaining the rotation matrix R such that C=R*(diagonal matrix of eigenvalues)*inv(R)
Regards,
Ashu
Hello Jorgen,
I have figured out the answer to my question in the last posting. Sorry about posting that message.
But there is one more thing I would wanna ask. How is one supposed to use PNEWDT ?
If I want to decrease the time increment, I update PNEWDT to a value less than 1. After that, am I supposed to exit from UMAT (tranfer the control to the END) or just keep continuing? Where does the control go once PNEWDT is updated?
Regards,
Ashu
Hello Ashu,
My understanding is that if your UMAT runs into a numerical problem that it cannot fix, then you should set PNEWDT < 1 and directly transfer to the end of the UMAT subroutine. I.e. there is no need to update the stress, state variables, or the Jacobian. What ABAQUS will do is to start over the increment with a smaller time increment size.
- Jorgen
Hello,
I encounter a very strange problem in my UMAT. When I have a look at the DFGRD1 values passed by ABAQUS, I see that it starts with DFGRD1=I. The stress as per the model comes out to be zero but then the message file says that the residual forces and displacements are also zero. Thus the consecutive DFGRD1 values ABAQUS provides are NANQ and the results get screwed. Can you tell me what could probbly be going on that is making the residual forces zero as well??
Ashu
I resolved the problem I referred to in the last posting. Sorry for the posting earlier.
Cheers,
Ashu
Hi Jorgen,
I had a couple of questions regarding my UMAT.
1. When I see the DFGRD1 values supplied by ABAQUS, I find that in the first iteration DFGRD1= I, and some other floating point numbers in second iteration. But after second iteration I get following error message in .msg file:
INCREMENT 1 STARTS. ATTEMPT NUMBER 1, TIME INCREMENT 3.125E-04
***NOTE: USER SUBROUTINE REQUESTS A SMALLER TIME INCREMENT.
***ERROR: USER SUBROUTINE HAS GIVEN A NEGATIVE TIME INCREMENT FACTOR OF 0.00000E+00
I dont understand this error and my program doesnt run coz of the same.
Have you ever encountered such a problem or do you have any idea what could be causing it ?
2. The DFGRD1 values I get in the second iteration are of the order of 1, and so the stress calculated after taking log is of the order of 0.1 MPa. So is it the case that once the above problem gets resolved, I will start gettiing larger values of DFGRD1, so that stress obatined is of the order of 100MPa ?
I would be highly obliged if you can give some comments/suggestion regarding the same.
(1) It sounds like your UMAT subroutine is setting the variable pnewdt=0 before returning from the first increment.
(2) Typically, the variable DFGRD1 will gradually grow with each increment until the final increment size is obtained. Have you checked that the Jacobian that your UMAT calculates is consistent with the material model?
- Jorgen
Hello Jorgen!
As per your advice, I printed out PNEWDT few times and it appears that suddenly after one matrix multiplication function call, it changes its value from 1 to 0. What I dont understand is that I am not doing anything with PNEWDT as such in that function, so how come is it changing? Does ABAQUS itself changes it at any point of time it wants?
Regards,
Ashu
My guess is that you unintentionally overwrite the PNEWDT variable somewhere in your code. For example, if you write to an array past the declared size of the array, then you can create all kinds of strange problems, including overwriting other variables.
- Jorgen
Hello Jorgen,
As you know, the Mulliken Boyce model gives stress for alpha and beta components as:
T=function (V_elastic)
where F_elastic=V_elastic*R_elastic
Dp=gammapdot*N
N=T*'/(tau*sqrt(2))
tau=sqrt(0.5*trace(T'T'))
Total stress sigma=T_alpha+T_beta+T_backstress
My question is that for finding Jacobian, I need:
d(delta_sigma)/d(delta_epsilon)
But in the whole model, epsilon doesnt appear and I am not able to understand, what should I take as epsilon and carry out the partial differentiation? I would appraciate if you can throw some light on how should I go about finding the jacobian for the model.
Regards,
Ashu
Good question. You have a few options when it comes to find the Jacobian matrix for the model:
(1) derive the exact expression for the Jacobian using tensor algebra
(2) derive an approximate expression for the Jacobian using tensor algebra and experience
(3) calculate a numerical approximation of the Jacobian
Note, for many models, including the Mulliken-Boyce model, it is not easy to derive an exact expression for the Jacobian. Most people therefore end up using either option 2 or 3 when determining the Jacobian.
Also note that it is not uncommon for graduate students to simply work with explicit simulations using VUMATs, since these subroutines do not require the Jacobian.
- Jorgen
Hello Jorgen,
Thanks for the response. To put the problem more precisely, inorder to calculate Jacobian, one needs to find out d(delta_sigma)/d(delta_strain).
My question is what is strain in the model. I could see the expression for stress=function(V_elastic), but what is strain??
Should I take V_elastic as strain?? or V_total as strain. I need to figure this out first before I do the differentiation. I would be obliged if you can answer this.
Looking forward to your reply.
Regards,
Ashu
The strain is the total logarithmic (true) strain defined by: E = ln[V]
- Jorgen
Hello Jorgen,
Thanks for the reply. I have successfully written a UMAT (explicit scheme) for only alpha component of the model (beta component still has some problems) but what I observe in the results is that the final stress that comes out of the model is somewhat higher than the applied stress and this error increases with the applied stress. For a stress of 100MPa in the .inp file, the stress calculated by ABAQUS comes out to be 102MPa. The question is whether this discrepency is due to any bug in UMAT or it is due to the way ABAQUS calculates the stresses ?
In the first increment I get a large value of gammapdot_beta (~300) while gammapdot_alf ~ 1e-12. This gives a high value of Fpdot_beta which inturn gives a very small value of Fe_beta (as Fp_beta is large), and thus stresses calculated after taking log come out to be highly negative for beta component. This implies large tau_beta, and inturn gammapdot_beta becomes INF. The alpha component runs fine but beta gets stuck, although I apply the same methodology to both. Could you please give me some suggestion as to if there is some check or scaling to be done??
Awaiting your response..
Ashu
Cheers,
Ashu
Ashu,
What do you mean when you say that the final stress is higher than the applied stress? Can you explain your simulation model, eg. how many elements? what BC and applied loads did you use?
The difference that you describe sounds quite large.
I don't quite follow your description of the different gamma_dots. If you are using an explicit simulation then the time increments are likely very small, and hence I am surprised to see that you got gammapdot_beta about 300. That sounds too large (?) Do you have reasonable material parameters?
- Jorgen
Hello Jorgen!
The way I programmed my UMAT is at any point of time, I get Fp,Fp_dot, s and gammapdot from state variables. At the current instant I evaluate Fe using Fe=F*inv(Fp). Then I calculate stresses, and gammapdot at new time and if the two gammapdot values differ by 20%, I use PNEWDT to reduce the dtime and go back. The equations for the alpha component and beta component are similar.
Now first task was to initialize these quantities at zero time, which I did by assuming Fp=I and gammapdot as their values at zero stress:
gammapdot=gammadot0*exp((-delg/(k*theta))(1-(tau/(s+alpha*p))))
For beta component, from Mulliken's paper I get
gammadot0=3.39E5 k= 1.38E-23 theta=400K s=s0=0.247127,delg=3.769E-20
alpha=0.245 p=tau=0 (as stress is zero at starting)
-->gammapdot = gammadot0*exp(-delg/k*theta) = 367.1867
For alpha component it comes out to be ~1e-12 and runs fine.
The problem is that with this high gammapdot_beta, one gets high Dp=gammapdot_beta*N. This inturn gives very high values of Fp_beta in the next iteration-->low Fe-->highly negative stress_beta-->large tau_beta-->exponentially large gammapdot_beta--> request for a smaller time increment (and this keeps repeating until I get error -->Too may attempts for this increment
The time increment is as low as 0.001.
If I don’t take beta component into account at all, the program runs fine (probably coz gammapdot_alf’s initial value is ~1E-12)
I have a single element (10x10x1 in size) and apply
*Boundary
1, 1,2
4, 1,1
*Cload
2, 1, 100
3, 1, 100
Thus I know what final stress should I get (Total Force/Area = 20 in this case).
I think this should explain my approach and would be obliged to seek your opinion on why the beta component doesn’t run. What is still wrong in my approach ?
Wating for your response...
Regards,
Ashu
The typical elastic modulii taken by running code in Mulliken's thesis are:
E_alf = 1286.7
E_beta = 5.492
These are at strain rates = 0.01 /sec and temperature 400 K
Ashu
Hello Jorgen!
My UMAT doesnt seem to stop at the stress level it should and keeps continuing. I have just one single element and apply a uniaxial tensile force, thus I know the stress I should expect. Instead the UMAT always shows me final stress greater than that value and this keeps on increasing with larger forces. Could you throw some light on as to why such a thing is happening in my code? I have posted the specifications of my .inp and umat file in previous postings, in case u need to look at them.
Regards,
Ashu
I didn't check your values, but I am concerned with the high gammapdot=367. That value suggests that the material is does not have any plastic flow resistane (in network B). Perhaps that is valid at the temperature that you are studying.
Now, if the material parameters are indeed correct (ie no typing error), then you have 2 choices: (1) switch over to an explicit simulation (which is what I think Mulliken did), or (2) implement a more robust integration scheme that works for stiff ODEs.
I guess that your other questions are all related to this issue as well.
- Jorgen
Hello Jorgen!
Thanks for the reply. I would check the E values I am using and get back to you on it later. Meanwhile, I was trying to calculate the Jacobian but am confused as to what should I take as delta(sigma) ?
I would be obliged if you can give a few comments as to how should I approach the calculation.
Regards,
Ashu
At short times and small deformations, and for debugging purposes, you should be able to use an elastic Jacobian.
- Jorgen
Hi!
When one runs a UMAT, how does one come to know the strain rate at which one gets that particular stress strain response is produced. Suppose I run the Mullike's model with one element :
*Node
1, 0., 0.
2, 10., 0.
3, 10., 10.
4, 0., 10.
*Boundary
1, 1,2
4, 1,1
2, 1,1, 16
3, 1,1, 16
*Static
0.001, 1
What is the strain rate corresponding to the response I obtain?
Ashu
Powered by vBulletin® Version 4.2.0 Copyright © 2013 vBulletin Solutions, Inc. All rights reserved.