v0.14.0
RVE_Scida.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 
5 cubit.cmd('new')
6 #=============================================================
7 #Geometry Input parametrs
8 #=============================================================
9 a1=0.6
10 a2=0.6
11 a3=0.05
12 ag1=0.02
13 ag2=0.02
14 
15 W_warp=a2-ag2/2
16 W_weft=a1-ag1/2
17 
18 L_RVE=4*a2
19 W_RVE=4*a1
20 H_RVE=2.28*a3
21 
22 dy_w=W_warp/100
23 dy_g=ag2/50
24 dx_w=W_weft/100
25 dx_g=ag1/50
26 
27 autofactor=7;
28 
29 #=============================================================
30 #Create vartices for splines to create splines
31 #=============================================================
32 
33 import math
34 da3=0.003;
35 a3=a3+da3;
36 #---------------------
37 # Fill 1 path (14) top
38 #---------------------
39 coordx=-a1
40 coordy=-3*a2-dy_w
41 for 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 #---------------------
54 coordx=a1
55 coordy=-3*a2-dy_w
56 for 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 #---------------------
68 coordy=-a2
69 coordx=-3*a1-dx_w
70 for 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 #---------------------
82 coordy=a2
83 coordx=-3*a1-dx_w
84 for 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 
97 Tvertices=751;
98 vertices=range(1,Tvertices+1); count1=0;
99 for i in range(0, 4):
100  str1='create curve spline vertex '
101  for j in range(0, 751):
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 
109 a3=a3-da3;
110 #-------------------------------------
111 # Fill 1 front surface top profile (1)
112 #-------------------------------------
113 coordy=3*a2
114 coordx=-2*a1+ag1/2-dx_w
115 A_coef=-a3/2*math.sin(math.pi*ag1/4/a1)-a3/2
116 B_coef=math.pi/(ag1-2*a1)
117 C_coef=-math.pi*(ag1-4*a1)/2/(ag1-2*a1)
118 D_coef=-a3/2*math.sin(math.pi*ag1/4/a1)+a3/2
119 for 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 
126 #----------------------------------------
127 # Fill 1 front surface bottom profile (2)
128 #----------------------------------------
129 A_coef=a3/2
130 B_coef=math.pi/(2*a1)
131 C_coef=0
132 D_coef=a3/2
133 coordy=3*a2
134 coordx=-2*a1+ag1/2-dx_w
135 for i in range(0, 201):
136  coordx=coordx+dx_w;
137  coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef
138  coordz=coordz+da3;
139  str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
140  cubit.cmd(str1)
141 
142 Tvertices=201;
143 vertices=range(1,Tvertices+1); count1=3004;
144 for i in range(0, 2):
145  str1='create curve spline vertex '
146  for j in range(0, 201):
147  str1=str1+' '+ str(count1+vertices[j])
148  cubit.cmd(str1); count1=count1+201;
149 
150 #-------------------------------------
151 # Fill 2 front surface top profile (5)
152 #-------------------------------------
153 coordy=3*a2
154 coordx=ag1/2-dx_w
155 A_coef=a3/2
156 B_coef=math.pi/(2*a1)
157 C_coef=0
158 D_coef=-a3/2
159 for i in range(0, 201):
160  coordx=coordx+dx_w;
161  coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef;
162  coordz=coordz;
163  str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
164  cubit.cmd(str1)
165 
166 #----------------------------------------
167 # Fill 2 front surface bottom profile (6)
168 #----------------------------------------
169 coordy=3*a2
170 coordx=ag1/2-dx_w
171 A_coef=-a3/2*math.sin(math.pi*ag1/4/a1)-a3/2
172 B_coef=math.pi/(ag1-2*a1)
173 C_coef=math.pi*(ag1-4*a1)/2/(ag1-2*a1)
174 D_coef=a3/2*math.sin(math.pi*ag1/4/a1)-a3/2
175 for i in range(0, 201):
176  coordx=coordx+dx_w;
177  coordz=A_coef*math.sin(B_coef*coordx+C_coef)+D_coef;
178  coordz=coordz;
179  str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
180  cubit.cmd(str1)
181 
182 Tvertices=201;
183 vertices=range(1,Tvertices+1); count1=3406;
184 for i in range(0, 2):
185  str1='create curve spline vertex '
186  for j in range(0, 201):
187  str1=str1+' '+ str(count1+vertices[j])
188  cubit.cmd(str1); count1=count1+201;
189 
190 #-------------------------------------
191 # Wrap 1 front surface top profile (9)
192 #-------------------------------------
193 
194 coordx=3*a1
195 coordy=-2*a2+ag2/2-dy_w
196 
197 A_coef=a3/2
198 B_coef=-math.pi/(2*a2)
199 C_coef=0
200 D_coef=-a3/2
201 for i in range(0, 201):
202  coordy=coordy+dy_w;
203  coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef
204  coordz=coordz;
205  str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
206  cubit.cmd(str1)
207 #-----------------------------------------
208 # Warp 1 front surface bottom profile (10)
209 #-----------------------------------------
210 coordx=3*a1
211 coordy=-2*a2+ag2/2-dy_w
212 A_coef=-a3/2*math.sin(math.pi*ag2/4/a2)-a3/2
213 B_coef=-math.pi/(ag2-2*a2)
214 C_coef=math.pi*(ag2-4*a2)/2/(ag2-2*a2)
215 D_coef=a3/2*math.sin(math.pi*ag2/4/a2)-a3/2
216 for i in range(0, 201):
217  coordy=coordy+dy_w;
218  coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef;
219  coordz=coordz;
220  str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
221  cubit.cmd(str1)
222 
223 Tvertices=201;
224 vertices=range(1,Tvertices+1); count1=3808;
225 for i in range(0, 2):
226  str1='create curve spline vertex '
227  for j in range(0, 201):
228  str1=str1+' '+ str(count1+vertices[j])
229  cubit.cmd(str1); count1=count1+201;
230 
231 
232 #-------------------------------------
233 # Wrap 2 front surface top profile (13)
234 #-------------------------------------
235 coordx=3*a1
236 coordy=ag2/2-dy_w
237 A_coef=-a3/2*math.sin(math.pi*ag2/4/a2)-a3/2
238 B_coef=-math.pi/(ag2-2*a2)
239 C_coef=-math.pi*(ag2-4*a2)/2/(ag2-2*a2)
240 D_coef=-a3/2*math.sin(math.pi*ag2/4/a2)+a3/2
241 for i in range(0, 201):
242  coordy=coordy+dy_w;
243  coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef
244  coordz=coordz+da3;
245  str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
246  cubit.cmd(str1)
247 
248 #-----------------------------------------
249 # Warp 2 front surface bottom profile (14)
250 #-----------------------------------------
251 coordx=3*a1
252 coordy=ag2/2-dy_w
253 A_coef=a3/2
254 B_coef=-math.pi/(2*a2)
255 C_coef=0
256 D_coef=a3/2
257 for i in range(0, 201):
258  coordy=coordy+dy_w;
259  coordz=A_coef*math.sin(B_coef*coordy+C_coef)+D_coef
260  coordz=coordz+da3;
261  str1='create vertex ' + str(coordx) +' '+str(coordy)+' '+str(coordz)+' color' ;
262  cubit.cmd(str1)
263 
264 Tvertices=201;
265 vertices=range(1,Tvertices+1); count1=4210;
266 for i in range(0, 2):
267  str1='create curve spline vertex '
268  for j in range(0, 201):
269  str1=str1+' '+ str(count1+vertices[j])
270  cubit.cmd(str1); count1=count1+201;
271 
272 
273 # Create yarn cross-section and
274 # sweep along guideline/path to generate yarn volume
275 #
276 
277 # Fill/weft 1-2
278 cubit.cmd('create surface curve 5 6')
279 cubit.cmd('create surface curve 7 8')
280 # Warp 1-2
281 cubit.cmd('create surface curve 9 10')
282 cubit.cmd('create surface curve 11 12')
283 #
284 cubit.cmd('sweep surface 1 along curve 1')
285 cubit.cmd('sweep surface 2 along curve 2')
286 cubit.cmd('sweep surface 3 along curve 3')
287 cubit.cmd('sweep surface 4 along curve 4')
288 
289 #=============================================================
290 #Deletion of remaining free vertices, curves and bodies
291 #=============================================================
292 
293 cubit.cmd('delete vertex all')
294 cubit.cmd('delete curve all')
295 
296 #=============================================================
297 #Create matrix and imprint and merge fibres and matrix to create common surfaces
298 #=============================================================
299 
300 str1='brick x '+str(L_RVE)+' y '+str(W_RVE)+' z '+str(H_RVE); cubit.cmd(str1)
301 
302 cubit.cmd('intersect volume all keep')
303 cubit.cmd('subtract volume 6 7 8 9 from volume 5 keep')
304 cubit.cmd('delete volume 1 2 3 4 5')
305 cubit.cmd('imprint volume 6 7 8 9 10')
306 cubit.cmd('merge volume 6 7 8 9 10')
307 
308 #=============================================================
309 #Meshing negative x, y and z face of the RVE and coping the same to its positive faces
310 #=============================================================
311 # 1: +x 2: +y 3: +z
312 # 4: -x 5: -y 6: -z
313 cubit.cmd('group "g1" add surface 34 38 52')
314 cubit.cmd('group "g2" add surface 26 30 51')
315 cubit.cmd('group "g3" add surface 47')
316 cubit.cmd('group "g4" add surface 33 37 50')
317 cubit.cmd('group "g5" add surface 25 29 49')
318 cubit.cmd('group "g6" add surface 48')
319 
320 # --------
321 # x plane
322 # --------
323 cubit.cmd('surface 52 scheme trimesh')
324 str1='surface 52 size auto factor '+str(autofactor); cubit.cmd(str1)
325 cubit.cmd('mesh surface 52')
326 cubit.cmd('surface 50 scheme copy source surface 52 source vertex 4661 target vertex 4664 source curve 89 target curve 91 nosmoothing')
327 cubit.cmd('mesh surface 50')
328 
329 cubit.cmd('surface 34 scheme trimesh')
330 str1='surface 34 size auto factor '+str(autofactor); cubit.cmd(str1)
331 cubit.cmd('mesh surface 34')
332 cubit.cmd('surface 33 scheme copy source surface 34 source vertex 4640 target vertex 4637 source curve 55 target curve 53 nosmoothing')
333 cubit.cmd('mesh surface 33')
334 
335 cubit.cmd('surface 38 scheme trimesh')
336 str1='surface 38 size auto factor '+str(autofactor); cubit.cmd(str1)
337 cubit.cmd('mesh surface 38')
338 cubit.cmd('surface 37 scheme copy source surface 38 source vertex 4644 target vertex 4641 source curve 61 target curve 59 nosmoothing')
339 cubit.cmd('mesh surface 37')
340 
341 # --------
342 # y plane
343 # --------
344 #
345 cubit.cmd('surface 51 scheme trimesh')
346 str1='surface 51 size auto factor '+str(autofactor); cubit.cmd(str1)
347 cubit.cmd('mesh surface 51')
348 cubit.cmd('surface 49 scheme copy source surface 51 source vertex 4662 target vertex 4661 source curve 90 target curve 92 nosmoothing')
349 cubit.cmd('mesh surface 49')
350 
351 cubit.cmd('surface 26 scheme trimesh')
352 str1='surface 26 size auto factor '+str(autofactor); cubit.cmd(str1)
353 cubit.cmd('mesh surface 26')
354 cubit.cmd('surface 25 scheme copy source surface 26 source vertex 4631 target vertex 4630 source curve 43 target curve 41 nosmoothing')
355 cubit.cmd('mesh surface 25')
356 
357 cubit.cmd('surface 30 scheme trimesh')
358 str1='surface 30 size auto factor '+str(autofactor); cubit.cmd(str1)
359 cubit.cmd('mesh surface 30')
360 cubit.cmd('surface 29 scheme copy source surface 30 source vertex 4635 target vertex 4634 source curve 49 target curve 47 nosmoothing')
361 cubit.cmd('mesh surface 29')
362 
363 # --------
364 # z plane
365 # --------
366 cubit.cmd('surface 47 scheme trimesh')
367 str1='surface 47 size auto factor '+str(autofactor); cubit.cmd(str1)
368 cubit.cmd('mesh surface 47')
369 cubit.cmd('surface 48 scheme copy source surface 47 source vertex 4664 target vertex 4667 source curve 92 target curve 94 nosmoothing')
370 cubit.cmd('mesh surface 48')
371 
372 #=============================================================
373 # Mesh volume
374 #=============================================================
375 
376 cubit.cmd('volume all scheme Tetmesh')
377 cubit.cmd('mesh volume all')
378 
379 #================================================================================
380 # Defining blocks for elastic, transversely-isotropic and potential flow problems
381 #================================================================================
382 
383 vol=['10', '6,7,8,9', '6', '7', '8', '9']
384 mat=['MAT_ELASTIC_1','MAT_ELASTIC_TRANSISO_1','PotentialFlow_1','PotentialFlow_2','PotentialFlow_3','PotentialFlow_4']
385 for i in range(0, 6):
386  cubit.cmd('set duplicate block elements on')
387  str1='block ' + str(i+1) +' volume '+vol[i]; cubit.cmd(str1)
388  str1='block ' + str(i+1) +' name "'+mat[i] + '"'; cubit.cmd(str1)
389 
390 #=============================================================
391 # Material properties for matrix part
392 #=============================================================
393 
394 cubit.cmd('block 1 attribute count 2')
395 Em=3.4e3; Enu=0.35; #gig to mega as we used dimension in mm
396 
397 Elastic=[str(Em), str(Enu)]
398 for i in range(0, 2):
399  str1='block 1 attribute index ' + str(i+1) +' '+Elastic[i]; cubit.cmd(str1)
400 
401 #=============================================================
402 # Material properties for fibres
403 #=============================================================
404 
405 #to use as isotropic
406 #cubit.cmd('block 2 attribute count 5')
407 #Ep=Em; Ez=Em; nup=Enu; nupz=Enu; Gzp=Em/(2*(1+Enu));
408 
409 cubit.cmd('block 2 attribute count 5')
410 Ep=20.865e3; Ez=58.397e3; nup=0.386; nupz=0.08611; Gzp=8.465e3;
411 
412 TransIso=[str(Ep), str(Ez), str(nup), str(nupz), str(Gzp)]
413 for i in range(0, 5):
414  str1='block 2 attribute index ' + str(i+1) +' '+TransIso[i]; cubit.cmd(str1)
415 
416 #=============================================================
417 # Material properties for interface between fibres and matrix
418 #=============================================================
419 
420 alpha_interf=500
421 cubit.cmd('set duplicate block elements on')
422 str1='block 7 surface 23 24 27 28 31 32 35 36'; cubit.cmd(str1)
423 str1='block 7 name "MAT_INTERF_1"'; cubit.cmd(str1)
424 cubit.cmd('block 7 attribute count 4')
425 str1='block 7 attribute index 1 '+str(alpha_interf); cubit.cmd(str1) #now we use 4 parameters for interface
426 str1='block 7 attribute index 2 '+str(0.0); cubit.cmd(str1)
427 str1='block 7 attribute index 3 '+str(0.0); cubit.cmd(str1)
428 str1='block 7 attribute index 4 '+str(0.0); cubit.cmd(str1)
429 
430 #=============================================================
431 # Defining interfaces
432 #=============================================================
433 
434 Interface=['23', '24','27','28', '31', '32','35','36']
435 for i in range(0, 8):
436  str1='sideset ' + str(i+1) +' surface '+Interface[i]; cubit.cmd(str1)
437  str1='sideset ' + str(i+1) +' name "interface'+str(i+1); cubit.cmd(str1)
438 
439 #=============================================================
440 # Defining pressures for potential flow problem
441 #=============================================================
442 
443 Pres=['25','26','29','30','33','34','37','38']; count=0; count1=len(Interface);
444 for i in range(0, 4):
445  str1='create pressure '+str(count+1)+' on surface '+str(Pres[count])+' magnitude 1'; cubit.cmd(str1)
446  str1='create pressure '+str(count+2)+' on surface '+str(Pres[count+1])+' magnitude -1'; cubit.cmd(str1)
447  str1='sideset '+str(count1+1)+' name "PressureIO_' + str(i+1) + '_1"'; cubit.cmd(str1)
448  str1='sideset '+str(count1+2)+' name "PressureIO_' + str(i+1) + '_2"'; cubit.cmd(str1)
449  count=count+2; count1=count1+2;
450 
451 #=============================================================
452 # Definign zero proessrues for potential flow problem (This should be of the same order as PotentialFlow blocks )
453 #=============================================================
454 
455 #zeroPressureNode=[545, 517, 30, 2]
456 #zeroPressureNode=[660, 618, 157, 115]
457 #zeroPressureNode=[561, 519, 44, 2]
458 # da3=0.003; H_RVE=2.28*a3
459 zeroPressureNode=[407, 389, 20, 2]
460 # da3=0.002; H_RVE=2.22*a3
461 #zeroPressureNode=[409, 391, 20, 2]
462 #zeroPressureNode=[331, 319, 14, 2]
463 for i in range(0, 4):
464  str1='nodeset ' + str(i+1) + ' node ' + str(zeroPressureNode[i]); cubit.cmd(str1)
465  str1='nodeset ' + str(i+1)+' name "ZeroPressure_' + str(i+1)+ '"'; cubit.cmd(str1)
466 
467 #=============================================================
468 # Defining surfaces for dispacement, traction and periodic boundary conditions
469 #=============================================================
470 
471 cubit.cmd('sideset 101 surface 33 37 50 25 29 49 48') # all -ve boundary surfaces for periodic boundary conditions
472 cubit.cmd('sideset 102 surface 34 38 52 26 30 51 47') # all +ve boundary surfaces for periodic boundary conditions
473 cubit.cmd('sideset 103 surface 34 38 52 26 30 51 47 33 37 50 25 29 49 48') # all boundary surfaces
Generated by Doxygen 1.8.17 and hosted at University of Glasgow