v0.15.0
Loading...
Searching...
No Matches
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
5cubit.cmd('new')
6#=============================================================
7#Geometry Input parametrs
8#=============================================================
9a1=0.6
10a2=0.6
11a3=0.05
12ag1=0.02
13ag2=0.02
14
15W_warp=a2-ag2/2
16W_weft=a1-ag1/2
17
18L_RVE=4*a2
19W_RVE=4*a1
20H_RVE=2.28*a3
21
22dy_w=W_warp/100
23dy_g=ag2/50
24dx_w=W_weft/100
25dx_g=ag1/50
26
27autofactor=7;
28
29#=============================================================
30#Create vartices for splines to create splines
31#=============================================================
32
33import math
34da3=0.003;
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(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
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
126#----------------------------------------
127# Fill 1 front surface bottom profile (2)
128#----------------------------------------
129A_coef=a3/2
130B_coef=math.pi/(2*a1)
131C_coef=0
132D_coef=a3/2
133coordy=3*a2
134coordx=-2*a1+ag1/2-dx_w
135for 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
142Tvertices=201;
143vertices=range(1,Tvertices+1); count1=3004;
144for 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#-------------------------------------
153coordy=3*a2
154coordx=ag1/2-dx_w
155A_coef=a3/2
156B_coef=math.pi/(2*a1)
157C_coef=0
158D_coef=-a3/2
159for 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#----------------------------------------
169coordy=3*a2
170coordx=ag1/2-dx_w
171A_coef=-a3/2*math.sin(math.pi*ag1/4/a1)-a3/2
172B_coef=math.pi/(ag1-2*a1)
173C_coef=math.pi*(ag1-4*a1)/2/(ag1-2*a1)
174D_coef=a3/2*math.sin(math.pi*ag1/4/a1)-a3/2
175for 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
182Tvertices=201;
183vertices=range(1,Tvertices+1); count1=3406;
184for 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
194coordx=3*a1
195coordy=-2*a2+ag2/2-dy_w
196
197A_coef=a3/2
198B_coef=-math.pi/(2*a2)
199C_coef=0
200D_coef=-a3/2
201for 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#-----------------------------------------
210coordx=3*a1
211coordy=-2*a2+ag2/2-dy_w
212A_coef=-a3/2*math.sin(math.pi*ag2/4/a2)-a3/2
213B_coef=-math.pi/(ag2-2*a2)
214C_coef=math.pi*(ag2-4*a2)/2/(ag2-2*a2)
215D_coef=a3/2*math.sin(math.pi*ag2/4/a2)-a3/2
216for 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
223Tvertices=201;
224vertices=range(1,Tvertices+1); count1=3808;
225for 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#-------------------------------------
235coordx=3*a1
236coordy=ag2/2-dy_w
237A_coef=-a3/2*math.sin(math.pi*ag2/4/a2)-a3/2
238B_coef=-math.pi/(ag2-2*a2)
239C_coef=-math.pi*(ag2-4*a2)/2/(ag2-2*a2)
240D_coef=-a3/2*math.sin(math.pi*ag2/4/a2)+a3/2
241for 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#-----------------------------------------
251coordx=3*a1
252coordy=ag2/2-dy_w
253A_coef=a3/2
254B_coef=-math.pi/(2*a2)
255C_coef=0
256D_coef=a3/2
257for 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
264Tvertices=201;
265vertices=range(1,Tvertices+1); count1=4210;
266for 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
278cubit.cmd('create surface curve 5 6')
279cubit.cmd('create surface curve 7 8')
280# Warp 1-2
281cubit.cmd('create surface curve 9 10')
282cubit.cmd('create surface curve 11 12')
283#
284cubit.cmd('sweep surface 1 along curve 1')
285cubit.cmd('sweep surface 2 along curve 2')
286cubit.cmd('sweep surface 3 along curve 3')
287cubit.cmd('sweep surface 4 along curve 4')
288
289#=============================================================
290#Deletion of remaining free vertices, curves and bodies
291#=============================================================
292
293cubit.cmd('delete vertex all')
294cubit.cmd('delete curve all')
295
296#=============================================================
297#Create matrix and imprint and merge fibres and matrix to create common surfaces
298#=============================================================
299
300str1='brick x '+str(L_RVE)+' y '+str(W_RVE)+' z '+str(H_RVE); cubit.cmd(str1)
301
302cubit.cmd('intersect volume all keep')
303cubit.cmd('subtract volume 6 7 8 9 from volume 5 keep')
304cubit.cmd('delete volume 1 2 3 4 5')
305cubit.cmd('imprint volume 6 7 8 9 10')
306cubit.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
313cubit.cmd('group "g1" add surface 34 38 52')
314cubit.cmd('group "g2" add surface 26 30 51')
315cubit.cmd('group "g3" add surface 47')
316cubit.cmd('group "g4" add surface 33 37 50')
317cubit.cmd('group "g5" add surface 25 29 49')
318cubit.cmd('group "g6" add surface 48')
319
320# --------
321# x plane
322# --------
323cubit.cmd('surface 52 scheme trimesh')
324str1='surface 52 size auto factor '+str(autofactor); cubit.cmd(str1)
325cubit.cmd('mesh surface 52')
326cubit.cmd('surface 50 scheme copy source surface 52 source vertex 4661 target vertex 4664 source curve 89 target curve 91 nosmoothing')
327cubit.cmd('mesh surface 50')
328
329cubit.cmd('surface 34 scheme trimesh')
330str1='surface 34 size auto factor '+str(autofactor); cubit.cmd(str1)
331cubit.cmd('mesh surface 34')
332cubit.cmd('surface 33 scheme copy source surface 34 source vertex 4640 target vertex 4637 source curve 55 target curve 53 nosmoothing')
333cubit.cmd('mesh surface 33')
334
335cubit.cmd('surface 38 scheme trimesh')
336str1='surface 38 size auto factor '+str(autofactor); cubit.cmd(str1)
337cubit.cmd('mesh surface 38')
338cubit.cmd('surface 37 scheme copy source surface 38 source vertex 4644 target vertex 4641 source curve 61 target curve 59 nosmoothing')
339cubit.cmd('mesh surface 37')
340
341# --------
342# y plane
343# --------
344#
345cubit.cmd('surface 51 scheme trimesh')
346str1='surface 51 size auto factor '+str(autofactor); cubit.cmd(str1)
347cubit.cmd('mesh surface 51')
348cubit.cmd('surface 49 scheme copy source surface 51 source vertex 4662 target vertex 4661 source curve 90 target curve 92 nosmoothing')
349cubit.cmd('mesh surface 49')
350
351cubit.cmd('surface 26 scheme trimesh')
352str1='surface 26 size auto factor '+str(autofactor); cubit.cmd(str1)
353cubit.cmd('mesh surface 26')
354cubit.cmd('surface 25 scheme copy source surface 26 source vertex 4631 target vertex 4630 source curve 43 target curve 41 nosmoothing')
355cubit.cmd('mesh surface 25')
356
357cubit.cmd('surface 30 scheme trimesh')
358str1='surface 30 size auto factor '+str(autofactor); cubit.cmd(str1)
359cubit.cmd('mesh surface 30')
360cubit.cmd('surface 29 scheme copy source surface 30 source vertex 4635 target vertex 4634 source curve 49 target curve 47 nosmoothing')
361cubit.cmd('mesh surface 29')
362
363# --------
364# z plane
365# --------
366cubit.cmd('surface 47 scheme trimesh')
367str1='surface 47 size auto factor '+str(autofactor); cubit.cmd(str1)
368cubit.cmd('mesh surface 47')
369cubit.cmd('surface 48 scheme copy source surface 47 source vertex 4664 target vertex 4667 source curve 92 target curve 94 nosmoothing')
370cubit.cmd('mesh surface 48')
371
372#=============================================================
373# Mesh volume
374#=============================================================
375
376cubit.cmd('volume all scheme Tetmesh')
377cubit.cmd('mesh volume all')
378
379#================================================================================
380# Defining blocks for elastic, transversely-isotropic and potential flow problems
381#================================================================================
382
383vol=['10', '6,7,8,9', '6', '7', '8', '9']
384mat=['MAT_ELASTIC_1','MAT_ELASTIC_TRANSISO_1','PotentialFlow_1','PotentialFlow_2','PotentialFlow_3','PotentialFlow_4']
385for 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
394cubit.cmd('block 1 attribute count 2')
395Em=3.4e3; Enu=0.35; #gig to mega as we used dimension in mm
396
397Elastic=[str(Em), str(Enu)]
398for 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
409cubit.cmd('block 2 attribute count 5')
410Ep=20.865e3; Ez=58.397e3; nup=0.386; nupz=0.08611; Gzp=8.465e3;
411
412TransIso=[str(Ep), str(Ez), str(nup), str(nupz), str(Gzp)]
413for 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
420alpha_interf=500
421cubit.cmd('set duplicate block elements on')
422str1='block 7 surface 23 24 27 28 31 32 35 36'; cubit.cmd(str1)
423str1='block 7 name "MAT_INTERF_1"'; cubit.cmd(str1)
424cubit.cmd('block 7 attribute count 4')
425str1='block 7 attribute index 1 '+str(alpha_interf); cubit.cmd(str1) #now we use 4 parameters for interface
426str1='block 7 attribute index 2 '+str(0.0); cubit.cmd(str1)
427str1='block 7 attribute index 3 '+str(0.0); cubit.cmd(str1)
428str1='block 7 attribute index 4 '+str(0.0); cubit.cmd(str1)
429
430#=============================================================
431# Defining interfaces
432#=============================================================
433
434Interface=['23', '24','27','28', '31', '32','35','36']
435for 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
443Pres=['25','26','29','30','33','34','37','38']; count=0; count1=len(Interface);
444for 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
459zeroPressureNode=[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]
463for 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
471cubit.cmd('sideset 101 surface 33 37 50 25 29 49 48') # all -ve boundary surfaces for periodic boundary conditions
472cubit.cmd('sideset 102 surface 34 38 52 26 30 51 47') # all +ve boundary surfaces for periodic boundary conditions
473cubit.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.9.8 and hosted at University of Glasgow