[kornetka]

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

Regards,
kornetka

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
}

[blackmask]

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]

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

[blackmask]

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

F_R(f,tf)

has changed to "c" and "t" in

C_R(c,t)

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.

[kornetka]

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