CFD Online Logo CFD Online URL
Home > Forums > Fluent UDF and Scheme Programming

Mesh UDF problem

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By blackmask

LinkBack Thread Tools Display Modes
Old   July 21, 2013, 10:15
Default Mesh UDF problem
New Member
kornetka's Avatar
Join Date: Jun 2013
Posts: 15
Rep Power: 5
kornetka is on a distinguished road
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 DS-Partition 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 fabs-ed the flux. Resulting distribution is shown in DS-fabsflux. 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 DS-fabsflux. 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.99873e-01
- Maximum Aspect Ratio = 3.52643e+01


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) 
            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);

                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]);
                        if (NNULLP(THREAD_STORAGE(tf,SV_DENSITY)))
                            dens = F_R(f,tf); /* Set dens to face value if available */
                            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    */
Attached Images
File Type: png Partitions.png (18.8 KB, 12 views)
File Type: jpg DS-Partitions.jpg (62.9 KB, 11 views)
File Type: png DS-fabsflux.png (22.9 KB, 14 views)
File Type: png Mesh.png (93.9 KB, 10 views)
kornetka is offline   Reply With Quote

Old   July 21, 2013, 21:26
Senior Member
Join Date: Aug 2011
Posts: 315
Rep Power: 13
blackmask will become famous soon enough
In my opinion you should calculate the inner product of the strain-rate 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.
kornetka likes this.
blackmask is offline   Reply With Quote

Old   July 23, 2013, 05:42
New Member
kornetka's Avatar
Join Date: Jun 2013
Posts: 15
Rep Power: 5
kornetka is on a distinguished road
Thank you for the answer, blackmask.

What you mean is something like this:
First assign flux_vector - vector F_FLUX()*F_AREA() and then
for(ii=0; ii<3; ii++)
  for(jj=0; jj<3; jj++)        
     DS[ii][jj] += Sface[ii][jj]*flux_vector[jj]
In other words DS[3x3] = S_face[3x3]*Flux_vector[1x3], right? Your advice seems to solve this problem.
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 DS-fabsflux 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
kornetka is offline   Reply With Quote

Old   July 23, 2013, 20:11
Senior Member
Join Date: Aug 2011
Posts: 315
Rep Power: 13
blackmask will become famous soon enough
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
has changed to "c" and "t" in
When a tensor contract with a vector, the result should be a vector. So it should be
	for(ii=0; ii<3; ii++)
 for(jj=0; jj<3; jj++)              DS[ii] += Sface[ii][jj]*flux_vector[jj]
I find it hard to understand the physical interpretation for DS in your original form but it does not have much meaning in my revised form. If you want to study the convection term of the Reynolds stress term (or something a like) then DS should be a third-order tensor, with

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.
blackmask is offline   Reply With Quote

Old   July 25, 2013, 06:54
New Member
kornetka's Avatar
Join Date: Jun 2013
Posts: 15
Rep Power: 5
kornetka is on a distinguished road
Of course, you're right, my tensor-vector multiplication was a total brainfart.
What I'm trying to do is calculate the substantial derivative of the strain rate tensor. Menter-Smirnov 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 flux-through-face/(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
kornetka is offline   Reply With Quote


mesh, partition, periodic bc, udf

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On

Similar Threads
Thread Thread Starter Forum Replies Last Post
Dynamic Mesh UDF Qureshi FLUENT 7 March 23, 2017 08:37
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
mesh problem and looking for udf sine wave'case wanghong FLUENT 0 March 22, 2005 20:43

All times are GMT -4. The time now is 00:21.