
[Sponsors] 
July 21, 2013, 10:15 
Mesh UDF problem

#1 
New Member
Join Date: Jun 2013
Posts: 15
Rep Power: 5 
Hello everyone,
I have strange problem with results of my UDF. As you can see, my UDF loops over a cell fluid cell thread and inside every cell does a loop over all faces. In the face loop it takes values from the adjacent cells (dudux1, ...) and performs simple operations (SRT1, Sface, multiplication by flux/density  still inside the face loop). The result of the face loop is tensor DS. As of now, I have 2 problems: 1) As you can see in thumbnail DSPartition there is a slice which has reversed sign. The contours have right shape, but opposite colours. I wondered whether it may be caused be a partition crossing periodic boundary condition (Partition1) so i repartitioned the mesh (Partition2) and repeat calculations. It turned out to be a dead end, since DS for Partition2 looks almost exactly the same. My next idea was as follows: As I understand it, it isn't precisely defined which cell is C0 and which one is C1 (adjacent cell). Thus the flux between them can be taken with opposite sign. So I just fabsed the flux. Resulting distribution is shown in DSfabsflux. It looks nice and smooth, but what are the chances that it is correct? Does anyone know reliable solution to a problem like this? 2) There is clear problem with my periodic boundary condition. It seems that the value is set to 0 there. Since my UDF is to set the value to zero on all wall faces, I thought that Fluent treats pbc as a wall and tried to include appropriate condition in the if sentence, based on the ID value of the pbc. But it failed horribly. Can somebody push me in the right direction as how to deal with this? Also, in the DS for both Partition1 and Partition2 there are horizontal lines of 0 value at the level of the rib. They seem to vanish in the DSfabsflux. May this all situation be connected somehow with the way how the model is meshed? (Thumbnail: Mesh). For the record, according to fluent:  Minimum Orthogonal Quality = 9.99873e01  Maximum Aspect Ratio = 3.52643e+01 Regards, kornetka Code:
DEFINE_ADJUST(adj_turb_var, domain) { int ii,jj,kk,n, zone_ID; real dens, dens1,flux; real rho,mu,tke,omega,S,R,sum,HI,f_beta,S_dash,a2,omega_tilde,mut; real dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz; real dudx1,dudy1,dudz1,dvdx1,dvdy1,dvdz1,dwdx1,dwdy1,dwdz1; real SRT[3][3], SRT1[3][3], Sface[3][3]; real DS[3][3], DS2[3][3]; /* Parallelization */ #if !RP_Host face_t f=1; cell_t c,c1=1; Thread *t,*t1=NULL,*tf=NULL; t = Lookup_Thread(domain, 2); thread_loop_c (t, domain) begin_c_loop(c,t) { dens = 0.0; dens1 = 0.0; flux = 0.0; DS[0][0] = 0.0; DS[0][1] = 0.0; DS[0][2] = 0.0; DS[1][0] = 0.0; DS[1][1] = 0.0; DS[1][2] = 0.0; DS[2][0] = 0.0; DS[2][1] = 0.0; DS[2][2] = 0.0; C_UDSI(c,t,TKE) = MAX(C_UDSI(c,t,TKE), 0.0); C_UDSI(c,t,TKE) = MIN(C_UDSI(c,t,TKE), 10000.0); C_UDSI(c,t,OM) = MAX(C_UDSI(c,t,OM), 0.0); C_UDSI(c,t,OM) = MIN(C_UDSI(c,t,OM), 1.0e15); rho = C_R(c,t); mu = C_MU_L(c,t); tke = MAX(C_UDSI(c,t,TKE),SNUM); omega = MAX(C_UDSI(c,t,OM),SNUM); S = Strainrate_Mag(c,t); R = Rotationrate_Mag(c,t); /* Strain rate tensor in cell centre */ dudx = C_U_G(c,t)[0]; dudy = C_U_G(c,t)[1]; dudz = C_U_G(c,t)[2]; dvdx = C_V_G(c,t)[0]; dvdy = C_V_G(c,t)[1]; dvdz = C_V_G(c,t)[2]; dwdx = C_W_G(c,t)[0]; dwdy = C_W_G(c,t)[1]; dwdz = C_W_G(c,t)[2]; SRT[0][0] = 0.5*(dudx + dudx); SRT[0][1] = 0.5*(dudy + dvdx); SRT[0][2] = 0.5*(dudz + dwdx); SRT[1][0] = 0.5*(dvdx + dudy); SRT[1][1] = 0.5*(dvdy + dvdy); SRT[1][2] = 0.5*(dvdz + dwdy); SRT[2][0] = 0.5*(dwdx + dudz); SRT[2][1] = 0.5*(dwdy + dvdz); SRT[2][2] = 0.5*(dwdz + dwdz); c_face_loop(c,t,n) { f = C_FACE(c,t,n); tf = C_FACE_THREAD(c,t,n); zone_ID = THREAD_ID(tf); /* Strain rate tensor in adjacent cell, if it exists */ if (!BOUNDARY_FACE_THREAD_P(tf)) && zone_ID != 4) { /* Get cell on other side of face */ c1 = F_C1(f,tf); t1 = F_C1_THREAD(f,tf); dens = C_R(c,t); dens1 = C_R(c1,t1); dens = (dens + dens1)/2.0; dudx1 = C_U_G(c1,t)[0]; dudy1 = C_U_G(c1,t)[1]; dudz1 = C_U_G(c1,t)[2]; dvdx1 = C_V_G(c1,t)[0]; dvdy1 = C_V_G(c1,t)[1]; dvdz1 = C_V_G(c1,t)[2]; dwdx1 = C_W_G(c1,t)[0]; dwdy1 = C_W_G(c1,t)[1]; dwdz1 = C_W_G(c1,t)[2]; SRT1[0][0] = 0.5*(dudx1 + dudx1); SRT1[0][1] = 0.5*(dudy1 + dvdx1); SRT1[0][2] = 0.5*(dudz1 + dwdx1); SRT1[1][0] = 0.5*(dvdx1 + dudy1); SRT1[1][1] = 0.5*(dvdy1 + dvdy1); SRT1[1][2] = 0.5*(dvdz1 + dwdy1); SRT1[2][0] = 0.5*(dwdx1 + dudz1); SRT1[2][1] = 0.5*(dwdy1 + dvdz1); SRT1[2][2] = 0.5*(dwdz1 + dwdz1); /* Strain rate tensor on the face */ Sface[0][0] = 0.5*(SRT[0][0] + SRT1[0][0]); Sface[0][1] = 0.5*(SRT[0][1] + SRT1[0][1]); Sface[0][2] = 0.5*(SRT[0][2] + SRT1[0][2]); Sface[1][0] = 0.5*(SRT[1][0] + SRT1[1][0]); Sface[1][1] = 0.5*(SRT[1][1] + SRT1[1][1]); Sface[1][2] = 0.5*(SRT[1][2] + SRT1[1][2]); Sface[2][0] = 0.5*(SRT[2][0] + SRT1[2][0]); Sface[2][1] = 0.5*(SRT[2][1] + SRT1[2][1]); Sface[2][2] = 0.5*(SRT[2][2] + SRT1[2][2]); } else { if (NNULLP(THREAD_STORAGE(tf,SV_DENSITY))) dens = F_R(f,tf); /* Set dens to face value if available */ else dens = C_R(c,t); /* else, set dens to cell value */ Sface[0][0] = 0.0; Sface[0][1] = 0.0; Sface[0][2] = 0.0; Sface[1][0] = 0.0; Sface[1][1] = 0.0; Sface[1][2] = 0.0; Sface[2][0] = 0.0; Sface[2][1] = 0.0; Sface[2][2] = 0.0; } /* Flux through the face */ flux = fabs(F_FLUX(f,tf)); /* Numerator of the formula for some other formula */ for(ii=0; ii<3; ii++) for(jj=0; jj<3; jj++) { DS[ii][jj] += Sface[ii][jj]*flux/dens; } } /* Here ends the face loop */ } end_c_loop(c,t) #endif } 

July 21, 2013, 21:26 

#2 
Senior Member
Join Date: Aug 2011
Posts: 315
Rep Power: 13 
In my opinion you should calculate the inner product of the strainrate tensor (Sface) and the flux vector, i.e., F_FLUX()*F_AREA(), rather than a scalar F_FLUX times the tensor Sface. For the record, C0 and C1 are well defined and can be unambiguously distinguished.
I do not quite understand your second question, but F_R and C_R are the same thing. So if the storage is not allocated you should receive the segmentation fault error or things alike. 

July 23, 2013, 05:42 

#3 
New Member
Join Date: Jun 2013
Posts: 15
Rep Power: 5 
Thank you for the answer, blackmask.
What you mean is something like this: First assign flux_vector  vector F_FLUX()*F_AREA() and then Code:
for(ii=0; ii<3; ii++) for(jj=0; jj<3; jj++) DS[ii][jj] += Sface[ii][jj]*flux_vector[jj] Good to know that C0 and C1 are well defined. I misunderstood some other thread when I was searching for the similiar topic. Thanks for clreaing it up. As for the second question: I copied the part with C_R and F_R from UDF manual, 2.7.3. DEFINE_UDS_FLUX example. As I understand, it is supposed to check whether density is defined on the face and read this face density (F_R) or density from cell centre (C_R) if otherwise. As for the periodic boundary condition stuff. If you look at the DSfabsflux thumbnail, strange (unphysical) drop of DS values appears at the very periodic bc, that is on the left, right and middle vertical lines (I didn't state it explicitly, but the pictures are /Views/Define Periodic Repeats twice, so oly half is computed actually). However, when I inspect present results I don't see this drop anymore. I don't really know why it behaves like this (I don't recall any changes that I have done), but it seems like the question is gone. Hope this makes some sense now. Regards, kornetka 

July 23, 2013, 20:11 

#4 
Senior Member
Join Date: Aug 2011
Posts: 315
Rep Power: 13 
Sorry, You are right on the second question. While "F_R" and "C_R" are identical, I did not notice that "f" and "tf" in
Code:
F_R(f,tf) Code:
C_R(c,t) Code:
for(ii=0; ii<3; ii++) for(jj=0; jj<3; jj++) DS[ii] += Sface[ii][jj]*flux_vector[jj] DS[ii][jj][kk] = Sface[ii][jj]*flux_vector[kk] Anyway, you must be better than me in this topic and you have a good explanation for the DS you calculated. 

July 25, 2013, 06:54 

#5 
New Member
Join Date: Jun 2013
Posts: 15
Rep Power: 5 
Of course, you're right, my tensorvector multiplication was a total brainfart.
What I'm trying to do is calculate the substantial derivative of the strain rate tensor. MenterSmirnov article on which I base my attempts explicitly states following formula for it: DS/Dt =(SumOverFaces( StrainRateTensor@face * VelocityNormal@face * face_area))/volume I'm trying to compute the Vel_normal as fluxthroughface/(face_area*density), as also stated in the source article, but this simple implementation leads to the different sign of the DS in some region, as pictured in my first post. During the conversation with my supervisor we came to the conclusion that it is most probably just the sign of the flux, and if it can be assigned properly, results are going to be fine. So for now I'm going to focus on this. Thank you for all your help, blackmask. Regards, kornetka 

Tags 
mesh, partition, periodic bc, udf 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
No layers in a small gap  bobburnquist  OpenFOAM Native Meshers: snappyHexMesh and Others  6  August 26, 2015 09:38 
CFX Mesh Deformation problem  Silmaril  CFX  7  October 19, 2010 10:00 
Dynamic Mesh Problem.  Tom Clark  FLUENT  9  July 7, 2010 07:56 
Dynamic Mesh UDF  Qureshi  FLUENT  1  December 2, 2009 01:27 
mesh problem and looking for udf sine wave'case  wanghong  FLUENT  0  March 22, 2005 20:43 