v0.14.0
Loading...
Searching...
No Matches
RVE_CERL_2.jou
Go to the documentation of this file.
1#!python
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.
4
5cubit.cmd('new')
6#=============================================================
7#Geometry Input parametrs
8#=============================================================
9a1=0.920
10a2=0.920
11a3=0.250
12ag1=0.17
13ag2=0.17
14
15W_warp=a2-ag2/2
16W_weft=a1-ag1/2
17
18L_RVE=4*a2
19W_RVE=4*a1
20H_RVE=2.12*a3
21
22dy_w=W_warp/100
23dy_g=ag2/50
24dx_w=W_weft/100
25dx_g=ag1/50
26
27autofactor=5;
28
29#=============================================================
30#Create vartices for splines to create splines
31#=============================================================
32
33import math
34da3=0.014;
35a3=a3+da3;
36#---------------------
37# Fill 1 path (14) top
38#---------------------
39coordx=-a1
40coordy=-3*a2-dy_w
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):
43 coordy=coordy+dy_w;
44 else:
45 coordy=coordy+dy_g;
46 coordz=a3/2*(math.sin(-math.pi*(coordy)/(2*a2)))+a3/2
47 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
48 cubit.cmd(str1)
49
50
51#---------------------
52# Fill 2 path (12) top
53#---------------------
54coordx=a1
55coordy=-3*a2-dy_w
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):
58 coordy=coordy+dy_w;
59 else:
60 coordy=coordy+dy_g;
61 coordz=a3/2*(math.sin(math.pi*(coordy)/(2*a2)))+a3/2
62 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
63 cubit.cmd(str1)
64
65#---------------------
66# Warp 1 path (2) top
67#---------------------
68coordy=-a2
69coordx=-3*a1-dx_w
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):
72 coordx=coordx+dx_w;
73 else:
74 coordx=coordx+dx_g;
75 coordz=a3/2*(math.sin(math.pi*(coordx)/(2*a1)))+a3/2
76 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
77 cubit.cmd(str1)
78
79#---------------------
80# Warp 2 path (8) top
81#---------------------
82coordy=a2
83coordx=-3*a1-dx_w
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):
86 coordx=coordx+dx_w;
87 else:
88 coordx=coordx+dx_g;
89 coordz=a3/2*(math.sin(-math.pi*(coordx)/(2*a1)))+a3/2
90 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
91 cubit.cmd(str1)
92
93#=============================================================
94#Joint vertices to create splines
95#=============================================================
96
97Tvertices=751;
98vertices=range(1,Tvertices+1); count1=0;
99for i in range(0, 4):
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;
104
105#=============================================================
106#Creating yarns/inclusions with sinusoidal cross sections
107#=============================================================
108
109a3=a3-da3;
110#-------------------------------------
111# Fill 1 front surface top profile (1)
112#-------------------------------------
113coordy=3*a2
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):
120 coordx=coordx+dx_w;
121 coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef
122 coordz=coordz+da3;
123 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
124 cubit.cmd(str1)
125
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
128dz_f1=dzmax-dzmin
129
130#----------------------------------------
131# Fill 1 front surface bottom profile (2)
132#----------------------------------------
133A_coef=a3/2
134B_coef=math.pi/(2*a1)
135C_coef=0
136D_coef=a3/2
137coordy=3*a2
138coordx=-2*a1+ag1/2-dx_w
139for i in range(0, 201):
140 coordx=coordx+dx_w;
141 coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef
142 coordz=coordz+da3;
143 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
144 cubit.cmd(str1)
145
146Tvertices=201;
147vertices=range(1,Tvertices+1); count1=3004;
148for i in range(0, 2):
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;
153
154#-------------------------------------
155# Fill 2 front surface top profile (5)
156#-------------------------------------
157coordy=3*a2
158coordx=ag1/2-dx_w
159A_coef=a3/2
160B_coef=math.pi/(2*a1)
161C_coef=0
162D_coef=-a3/2
163for i in range(0, 201):
164 coordx=coordx+dx_w;
165 coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef;
166 coordz=coordz;
167 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
168 cubit.cmd(str1)
169
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
172dz_f2=dzmax-dzmin
173
174#----------------------------------------
175# Fill 2 front surface bottom profile (6)
176#----------------------------------------
177coordy=3*a2
178coordx=ag1/2-dx_w
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):
184 coordx=coordx+dx_w;
185 coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef;
186 coordz=coordz;
187 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
188 cubit.cmd(str1)
189
190Tvertices=201;
191vertices=range(1,Tvertices+1); count1=3406;
192for i in range(0, 2):
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;
197
198#-------------------------------------
199# Wrap 1 front surface top profile (9)
200#-------------------------------------
201
202coordx=3*a1
203coordy=-2*a2+ag2/2-dy_w
204
205A_coef=a3/2
206B_coef=-math.pi/(2*a2)
207C_coef=0
208D_coef=-a3/2
209for i in range(0, 201):
210 coordy=coordy+dy_w;
211 coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef
212 coordz=coordz;
213 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
214 cubit.cmd(str1)
215
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
218dz_w1=dzmax-dzmin
219
220#-----------------------------------------
221# Warp 1 front surface bottom profile (10)
222#-----------------------------------------
223coordx=3*a1
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):
230 coordy=coordy+dy_w;
231 coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef;
232 coordz=coordz;
233 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
234 cubit.cmd(str1)
235
236Tvertices=201;
237vertices=range(1,Tvertices+1); count1=3808;
238for i in range(0, 2):
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;
243
244
245#-------------------------------------
246# Wrap 2 front surface top profile (13)
247#-------------------------------------
248coordx=3*a1
249coordy=ag2/2-dy_w
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):
255 coordy=coordy+dy_w;
256 coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef
257 coordz=coordz+da3;
258 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
259 cubit.cmd(str1)
260
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
263dz_w2=dzmax-dzmin
264#-----------------------------------------
265# Warp 2 front surface bottom profile (14)
266#-----------------------------------------
267coordx=3*a1
268coordy=ag2/2-dy_w
269A_coef=a3/2
270B_coef=-math.pi/(2*a2)
271C_coef=0
272D_coef=a3/2
273for i in range(0, 201):
274 coordy=coordy+dy_w;
275 coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef
276 coordz=coordz+da3;
277 str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
278 cubit.cmd(str1)
279
280Tvertices=201;
281vertices=range(1,Tvertices+1); count1=4210;
282for i in range(0, 2):
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;
287
288
289# Create yarn cross-section and
290# sweep along guideline/path to generate yarn volume
291#
292
293# Fill/weft 1-2
294cubit.cmd('create surface curve 5 6')
295cubit.cmd('create surface curve 7 8')
296# Warp 1-2
297cubit.cmd('create surface curve 9 10')
298cubit.cmd('create surface curve 11 12')
299
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);
304
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)
309
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)
314
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)
319
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')
328
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')
333
334#=============================================================
335#Deletion of remaining free vertices, curves and bodies
336#=============================================================
337
338cubit.cmd('delete vertex all')
339cubit.cmd('delete curve all')
340
341#=============================================================
342#Create matrix and imprint and merge fibres and matrix to create common surfaces
343#=============================================================
344
345str1='brick x '+str(L_RVE)+' y '+str(W_RVE)+' z '+str(H_RVE); cubit.cmd(str1)
346
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')
352
353#=============================================================
354#Meshing negative x, y and z face of the RVE and coping the same to its positive faces
355#=============================================================
356# 1: +x 2: +y 3: +z
357# 4: -x 5: -y 6: -z
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')
364
365# --------
366# x plane
367# --------
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')
373
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')
379
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')
385
386# --------
387# y plane
388# --------
389#
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')
395
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')
401
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')
407
408# --------
409# z plane
410# --------
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')
416
417#=============================================================
418# Mesh volume
419#=============================================================
420
421cubit.cmd('volume all scheme Tetmesh')
422cubit.cmd('mesh volume all')
423
424#================================================================================
425# Defining blocks for elastic, transversely-isotropic and potential flow problems
426#================================================================================
427
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']
430for i in range(0, 6):
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)
434
435#=============================================================
436# Material properties for matrix part
437#=============================================================
438
439cubit.cmd('block 1 attribute count 2')
440Em=3.4e3; Enu=0.35; #gig to mega as we used dimension in mm
441
442Elastic=[str(Em), str(Enu)]
443for i in range(0, 2):
444 str1='block 1 attribute index ' + str(i+1) +' '+Elastic[i]; cubit.cmd(str1)
445
446#=============================================================
447# Material properties for fibres
448#=============================================================
449
450#to use as isotropic
451#cubit.cmd('block 2 attribute count 5')
452#Ep=Em; Ez=Em; nup=Enu; nupz=Enu; Gzp=Em/(2*(1+Enu));
453
454cubit.cmd('block 2 attribute count 5')
455Ep=19.489e3; Ez=160.755e3; nup=0.415; nupz=0.03395; Gzp=7.393e3;
456
457TransIso=[str(Ep), str(Ez), str(nup), str(nupz), str(Gzp)]
458for i in range(0, 5):
459 str1='block 2 attribute index ' + str(i+1) +' '+TransIso[i]; cubit.cmd(str1)
460
461#=============================================================
462# Material properties for interface between fibres and matrix
463#=============================================================
464
465alpha_interf=500
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)
474
475#=============================================================
476# Defining interfaces
477#=============================================================
478
479Interface=['23', '24','27','28', '31', '32','35','36']
480for i in range(0, 8):
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)
483
484#=============================================================
485# Defining pressures for potential flow problem
486#=============================================================
487
488Pres=['25','26','29','30','33','34','37','38']; count=0; count1=len(Interface);
489for i in range(0, 4):
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;
495
496#=============================================================
497# Definign zero proessrues for potential flow problem (This should be of the same order as PotentialFlow blocks )
498#=============================================================
499
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
505for i in range(0, 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)
508
509#=============================================================
510# Defining surfaces for dispacement, traction and periodic boundary conditions
511#=============================================================
512
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