User Defines power-law viscosity model

MFiX
#1

Hello everyone.
Now I’m trying to simulate a non-Newtonian liquid-solid fluidized bed with 2D model, and power-law fluids are selected. And now I wanna replace constant viscosity with power-law viscosity, and I tried to put constitutive equation of power-law fluid into USR_PROP_Mug in usr-properties.f, but unfortunately it didn’t work.The constitutive equation of power-law fluid is as follows:
1638339549(1)
and for D:D, I used the codes in calc_mu_s.f
1638339812(1)
I hope someone to help me make the codes right!
By the way, I also want to know why solid phase is double defined as M, MM.usr_properties.f (38.4 KB)

#2

Hi Zihan. Please provide more details, you say this “didn’t work” but how did it fail? If you can attach all files, we can try running your case here. Use “Submit bug report” from the main menu to create a zipfile containing all files needed to run the case.

If you look at where the MM variable is used, it’s only used as a loop index:

                  DO MM = 1, SMAX
                     SUM_EpsGo =  SUM_EpsGo+EP_s(IJK,MM)*G_0(IJK,M,MM)
                  ENDDO

MM is a separate variable from M so we can use it as a loop index without changing the value of M

– Charles

#3

Thank you, Charles. Now I have uploaded all files. It can be built successfully, however it won’t run and calculate. I’ll appreciate it and look forward to your help.
MN.mfx (18.0 KB)
usr_properties.f (38.3 KB)

#4

When I run this I get an immediate error:

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  0x7fd6b00e66b0 in ???
#1  0x7fd6b00e5865 in ???
#2  0x7fd6fcc1931f in ???
#3  0x7fd6b0428c3e in usr_prop_mug_
	at /tmp/MN/usr_properties.f:235

This is an issue with your code, not with MFiX.

Line 235 of usr_properties.f looks like this:

                  Mu_sv =((2.d0+ALPHA)/3.d0)*((Mu_star/(Eta*(2.d0-Eta)*&
                     G_0(IJK,M,M)))*(ONE+1.6d0*Eta*SUM_EpsGo)&
                     *(ONE+1.6d0*Eta*(3d0*Eta-2d0)*&
                     SUM_EpsGo)+(0.6d0*Mu_bv*Eta))

There’s a lot of calculation going on there, and somewhere you are dividing by zero.

I suggest breaking up the complex calculation into several sub-calculations, and checking all denominators before doing any division, as described in this forum posting

#5

Okay, thanks for your kind apply, Charles. I’ll check something about it.