2#This scrip will create RVE including fibres and matrix, which can be used for all
3#three type of boundary conditions including dispalcement, traction and periodic.
6#=============================================================
7#Geometry Input parametrs
8#=============================================================
29#=============================================================
30#Create vartices for splines to create splines
31#=============================================================
41for i in range(0, 751):
42 if (i<101) or (i>150 and i<351) or (i>400 and i<601) or (i>650):
46 coordz=a3/2*(math.sin(-math.pi*(coordy)/(2*a2)))+a3/2
47 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
56for i in range(0, 751):
57 if (i<101) or (i>150 and i<351) or (i>400 and i<601) or (i>650):
61 coordz=a3/2*(math.sin(math.pi*(coordy)/(2*a2)))+a3/2
62 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
70for i in range(0, 751):
71 if (i<101) or (i>150 and i<351) or (i>400 and i<601) or (i>650):
75 coordz=a3/2*(math.sin(math.pi*(coordx)/(2*a1)))+a3/2
76 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
84for i in range(0, 751):
85 if (i<101) or (i>150 and i<351) or (i>400 and i<601) or (i>650):
89 coordz=a3/2*(math.sin(-math.pi*(coordx)/(2*a1)))+a3/2
90 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
93#=============================================================
94#Joint vertices to create splines
95#=============================================================
98vertices=range(1,Tvertices+1); count1=0;
100 str1='create curve spline vertex '
101 for j in range(125, 626):
102 str1=str1+' '+ str(count1+vertices[j])
103 cubit.cmd(str1); count1=count1+751;
105#=============================================================
106#Creating yarns/inclusions with sinusoidal cross sections
107#=============================================================
110#-------------------------------------
111# Fill 1 front surface top profile (1)
112#-------------------------------------
114coordx=-2*a1+ag1/2-dx_w
115A_coef=-a3/2*math.sin(math.pi*ag1/4/a1)-a3/2
116B_coef=math.pi/(ag1-2*a1)
117C_coef=-math.pi*(ag1-4*a1)/2/(ag1-2*a1)
118D_coef=-a3/2*math.sin(math.pi*ag1/4/a1)+a3/2
119for i in range(0, 201):
121 coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef
123 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
126dzmax=A_coef*math.sin(B_coef*(-a1)+C_coef)+D_coef
127dzmin=A_coef*math.sin(B_coef*(-ag1/2)+C_coef)+D_coef
130#----------------------------------------
131# Fill 1 front surface bottom profile (2)
132#----------------------------------------
138coordx=-2*a1+ag1/2-dx_w
139for i in range(0, 201):
141 coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef
143 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
147vertices=range(1,Tvertices+1); count1=3004;
149 str1='create curve spline vertex '
150 for j in range(0, 201):
151 str1=str1+' '+ str(count1+vertices[j])
152 cubit.cmd(str1); count1=count1+201;
154#-------------------------------------
155# Fill 2 front surface top profile (5)
156#-------------------------------------
163for i in range(0, 201):
165 coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef;
167 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
170dzmax=A_coef*math.sin(B_coef*a1+C_coef)+D_coef
171dzmin=A_coef*math.sin(B_coef*ag1/2+C_coef)+D_coef
174#----------------------------------------
175# Fill 2 front surface bottom profile (6)
176#----------------------------------------
179A_coef=-a3/2*math.sin(math.pi*ag1/4/a1)-a3/2
180B_coef=math.pi/(ag1-2*a1)
181C_coef=math.pi*(ag1-4*a1)/2/(ag1-2*a1)
182D_coef=a3/2*math.sin(math.pi*ag1/4/a1)-a3/2
183for i in range(0, 201):
185 coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef;
187 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
191vertices=range(1,Tvertices+1); count1=3406;
193 str1='create curve spline vertex '
194 for j in range(0, 201):
195 str1=str1+' '+ str(count1+vertices[j])
196 cubit.cmd(str1); count1=count1+201;
198#-------------------------------------
199# Wrap 1 front surface top profile (9)
200#-------------------------------------
203coordy=-2*a2+ag2/2-dy_w
206B_coef=-math.pi/(2*a2)
209for i in range(0, 201):
211 coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef
213 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
216dzmax=A_coef*math.sin(B_coef*(-a2)+C_coef)+D_coef
217dzmin=A_coef*math.sin(B_coef*(-ag2/2)+C_coef)+D_coef
220#-----------------------------------------
221# Warp 1 front surface bottom profile (10)
222#-----------------------------------------
224coordy=-2*a2+ag2/2-dy_w
225A_coef=-a3/2*math.sin(math.pi*ag2/4/a2)-a3/2
226B_coef=-math.pi/(ag2-2*a2)
227C_coef=math.pi*(ag2-4*a2)/2/(ag2-2*a2)
228D_coef=a3/2*math.sin(math.pi*ag2/4/a2)-a3/2
229for i in range(0, 201):
231 coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef;
233 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
237vertices=range(1,Tvertices+1); count1=3808;
239 str1='create curve spline vertex '
240 for j in range(0, 201):
241 str1=str1+' '+ str(count1+vertices[j])
242 cubit.cmd(str1); count1=count1+201;
245#-------------------------------------
246# Wrap 2 front surface top profile (13)
247#-------------------------------------
250A_coef=-a3/2*math.sin(math.pi*ag2/4/a2)-a3/2
251B_coef=-math.pi/(ag2-2*a2)
252C_coef=-math.pi*(ag2-4*a2)/2/(ag2-2*a2)
253D_coef=-a3/2*math.sin(math.pi*ag2/4/a2)+a3/2
254for i in range(0, 201):
256 coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef
258 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
261dzmax=A_coef*math.sin(B_coef*a2+C_coef)+D_coef
262dzmin=A_coef*math.sin(B_coef*ag2/2+C_coef)+D_coef
264#-----------------------------------------
265# Warp 2 front surface bottom profile (14)
266#-----------------------------------------
270B_coef=-math.pi/(2*a2)
273for i in range(0, 201):
275 coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef
277 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
281vertices=range(1,Tvertices+1); count1=4210;
283 str1='create curve spline vertex '
284 for j in range(0, 201):
285 str1=str1+' '+ str(count1+vertices[j])
286 cubit.cmd(str1); count1=count1+201;
289# Create yarn cross-section and
290# sweep along guideline/path to generate yarn volume
294cubit.cmd('create surface curve 5 6')
295cubit.cmd('create surface curve 7 8')
297cubit.cmd('create surface curve 9 10')
298cubit.cmd('create surface curve 11 12')
300str1='move surface 1 x 0 y '+str(-a1)+' z '+str(-(a3+da3)/2); cubit.cmd(str1);
301str1='move surface 2 x 0 y '+str(-a1)+' z '+str((a3+da3)/2); cubit.cmd(str1);
302str1='move surface 3 x '+str(-a2)+' y 0 z '+str((a3+da3)/2); cubit.cmd(str1);
303str1='move surface 4 x '+str(-a2)+' y 0 z '+str(-(a3+da3)/2); cubit.cmd(str1);
305str1='curve 1 copy move x '+str(-W_weft)+' y 0 z '+str(-dz_f1);cubit.cmd(str1)
306str1='curve 1 copy move x '+str(W_weft)+' y 0 z '+str(-dz_f1);cubit.cmd(str1)
307str1='curve 2 copy move x '+str(-W_weft)+' y 0 z '+str(-dz_f2);cubit.cmd(str1)
308str1='curve 2 copy move x '+str(W_weft)+' y 0 z '+str(-dz_f2);cubit.cmd(str1)
310str1='curve 3 copy move x 0 y '+str(-W_warp)+' z '+str(-dz_w1);cubit.cmd(str1)
311str1='curve 3 copy move x 0 y '+str(W_warp)+' z '+str(-dz_w1);cubit.cmd(str1)
312str1='curve 4 copy move x 0 y '+str(-W_warp)+' z '+str(-dz_w2);cubit.cmd(str1)
313str1='curve 4 copy move x 0 y '+str(W_warp)+' z '+str(-dz_w2);cubit.cmd(str1)
315str1='surface 1 copy move x 0 y '+str(-W_RVE)+' z 0';cubit.cmd(str1)
316str1='surface 2 copy move x 0 y '+str(-W_RVE)+' z 0';cubit.cmd(str1)
317str1='surface 3 copy move x '+str(-L_RVE)+' y 0 z 0';cubit.cmd(str1)
318str1='surface 4 copy move x '+str(-L_RVE)+' y 0 z 0';cubit.cmd(str1)
320cubit.cmd('create surface 5 22 13 14')
321cubit.cmd('create surface 6 21 13 14')
322cubit.cmd('create surface 7 24 15 16')
323cubit.cmd('create surface 8 23 15 16')
324cubit.cmd('create surface 9 26 17 18')
325cubit.cmd('create surface 10 25 17 18')
326cubit.cmd('create surface 11 28 19 20')
327cubit.cmd('create surface 12 27 19 20')
329cubit.cmd('create volume surface 1 5 9 10')
330cubit.cmd('create volume surface 2 6 11 12')
331cubit.cmd('create volume surface 3 7 13 14')
332cubit.cmd('create volume surface 4 8 15 16')
334#=============================================================
335#Deletion of remaining free vertices, curves and bodies
336#=============================================================
338cubit.cmd('delete vertex all')
339cubit.cmd('delete curve all')
341#=============================================================
342#Create matrix and imprint and merge fibres and matrix to create common surfaces
343#=============================================================
345str1='brick x '+str(L_RVE)+' y '+str(W_RVE)+' z '+str(H_RVE); cubit.cmd(str1)
347cubit.cmd('intersect volume all keep')
348cubit.cmd('subtract volume 18 19 20 21 from volume 17 keep')
349cubit.cmd('delete volume 10 12 14 16 17')
350cubit.cmd('imprint volume 18 19 20 21 22')
351cubit.cmd('merge volume 18 19 20 21 22')
353#=============================================================
354#Meshing negative x, y and z face of the RVE and coping the same to its positive faces
355#=============================================================
358cubit.cmd('group "g1" add surface 34 38 52')
359cubit.cmd('group "g2" add surface 26 30 51')
360cubit.cmd('group "g3" add surface 47')
361cubit.cmd('group "g4" add surface 33 37 50')
362cubit.cmd('group "g5" add surface 25 29 49')
363cubit.cmd('group "g6" add surface 48')
368cubit.cmd('surface 52 scheme trimesh')
369str1='surface 52 size auto factor '+str(autofactor); cubit.cmd(str1)
370cubit.cmd('mesh surface 52')
371cubit.cmd('surface 50 scheme copy source surface 52 source vertex 4697 target vertex 4700 source curve 113 target curve 115 nosmoothing')
372cubit.cmd('mesh surface 50')
374cubit.cmd('surface 34 scheme trimesh')
375str1='surface 34 size auto factor '+str(autofactor); cubit.cmd(str1)
376cubit.cmd('mesh surface 34')
377cubit.cmd('surface 33 scheme copy source surface 34 source vertex 4674 target vertex 4673 source curve 78 target curve 80 nosmoothing')
378cubit.cmd('mesh surface 33')
380cubit.cmd('surface 38 scheme trimesh')
381str1='surface 38 size auto factor '+str(autofactor); cubit.cmd(str1)
382cubit.cmd('mesh surface 38')
383cubit.cmd('surface 37 scheme copy source surface 38 source vertex 4678 target vertex 4677 source curve 84 target curve 86 nosmoothing')
384cubit.cmd('mesh surface 37')
390cubit.cmd('surface 51 scheme trimesh')
391str1='surface 51 size auto factor '+str(autofactor); cubit.cmd(str1)
392cubit.cmd('mesh surface 51')
393cubit.cmd('surface 49 scheme copy source surface 51 source vertex 4698 target vertex 4697 source curve 114 target curve 116 nosmoothing')
394cubit.cmd('mesh surface 49')
396cubit.cmd('surface 26 scheme trimesh')
397str1='surface 26 size auto factor '+str(autofactor); cubit.cmd(str1)
398cubit.cmd('mesh surface 26')
399cubit.cmd('surface 25 scheme copy source surface 26 source vertex 4666 target vertex 4665 source curve 68 target curve 66 nosmoothing')
400cubit.cmd('mesh surface 25')
402cubit.cmd('surface 30 scheme trimesh')
403str1='surface 30 size auto factor '+str(autofactor); cubit.cmd(str1)
404cubit.cmd('mesh surface 30')
405cubit.cmd('surface 29 scheme copy source surface 30 source vertex 4670 target vertex 4669 source curve 74 target curve 72 nosmoothing')
406cubit.cmd('mesh surface 29')
411cubit.cmd('surface 47 scheme trimesh')
412str1='surface 47 size auto factor '+str(autofactor); cubit.cmd(str1)
413cubit.cmd('mesh surface 47')
414cubit.cmd('surface 48 scheme copy source surface 47 source vertex 4700 target vertex 4703 source curve 116 target curve 118 nosmoothing')
415cubit.cmd('mesh surface 48')
417#=============================================================
419#=============================================================
421cubit.cmd('volume all scheme Tetmesh')
422cubit.cmd('mesh volume all')
424#================================================================================
425# Defining blocks for elastic, transversely-isotropic and potential flow problems
426#================================================================================
428vol=['22', '18,19,20,21', '18', '19', '20', '21']
429mat=['MAT_ELASTIC_1','MAT_ELASTIC_TRANSISO_1','PotentialFlow_1','PotentialFlow_2','PotentialFlow_3','PotentialFlow_4']
431 cubit.cmd('set duplicate block elements on')
432 str1='block ' + str(i+1) +' volume '+vol[i]; cubit.cmd(str1)
433 str1='block ' + str(i+1) +' name "'+mat[i] + '"'; cubit.cmd(str1)
435#=============================================================
436# Material properties for matrix part
437#=============================================================
439cubit.cmd('block 1 attribute count 2')
440Em=3.4e3; Enu=0.35; #gig to mega as we used dimension in mm
442Elastic=[str(Em), str(Enu)]
444 str1='block 1 attribute index ' + str(i+1) +' '+Elastic[i]; cubit.cmd(str1)
446#=============================================================
447# Material properties for fibres
448#=============================================================
451#cubit.cmd('block 2 attribute count 5')
452#Ep=Em; Ez=Em; nup=Enu; nupz=Enu; Gzp=Em/(2*(1+Enu));
454cubit.cmd('block 2 attribute count 5')
455Ep=19.489e3; Ez=160.755e3; nup=0.415; nupz=0.03395; Gzp=7.393e3;
457TransIso=[str(Ep), str(Ez), str(nup), str(nupz), str(Gzp)]
459 str1='block 2 attribute index ' + str(i+1) +' '+TransIso[i]; cubit.cmd(str1)
461#=============================================================
462# Material properties for interface between fibres and matrix
463#=============================================================
466cubit.cmd('set duplicate block elements on')
467str1='block 7 surface 23 24 27 28 31 32 35 36'; cubit.cmd(str1)
468str1='block 7 name "MAT_INTERF_1"'; cubit.cmd(str1)
469cubit.cmd('block 7 attribute count 4')
470str1='block 7 attribute index 1 '+str(alpha_interf); cubit.cmd(str1) #now we use 4 parameters for interface
471str1='block 7 attribute index 2 '+str(0.0); cubit.cmd(str1)
472str1='block 7 attribute index 3 '+str(0.0); cubit.cmd(str1)
473str1='block 7 attribute index 4 '+str(0.0); cubit.cmd(str1)
475#=============================================================
477#=============================================================
479Interface=['23', '24','27','28', '31', '32','35','36']
481 str1='sideset ' + str(i+1) +' surface '+Interface[i]; cubit.cmd(str1)
482 str1='sideset ' + str(i+1) +' name "interface'+str(i+1); cubit.cmd(str1)
484#=============================================================
485# Defining pressures for potential flow problem
486#=============================================================
488Pres=['25','26','29','30','33','34','37','38']; count=0; count1=len(Interface);
490 str1='create pressure '+str(count+1)+' on surface '+str(Pres[count])+' magnitude 1'; cubit.cmd(str1)
491 str1='create pressure '+str(count+2)+' on surface '+str(Pres[count+1])+' magnitude -1'; cubit.cmd(str1)
492 str1='sideset '+str(count1+1)+' name "PressureIO_' + str(i+1) + '_1"'; cubit.cmd(str1)
493 str1='sideset '+str(count1+2)+' name "PressureIO_' + str(i+1) + '_2"'; cubit.cmd(str1)
494 count=count+2; count1=count1+2;
496#=============================================================
497# Definign zero proessrues for potential flow problem (This should be of the same order as PotentialFlow blocks )
498#=============================================================
500#zeroPressureNode=[178, 172, 8, 2] # autofactor = 8
501#zeroPressureNode=[233, 223, 12, 2] # autofactor = 7
502#zeroPressureNode=[290, 276, 16, 2] # autofactor = 6
503zeroPressureNode=[366, 344, 24, 2] # autofactor = 5
504# zeroPressureNode=[816, 782, 36, 2] # autofactor = 4
506 str1='nodeset ' + str(i+1) + ' node ' + str(zeroPressureNode[i]); cubit.cmd(str1)
507 str1='nodeset ' + str(i+1)+' name "ZeroPressure_' + str(i+1)+ '"'; cubit.cmd(str1)
509#=============================================================
510# Defining surfaces for dispacement, traction and periodic boundary conditions
511#=============================================================
513cubit.cmd('sideset 101 surface 33 37 50 25 29 49 48') # all -ve boundary surfaces for periodic boundary conditions
514cubit.cmd('sideset 102 surface 34 38 52 26 30 51 47') # all +ve boundary surfaces for periodic boundary conditions
515cubit.cmd('sideset 103 surface 34 38 52 26 30 51 47 33 37 50 25 29 49 48') # all boundary surfaces