Greetings!
I have written a fish script the aim of which is to calculate the number and accumulated length of tension cracks. The generation of blocks, voronoi blocks, zoning, assigning properties is given in UCS11 and fish script is run after. The error which I am getting is shown below.
I understand that dividing by 0 is not possible. However, the contact length must not be 0, it can be very small like 10^-8 m. I dont know what is the error. UCS11 is attached and fish script is given below. UCS11.txt (2.5 KB)
;
fish define FUN1
n_tc = 0 ;set number of tension crack to 0
L_tc = 0 ;set length of tension cracks to 0
bi_c = block.contact.head ; obtain the index no. of first contact in the list of all contacts
loop while bi_c # 0 ; looping is done to check all the contacts one by one
f_n = block.contact.force.normal(bi_c) ; obtain the normal force for the given contact index
L = block.contact.length(bi_c) ; obtain the contact length
bl.contact.prop(bi_c,'tension') = T ; assign tensile strength of the contact to variable T
sigma_n = f_n / L ; calculate the normal stress at the contact
;CALCULATION OF TENSION CRACKS
if sigma_n >= T then
n_tc = n_tc + 1 ; a tension crack is generated and the count is increased by 1 everytime
L_tc = L_tc + L ; calculate the length of the tension crack
else
n_tc = n_tc ; if the criteria is not satisfied, then number of cracks remain unchanged
L_tc = L_tc ; length of cracks remains unchanged
end_if
bi_c = bl.contact.next(bi_c) ;move to the next contact in the list of all contacts
end_loop ; loop ends
end ; fish function ends
;
fish history @FUN1
fish history @n_tc
fish history @L_tc
fish history @L
history interval 1000
block cycle time 0.03
Is it because I took a bigger tolerance corner-round-length in my initial script and I should make it smaller? If I take a very small tolerance then overlap happens and the run stops.
Secondly, if I keep the formula sigma_n = f_n * L, then the error vanishes. However, the criteria is incorrect and hence it cant be done. So I know that the problem is with L value. Kindly help. Thank you.
I’ll give you two answers: one for your actual question and one for what you are actually doing (or not doing) with your evalution.
There is - numerically - a minimum contact length given by block contact length-minimum which defaults as 1/10 the average contact length. I suppose this is working correctly in the simulation, but maybe the associated block.contact.length() fish-function is buggy in this respect. You could specify your own contact length minimum and then use that one also in your fish-function whenever the contact length comes back as zero. If this really is a bug, it could also be helpful for you and others if you send this in as a bug report (maybe on simple test model).
You are trying to count the tensile fractures (and length), right? Then your evaluation makes little sense, because let’s say you assigned a tensile strength of 1 MPa and afterwards a residual strength of 0 MPa. Then the normal stress of that contact will never again be 1 MPa after it has failed in tension, it will always be zero. Your function however will not count this as a tensile crack, because sigma_n is not >= T.
It would make much more sense to use the state logic in UDEC (block.contact state) which will give you all tensile cracks by just checking their state:
0 elastic
1 shear failure now
2 tensile failure now
4 shear failure in past
8 tensile failure in past
This way you can also easily differentiate between shear and tensile cracks.
Thank you Markus. I appreciate the idea given. If I understand point no. 2 correctly, then what you are saying is this:
If for a given contact between a pair of blocks, sigma_n becomes >= T, then the normal stress will be equaled to zero. For the first time, this will be taken as tension crack and the length will get added to tensile crack length. When the normal stress equals 0, it physically implies that the blocks have separated. However, within a sample, the same pair of blocks may come in touch with each other again due to complex movement of particles (or blocks). Since sigma_n has been made zero, this time the calculation will exclude the count and length of any tension cracks formed between them.
Is this the right logic I have understood? Please let me know
Regarding point no 1, I will check and come back here.
As per your suggestion, I made a simpler model with bigger Voronoi blocks to test my commands. The model code is given here. TestModelCode.txt (1.8 KB)
Also, I made a very simple fish function as next step to be run after the above mentioned code
;
fish define fun1
bi_c= block.contact.head
loop while bi_c#0
cs = block.contact.state(bi_c)
if cs = 2
L = block.contact.length(bi_c)
n_tc = n_tc + 1
L_tc = L_tc + L
end_if
bi_c = block.contact.next(bi_c)
end_loop
end
fish set @n_tc = 0 @L_tc = 0
history @bi_c
history @fun1
history @cs
history @L
history @n_tc
history @L_tc
block cycle time 0.001
For now, I just kept cs = 2. I can keep cs = 2 or 8. However, I am getting this error Fish: Bad conversion from 1444159737 index to real as shown below
Initially, I thought if there is any contact index with id 1444159737. I used block contact list position range command to list all the contact indexes and found that such index number does not exist. So I dont have any clue about this error. Is the command block.contact.head working properly? Please help.
Correct, but (and I think this is the thing you are not yet understanding correctly) this has nothing to do with your Fish-function. If you set a certain strength and residual strength for the properties, UDEC will do this by itself.
No, because again your Fish-function has nothing to do with that. Your function will only (later) evaluate what it sees, but cannot correctly interpret what has actually happened based on the stresses alone. That is why I suggested that you use the state variables instead.
This might be because you are trying to make a history of a block pointer (@bi_c) which makes no sense.
Correct, but (and I think this is the thing you are not yet understanding correctly) this has nothing to do with your Fish-function. If you set a certain strength and residual strength for the properties, UDEC will do this by itself.
Thank you very much. I understand now that if I have chosen Mohr-Coulomb model for representing contacts between blocks of the sample, then UDEC is automatically calculating tension and shear cracks. And I just need to output it using the block.contact.state command.
This might be because you are trying to make a history of a block pointer (@bi_c) which makes no sense.
That’s right. I removed it and now the error is gone. Thanks a ton.
As per your suggestion, I made a fish function to calculate number and accumulated length of tension and shear cracks. Here is the code
;
fish define fun1
bi_c= block.contact.head
n_tc=0
L_tc=0
n_sc=0
L_sc=0
loop while bi_c#0
cs = block.contact.state(bi_c)
if cs = 2 or 8
La = block.contact.length(bi_c)
n_tc = n_tc + 1
L_tc = L_tc + La
else if cs = 1 or 4
Lb = block.contact.length(bi_c)
n_sc = n_sc + 1
L_sc = L_sc + Lb
end_if
bi_c = block.contact.next(bi_c)
end_loop
end
history @fun1
history @La
history @Lb
history @n_tc
history @n_sc
history @L_tc
history @L_sc
block mechanical history time-total
block cycle time 0.015
As far as I understand, I am getting the number and length of tension cracks correctly as shown below.
The plots look like your function is only evaluating the “now”-Failures.
Please note that mixed states will have “combined” numbers, e.g. a contact that has shear failed in the past (4) and is still failing now (1) will have a state of 4+1=5. So the shear-relevant states would be 1,4,5 and for tension 2,8,10.
Similarly a contact may have failed both in tension and shear, so you might want to think about how to handle that in your evaluation.
Thanks for providing more information on tension and shear relevant states. Accordingly, I changed the states to 1,4,5 for shear and 2,8,10 for tension. However, the number and length of shear cracks remains unchanged i.e., the same error continues. Any idea?
Are you sure that the function considers all the values attributed to the condition when you are using ‘or’ between your numerical values?
Maybe it is considering the first value only (1 and 2), and that’s why your function is only evaluating the “now” Failures?
One question here is regarding the number and length of shear cracks which decreases at few instances. Why would it decrease?
Secondly, is there a way to make the script shorter by reducing the ‘if’ statement?
Thank you.
Regarding the algorithm, I have also used the same method to take into consideration the different values… I don’t know if there is another shorter method, but this one should work.
Concerning the decrease in the number and length of shear cracks, I also have the same problem.
Maybe because the number of contacts is changing after a fracture. We should try to get the normalized number of cracks (with respect to the number of contacts).
Regarding the contact state, you need to know a little binary math.
def check_state
cs_shear = 1+4
cs_tensile = 2+8
; cs = block.contact.state(bi_c)
cs = 10 ; to demonstrate how it works
if math.and(cs,cs_shear) then
io.out(‘shear’)
endif
if math.and(cs, cs_tensile) then
io.out(‘tensile’)
endif
end @check_state
The number of contacts may change as the blocks move. This can be affected by the command block contact delete-open command. The default setting is to not delete contacts as the blocks separate. The length of a contact can also change as the blocks slide or separate.
It could be a matter of contact deletion, but since deleting is not the default setting it is probably not the problem in your case. A simple way to verify would be to also make a history of the total contact number in your function and watch whether it changes.
I would suspect that the problem lies in those shear-failed contacts that later on also experience tensile failure. These mixed states typically occur later on at higher deformation and shear displacement and would mess with your evaluation logic. Another brute-force method would be to monitor all possible states from 0 to 15 and monitor them separately.
Yes, you are correct. One way which has been followed by Gao and Stead (2014) is to divide the length of failed contacts to the total length of all the contacts and this was introduced as a damage parameter (section 3.1 Page 6 and Fig 7). I am attaching this paper here. I find this paper very useful. http://dx.doi.org/10.1016/j.ijrmms.2014.02.003
If you dont have full access, then let me know. I will send the paper.
Thank you @udecsupport for the code. I will check how it can be applied to my case. @markus I checked all possible states i.e. 1, 4, 5, 2, 8, 10 except 0. I obtained number and length of cracks for each case. I found that the shear cracks corresponding to 4 and 5 decrease with number of cycles. Figure shown below.
As I tried to explain, those are not all possible states. For example, a contact with shear-past and tension-past would have a state of 12 etc.
Also you could use mine and @MariamJ 's hint to observe the total number of contacts during simulation. If the count is constant, then you are simply missing those mixed tension/shear states. If the count is not constant, you could try MariamJs suggestion of a normalized count of fractures.
For those states that have both failed in shear and tension it will however be impossible to know, which failure came first (in order to count it as predominantly shear or tension). I remember modifying the constitutive model, adding some extra variables to store that information.
But after all, it might not be that decisive for the overall evaluation. What you should (probably observe) is predominantly tensile fracturing at the onset of loading and then more and more shear failure later on (at least from the materials I’m familiar with).
Regarding the number of contacts, the count is constant. So the function we are using is just missing some combination of states.
@venkateshmd if you add the values of 6 and 12 to both the tensile and shear failures incrementation, the count seems to be logical, and the curves will no more represent any decrease in the number of failures.
‘6’=2+4 represents the past shear + tensile now failures;
and ‘12’=4+8 represents the past shear + past tensile failures.
For the remaining value of 8+1=9 = tensile past + shear now, I didn’t take it into consideration for now, because I’m not sure if it exists. But we can take it into consideration in all cases, if it doesn’t exist, it won’t affect the results.
