When trying to run some simulations with gravity effects in 2D models with 1D fracture elements I stumbled across the fact that the velocities in the fractures where always extremely (read m/sec) high even if the system was completely hydrostatic, i.e. no fluid should move up or down. I included the correction for 1D fracture velocities under the influence of gravity as programmed in the CSMP reservoir simulator, which made things even worse although velocities were nicely aligned with the fracture (BTW, I don’t quite follow the logic what is done here). I then went on to reduce the problem to a homogeneous matrix with a single oblique fracture, which still gave ridiculously high velocities in the fracture at a perfectly homogeneous and uniform hydrostatic pressure where no fluid must move up or down. Further invistigations showed that there is something odd with the gradient calculations: The one I obtain using the interpolation function derivatives differ, sometimes a factor 10, from the ones I would calculate by hand.

I stripped down the code to the following stand-along test:

stl_index DIM(2); IsoparametricLinearLineElement elm(DIM,1); Element<csp_float,2> e(&elm); Node<csp_float,2> n0, n1; n0.x(0.0); n0.y(0.0); n1.x(1.0); n1.y(1.0); e.ConnectTo( 0, n0 ); e.ConnectTo( 1, n1 ); csp_float l = e.Volume(); csp_float p[2], v[2]; p[0] = 0.0; p[1] = 1.0; v[0] = v[1] = 0.0; DenseMatrix<csp_float,DM_MIN> DN, XY; e.CoordinateMatrix(XY); XY.Out(); for ( stl_index i=0; i<e.IntegrationPoints(); i++ ) { e.dN_AtIntegrationPoint( DN, i, 1 ); for ( stl_index n=0; n<e.Nodes(); n++ ) for ( stl_index j=0; j<DIM; j++ ) v[j] += p[n] * DN(j,n); } cout << "\n\nFE Gradient X: " << v[0] << ", should be " << (p[1]-p[0])/(n1.x()-n0.x()); cout << "\n\nFE Gradient Y: " << v[1] << ", should be " << (p[1]-p[0])/(n1.y()-n0.y()); cout.flush();

This is a 1D line element with two points at (0,0) and (1,1). There is a pressure of 0 at node 0 and a pressure of 1 at node 1. The length of the element is correctly calculated as square root of 2 and the pressure gradient in X and Y direction must be 1 (with the length of the vector between point 0 and 1 is square root of 2). Now, this is not what using the interpolation function derivatives yield. They yield gradients of 0.5 in X and Y directions. If I use coordinates of the nodes and the pressures from the actual mesh, I can recomputed the gradients that the “real” CSMP gives me, which differ from what I would think are the right ones. It only works if the line element is oriented horizontally or vertically.

Now, aside from the gravity effects, this could be very worrying because it implies that all fluid velocities in 2D DFM models where the fractures are 1D line elements are potentially not correct.

Am I doing something entirely stupid or is there an issue with the 1D line elements in 2D space?

Cheers, Sebastian