Finite Analytic Method with Anisotropy

The permeabilities are expressed in matrix form, as can be seen in Figure 5.1 where the permeabilities

are shown with the tensor components.

Finite Analytic Method Extension: Anisotropic Solution

29 | P a g e

Figure 5.1: Permeability tensor components shown alongside their respective permeability plug and the features they represent (Nelson, 2001)

Starting with the full symmetric tensor representation of the permeability

Finite Analytic Method Extension: Anisotropic Solution

30 | P a g e 𝑘

_𝑖

= ( 𝑘

_𝑖11

𝑘

_𝑖12

𝑘

_𝑖21

𝑘

_𝑖22

)

Equation 41

Where 𝑖 refers to the quadrants, and 𝑘

_𝑖12

and 𝑘

_𝑖21

can be set equal to each other due to orientation, as is seen in Figure 5.1.

The eigenvalues of the permeability tensor 𝑘

_𝑖

can then be calculated from the following expression 𝑑𝑒𝑡(𝜆𝐼 − 𝑘

_𝑖

) = | 𝜆 − 𝑘

_𝑖11

𝑘

_𝑖12

𝑘

_𝑖12

𝜆 − 𝑘

_𝑖22

| = (𝜆 − 𝑘

𝑖₁₁

)(𝜆 − 𝑘

_𝑖22

) − 𝑘

_𝑖12²

= 0 𝜆

²

− (𝑘

𝑖₁₁

+ 𝑘

_𝑖22

)𝜆 + 𝑘

_𝑖11

𝑘

_𝑖22

− 𝑘

_𝑖12²

= 0

𝜆 = 1

2 [(𝑘

_𝑖11

+ 𝑘

_𝑖22

) ± √(𝑘

_𝑖11

+ 𝑘

_𝑖22

)

²

− 4(𝑘

_𝑖11

𝑘

_𝑖22

− 𝑘

_𝑖12²

)]

𝜆 = 1

2 [(𝑘

𝑖₁₁

+ 𝑘

_𝑖22

) ± √𝑘

_𝑖11²

+ 𝑘

_𝑖22²

− 2𝑘

_𝑖11

𝑘

_𝑖22

+ 4𝑘

_𝑖12²

]

Equation 42

For a given eigenvector, the following applies:

( 𝜆 − 𝑘

_𝑖11

𝑘

_𝑖12

𝑘

_𝑖12

𝜆 − 𝑘

_𝑖22

) ( 𝑥

𝑦) = ( 0 0 ) 𝑥(𝜆 − 𝑘

_𝑖11

) − 𝑦𝑘

_𝑖12

= 0 𝑦(𝜆 − 𝑘

𝑖₂₂

) − 𝑥𝑘

_𝑖12

= 0

Equation 43

Therefore, using (𝜆 − 𝑘

𝑖₁₁

)(𝜆 − 𝑘

_𝑖22

) = 𝑘

_𝑖12²

from Equation 42, the two eigen vectors can then be shown to be:

( 𝑥 𝑦) = (

𝜆

₁

− 𝑘

_𝑖22

𝑘

_𝑖12

) & ( 𝜆

₂

− 𝑘

_𝑖22

𝑘

_𝑖12

)

Equation 44

This is well known from linear algebra (Anton, 2010).

A similarity transformation can now be performed on the matrix 𝑘

_𝑖

to obtain a diagonal matrix 𝑘̃

_𝑖

, where 𝑘̃

_𝑖

= 𝑄

⁻¹

𝑘

_𝑖

𝑄 , 𝑄

^𝜏

and 𝑄 is a matrix consisting of the eigenvectors for 𝑘

_𝑖

.

𝑘̃

_𝑖

= ( 𝑘̃

_𝑖1

0 0 𝑘̃

_𝑖2

)

Equation 45

Finite Analytic Method Extension: Anisotropic Solution

31 | P a g e Furthermore, it is well known that this similarity transformation corresponds to a rotation of the original standard Cartesian coordinate system to a coordinate system where the new axis are aligned with the eigenvectors of 𝑘

_𝑖

. The rotation, which is denoted as 𝜃̃

_𝑖

, is given by the expression:

Returning to the pressure equation with the full tensor permeability, it can now be written as:

∇ ∙ [(𝑘

_𝑖

(𝑥, 𝑦)∇𝑃

_𝑖

)] = 𝑘

_𝑖11

𝜕

²

𝑃

_𝑖

𝜕𝑥

²

+ 2𝑘

_𝑖12

𝜕

²

𝑃

_𝑖

𝜕𝑥𝜕𝑦 + 𝑘

_𝑖22

𝜕

²

𝑃

_𝑖

𝜕𝑦

²

= 0

Equation 47

Furthermore, it is also noted that:

∇ ∙ [(𝑘

_𝑖

(𝑥, 𝑦)∇𝑃

_𝑖

)] ∙ 𝑛⃗ =

𝑘̃

_𝑖

= 𝑄

⁻¹

𝑘

_𝑖

𝑄 shows how with a repositioned coordinate system, one can eliminate the off diagonal elements. In the rotated coordinate system, denoted by 𝑥̃, 𝑦̃, Equation 47 can now be written as:

𝑘̃

_𝑖1

𝜕

²

𝑃̃

_𝑖

Such that there is a rotation 𝜃̃

_𝑖

between the 𝑥, 𝑦 coordinate system and the 𝑥̃, 𝑦̃ coordinate system.

Moreover, an additional transformation is introduced (stretching transformation), such that:

Finite Analytic Method Extension: Anisotropic Solution

32 | P a g e 𝑥̃ = 𝑥̂√𝑘̃

_𝑖1

= 𝑥 cos 𝜃̃

_𝑖

+ 𝑦 sin 𝜃̃

_𝑖

→ 𝜕𝑥̂

𝜕𝑥 = cos 𝜃̃

_𝑖

√𝑘̃

_𝑖1

𝑦̃ = 𝑦̂√𝑘̃

_𝑖2

= 𝑦 cos 𝜃̃

_𝑖

− 𝑥 sin 𝜃̃

_𝑖

→ 𝜕𝑦̂

𝜕𝑦 = cos 𝜃̃

_𝑖

√𝑘̃

_𝑖2

Equation 50

In the new 𝑥̂, 𝑦̂ coordinate system, the pressure equation can finally be written as:

𝜕

²

𝑃

_𝑖

𝜕𝑥̂

²

+ 𝜕

²

𝑃

₂

𝜕𝑦̂

²

= 0

Equation 51

In other words, the Laplace equation for the pressure can is obtained in the 𝑥̂, 𝑦̂ coordinate system.

Introducing the variable 𝜏

_𝑖

, which represents the rotation of the axis and is described in complex form as 𝜏

_𝑖

= 𝑒

^𝑖𝜃̃

𝑧 = 𝑥 + 𝑖𝑦 , 𝑧̅ = 𝑥 − 𝑖𝑦 𝜏𝑧 = 𝑧̂ , 𝜏𝑧 ̅ = 𝑧̂̅

𝑧̂ = 𝑥̂ + 𝑖𝑦̂ , 𝑧̂̅ = 𝑥̂ − 𝑖𝑦̂

𝑧̂ = 𝑥 cos 𝜃̃

_𝑖

+ 𝑦 sin 𝜃̃

_𝑖

√𝑘̃

_𝑖1

+ 𝑖 𝑦 cos 𝜃̃

_𝑖

− 𝑥 sin 𝜃̃

_𝑖

√𝑘̃

_𝑖2

𝑧̂̅ = 𝑥 cos 𝜃̃

_𝑖

+ 𝑦 sin 𝜃̃

_𝑖

√𝑘̃

_𝑖1

− 𝑖 𝑦 cos 𝜃̃

_𝑖

− 𝑥 sin 𝜃̃

_𝑖

√𝑘̃

_𝑖2

Equation 52

By following to a certain degree the same methodology as with the isotropic case, the analytical nodal solution for the anisotropic case can be determined. However, the computations in the following will naturally be much more involved. The solutions for the individual pressure equation will be using the coordinate system described above (𝑥̂, 𝑦̂).

Assuming the same power law behavior as described in Equation 17, (but in the 𝑥̂, 𝑦̂ coordinate system) the pressure equations can then be expressed as follows:

𝑃

_𝑖

= 𝑓

_𝑖

(𝜏𝑧) + 𝑔

_𝑖

(𝜏𝑧 ̅ ) , 𝑜𝑟 𝑃

_𝑖

(𝑥̂, 𝑦̂) = 𝑓

_𝑖

(𝑧̂) + 𝑔

_𝑖

(𝑧̂̅)

Equation 53

As before, 𝑓

_𝑖

and 𝑔

_𝑖

are complex analytic functions. Thus, for the first quadrant, the following is valid:

Finite Analytic Method Extension: Anisotropic Solution

33 | P a g e 𝑓

₁

(𝑥̂) + 𝑔

₁

(𝑥̂) = 𝐴

₁

𝑥̂

^1−𝛼1

𝑖𝑓

^′₁

(𝑥̂) − 𝑖𝑔

^′₁

(𝑥̂) = (1 − 𝛽

₁

)𝐵

₁

𝑥̂

^−𝛽1

Equation 54

The details of the solution is presented in the following:

𝑓

^′₁

(𝑥̂) + 𝑔

^′₁

(𝑥̂) = (1 − 𝛼

₁

)𝐴

₁

𝑥̂

^−𝛼1

, |× 𝑖 𝑎𝑛𝑑 𝑠𝑢𝑚 2𝑖𝑓

^′₁

(𝑥̂) = (1 − 𝛼

₁

)𝑖𝐴

₁

𝑥̂

^−𝛼1

+ (1 − 𝛽

₁

)𝐵

₁

𝑥̂

^−𝛽1

, |× 1

2𝑖 𝑓

^′₁

(𝑥̂) = 1

2 {(1 − 𝛼

₁

)𝐴

₁

𝑥̂

^−𝛼1

− (1 − 𝛽

₁

)𝑖𝐵

₁

𝑥̂

^−𝛽1

} 𝑡ℎ𝑢𝑠: 𝑓

₁

(𝑥̂) = 1

2 (𝐴

₁

𝑥̂

^1−𝛼1

− 𝑖𝐵

₁

𝑥̂

^−𝛽1

+ 𝐶) 𝑔

₁

(𝑥̂) = 1

2 (𝐴

₁

𝑥̂

^1−𝛼1

+ 𝑖𝐵

₁

𝑥̂

^−𝛽1

− 𝐶) 𝑃

₁

(𝑥̂, 𝑦̂) = 1

2 {𝐴

₁

[𝑧̂

^1−𝛼1

+ 𝑧̂̅

^1−𝛼1

] − 𝑖𝐵

₁

[𝑧̂

^1−𝛽1

− 𝑧̂̅

^1−𝛽1

]}

Equation 55

Performing the same procedure on all four quadrants, the pressure equations now read:

𝑃

₁

(𝑥̂, 𝑦̂) = 1

2 {𝐴

₁

[(𝑧̂)

^1−𝛼1

+ (𝑧̂̅)

^1−𝛼1

] − 𝑖𝐵

₁

[(𝑧̂)

^1−𝛽1

− (𝑧̂̅)

^1−𝛽1

]}

𝑃

₂

(𝑥̂, 𝑦̂) = 1

2 {𝐴

₂

[(−𝑖𝑧̂)

^1−𝛼2

+ (𝑖𝑧̂̅)

^1−𝛼2

] + 𝑖𝐵

₂

[(−𝑖𝑧̂)

^1−𝛽2

− (𝑖𝑧̂̅)

^1−𝛽2

]}

𝑃

₃

(𝑥̂, 𝑦̂) = 1

2 {𝐴

₃

[(−𝑧̂)

^1−𝛼3

+ (−𝑧̂̅)

^1−𝛼3

] + 𝑖𝐵

₃

[(−𝑧̂)

^1−𝛽3

− (−𝑧̂̅)

^1−𝛽3

]}

𝑃

₄

(𝑥̂, 𝑦̂) = 1

2 {𝐴

₄

[(𝑖𝑧̂)

^1−𝛼4

+ (−𝑖𝑧̂̅)

^1−𝛼4

] − 𝑖𝐵

₄

[(𝑖𝑧̂)

^1−𝛽4

− (−𝑖𝑧̂̅)

^1−𝛽4

]}

Equation 56

Their respective derivatives are then deduced to be:

Finite Analytic Method Extension: Anisotropic Solution

34 | P a g e

𝜕𝑃

₁

𝜕𝑥̂ = 1

2 {(1 − 𝛼

₁

)𝐴

₁

[(𝑧̂)

^−𝛼1

+ (𝑧̂̅)

^−𝛼1

] − (1 − 𝛽

₁

)𝑖𝐵

₁

[(𝑧̂)

^−𝛽1

− (𝑧̂̅)

^−𝛽1

]}

𝜕𝑃

₁

𝜕𝑦̂ = 1

2 𝑖 {(1 − 𝛼

₁

)𝐴

₁

[(𝑧̂)

^−𝛼1

− (𝑧̂̅)

^−𝛼1

] − (1 − 𝛽

₁

)𝑖𝐵

₁

[(𝑧̂)

^−𝛽1

+ (𝑧̂̅)

^−𝛽1

]}

𝜕𝑃

₂

𝜕𝑥̂ = 1

2 {(1 − 𝛼

₂

)𝐴

₂

[(−𝑖𝑧̂)

^−𝛼2

+ (𝑖𝑧̂̅)

^−𝛼2

] + (1 − 𝛽

₂

)𝑖𝐵

₂

[(−𝑖𝑧̂)

^−𝛽2

− (𝑖𝑧̂̅)

^−𝛽2

]}

𝜕𝑃

₂

𝜕𝑦̂ = 1

2 𝑖 {(1 − 𝛼

2

)𝐴

₂

[(−𝑖𝑧̂)

^−𝛼2

− (𝑖𝑧̂̅)

^−𝛼2

] + (1 − 𝛽

₂

)𝑖𝐵

₂

[(−𝑖𝑧̂)

^−𝛽2

+ (𝑖𝑧̂̅)

^−𝛽2

]}

𝜕𝑃

₃

𝜕𝑥̂ = 1

2 {(1 − 𝛼

₃

)𝐴

₃

[(−𝑧̂)

^−𝛼3

+ (−𝑧̂̅)

^−𝛼3

] + (1 − 𝛽

₃

)𝑖𝐵

₃

[(−𝑧̂)

^−𝛽3

− (−𝑧̂̅)

^−𝛽3

]}

𝜕𝑃

₃

𝜕𝑦̂ = 1

2 𝑖 {(1 − 𝛼

₃

)𝐴

₃

[(−𝑧̂)

^−𝛼3

− (−𝑧̂̅)

^−𝛼3

] + (1 − 𝛽

₃

)𝑖𝐵

₃

[(−𝑧̂)

^−𝛽3

+ (−𝑧̂̅)

^−𝛽3

]}

𝜕𝑃

₄

𝜕𝑥̂ = 1

2 {(1 − 𝛼

₄

)𝐴

₄

[(𝑖𝑧̂)

^−𝛼4

+ (−𝑖𝑧̂̅)

^−𝛼4

] − (1 − 𝛽

₄

)𝑖𝐵

₄

[(𝑖𝑧̂)

^−𝛽4

− (−𝑖𝑧̂̅)

^−𝛽4

]}

𝜕𝑃

₄

𝜕𝑦̂ = 1

2 𝑖 {(1 − 𝛼

4

)𝐴

₄

[(𝑖𝑧̂)

^−𝛼4

− (−𝑖𝑧̂̅)

^−𝛼4

] − (1 − 𝛽

₄

)𝑖𝐵

₄

[(𝑖𝑧̂)

^−𝛽4

+ (−𝑖𝑧̂̅)

^−𝛽4

]}

Equation 57

For the same reason as mentioned in the

Finite Analytic Method section; 𝛼

_𝑖

= 𝛽

_𝑖

= 𝛼 . This is due to the fact that if unique solutions are to be possible, this has to occur.

From Equation 48, the various flux expressions for all the quadrants are as follows:

Finite Analytic Method Extension: Anisotropic Solution

The final pressure equations must abide by the same criteria as was set in Equation 17, however, using

the new coordinate system. Thus, the criteria can be expressed as follows:

Finite Analytic Method Extension: Anisotropic Solution

36 | P a g e 1. Along the positive x-axis:

𝑃

₁

(𝑥̂, 𝑦̂) = 𝑃

₄

(𝑥̂, 𝑦̂) 2. Along the positive y-axis:

𝑃

₂

(𝑥̂, 𝑦̂) = 𝑃

₁

(𝑥̂, 𝑦̂) 3. Along the negative x-axis:

𝑃

₃

(𝑥̂, 𝑦̂) = 𝑃

₂

(𝑥̂, 𝑦̂) 4. Along the negative y-axis:

𝑃

₄

(𝑥̂, 𝑦̂) = 𝑃

₃

(𝑥̂, 𝑦̂)

Using the polar form of complex numbers, the pressure equations set up in Equation 56 can be written as:

With this, the analytical nodal solution for the anisotropic case is expressed. However, the unknowns 𝐴

_𝑖

, 𝐵

_𝑖

, 𝛼 are yet to be determined. This is delayed till after the numerical fragment is examined (see end of section 3.3). Below is the necessary approach to achieve the correct numerical expressions for the fluxes needed in the implementation of the anisotropic extension for the finite analytic method.

Similar to before, the fluxes need to be expressed at the first and third quadrant, and to be included in the final solution for the inter nodal transmissivity. This is done to uphold the continuity critera described above. Taking basis in Equation 58, the fluxes can be denoted as such:

Note: in following equations; ∆𝑥 =

¹

2

∆𝑥 and ∆𝑦 =

¹

2

∆𝑦. This was done again to make the transition into

the program easier as the foundation of the code that was used in this study applied this notation.

Finite Analytic Method Extension: Anisotropic Solution

Using the following facts:

𝜕𝑥̂

Finite Analytic Method Extension: Anisotropic Solution

Finite Analytic Method Extension: Anisotropic Solution

39 | P a g e And thus, the final solution for the inter nodal transmissivity can be expressed as such:

𝛾

_𝑥1

= 𝜆̂

_𝑥1

𝜆̂

₂

− 𝜆̂

₁

, 𝛾

_𝑦1

= 𝜆̂

_𝑦1

𝜆̂

₄

− 𝜆̂

₁

𝛾

_𝑥2

= 𝜆̂

_𝑥1

𝜆̂

₂

− 𝜆̂

₁

, 𝛾

_𝑦2

= 𝜆̂

_𝑦3

𝜆̂

₃

− 𝜆̂

₂

𝛾

_𝑥3

= 𝜆̂

_𝑥3

𝜆̂

₃

− 𝜆̂

₄

, 𝛾

_𝑦3

= 𝜆̂

_𝑦3

𝜆̂

₃

− 𝜆̂

₂

𝛾

_𝑥4

= 𝜆̂

_𝑥3

𝜆̂

₃

− 𝜆̂

₄

, 𝛾

_𝑦4

= 𝜆̂

_𝑦1

𝜆̂

₄

− 𝜆̂

₁

Equation 63

Lastly, to solve for the unknowns 𝐴

_𝑖

and 𝐵

_𝑖

, which are functions of 𝛼, a backward substitution needs to be applied similar to that described in the previous section. 𝐴

_𝑖

and 𝐵

_𝑖

form a cyclical behavior (see Equation 66), where the solution of one is determined by the other. As such, by expressing equations for these unknowns, it is not only possible to set up a solution for each, but also solve for 𝛼 .

Seeing as how these sets of equations are much more involved than before, several temporary variables and changes in notation need to be established so as to keep the expressions more contained.

Starting by first expressing the following quantities with these notations:

𝑥̂

_𝑖

= 𝑥 cos 𝜃̃

_𝑖

+ 𝑦 sin 𝜃̃

_𝑖

√𝑘̃

_𝑖1

, 𝑦̂

_𝑖

= 𝑦 cos 𝜃̃

_𝑖

− 𝑥 sin 𝜃̃

_𝑖

√𝑘̃

_𝑖2

𝜕𝑥̂

_𝑖

𝜕𝑥 = cos 𝜃̃

_𝑖

√𝑘̃

_𝑖1

, 𝜕𝑦̂

_𝑖

𝜕𝑦 = cos 𝜃̃

_𝑖

√𝑘̃

_𝑖2

𝑟̂

_𝑖

= √𝑥̂

_𝑖²

+ 𝑦̂

_𝑖²

, 𝜃̂

_𝑖

= 𝜃 − 𝜃̃

_𝑖

𝜕𝑟̂

_𝑖

𝜕𝑥̂

_𝑖

= 𝑥̂

_𝑖

𝑟̂

_𝑖

, 𝜕𝑟̂

_𝑖

𝜕𝑦̂

_𝑖

= 𝑦̂

_𝑖

𝑟̂

_𝑖

, 𝜕𝜃̂

_𝑖

𝜕𝑥̂

_𝑖

= − 𝑦̂

_𝑖

𝑟̂

_𝑖²

, 𝜕𝜃̂

_𝑖

𝜕𝑦̂

_𝑖

= 𝑥̂

_𝑖

𝑟̂

_𝑖²

Equation 64

Where 𝜃 is the angle between the vertex and the grid node in the standard 𝑥, 𝑦 Cartesian coordinate

system. Moreover, the following set of definitions are introduced:

Finite Analytic Method Extension: Anisotropic Solution

Finite Analytic Method Extension: Anisotropic Solution

Finite Analytic Method Extension: Anisotropic Solution

Finite Analytic Method Extension: Anisotropic Solution

Finite Analytic Method Extension: Anisotropic Solution

44 | P a g e 𝑎

₃₂

= 𝑘

₃₁₂

cos 𝜃̃

₃

√𝑘̃

₃₁

(1 − 𝛼) { 𝜕𝑟̂

₃

𝜕𝑥̂

₃

𝑟̂

₃^−𝛼

cos[(𝜃̂

₃

+ 𝜋)(1 − 𝛼)]

− 𝜕𝜃̂

₃

𝜕𝑥̂

₃

𝑟̂

₃^1−𝛼

sin[(𝜃̂

3

+ 𝜋)(1 − 𝛼)]}

+ 𝑘

₃₂₂

cos 𝜃̃

₃

√𝑘̃

₃₂

(1 − 𝛼) { 𝜕𝑟̂

₃

𝜕𝑦̂

₃

𝑟̂

₃^−𝛼

cos[(𝜃̂

₃

+ 𝜋)(1 − 𝛼)]

− 𝜕𝜃̂

₃

𝜕𝑦̂

₃

𝑟̂

₃^1−𝛼

sin[(𝜃̂

₃

+ 𝜋)(1 − 𝛼)]}

𝑏

₃₂

= −𝑘

₃₁₂

cos 𝜃̃

₃

√𝑘̃

₃₁

(1 − 𝛼) { 𝜕𝑟̂

₃

𝜕𝑥̂

₃

𝑟̂

₃^−𝛼

sin[(𝜃̂

₃

+ 𝜋)(1 − 𝛼)]

+ 𝜕𝜃̂

₃

𝜕𝑥̂

₃

𝑟̂

₃^1−𝛼

cos[(𝜃̂

3

+ 𝜋)(1 − 𝛼)]}

− 𝑘

₃₂₂

cos 𝜃̃

₃

√𝑘̃

₃₂

(1 − 𝛼) { 𝜕𝑟̂

₃

𝜕𝑦̂

₃

𝑟̂

₃^−𝛼

sin[(𝜃̂

₃

+ 𝜋)(1 − 𝛼)]

+ 𝜕𝜃̂

₃

𝜕𝑦̂

₃

𝑟̂

₃^1−𝛼

cos[(𝜃̂

₃

+ 𝜋)(1 − 𝛼)]}

Equation 65

With this simplification, both the pressure continuity and the flux continuity can be systematically written as:

1. Along the positive x-axis (𝑦 = 0):

𝐴

₄

𝑎

₄₁

+ 𝐵

₄

𝑏

₄₁

= 𝐴

₁

𝑎

₁₁

+ 𝐵

₁

𝑏

₁₁

𝐴

₄

𝑎

₄₂

+ 𝐵

₄

𝑏

₄₂

= 𝐴

₁

𝑎

₁₂

+ 𝐵

₁

𝑏

₁₂

2. Along the positive y-axis (𝑥 = 0):

𝐴

₁

𝑎

₁₁

+ 𝐵

₁

𝑏

₁₁

= 𝐴

₂

𝑎

₂₁

+ 𝐵

₂

𝑏

₂₁

𝐴

₁

𝑎

₁₂

+ 𝐵

₁

𝑏

₁₂

= 𝐴

₂

𝑎

₂₂

+ 𝐵

₂

𝑏

₂₂

3. Along the negative x-axis (𝑦 = 0):

𝐴

₂

𝑎

₂₁

+ 𝐵

₂

𝑏

₂₁

= 𝐴

₃

𝑎

₃₁

+ 𝐵

₃

𝑏

₃₁

𝐴

₂

𝑎

₂₂

+ 𝐵

₂

𝑏

₂₂

= 𝐴

₃

𝑎

₃₂

+ 𝐵

₃

𝑏

₃₂

4. Along the negative y-axis (𝑥 = 0):

𝐴

₃

𝑎

₃₁

+ 𝐵

₃

𝑏

₃₁

= 𝐴

₄

𝑎

₄₁

+ 𝐵

₄

𝑏

₄₁

𝐴

₃

𝑎

₃₂

+ 𝐵

₃

𝑏

₃₂

= 𝐴

₄

𝑎

₄₂

+ 𝐵

₄

𝑏

₄₂

Equation 66

By using Cramer’s rule, one can acquire solutions for the individual unknowns. These can be expressed

as follows:

Finite Analytic Method Extension: Anisotropic Solution

45 | P a g e 𝐴

₁

= (𝐴

₄

𝑎

₄₁

+ 𝐵

₄

𝑏

₄₁

)𝑏

₁₂

− (𝐴

₄

𝑎

₄₂

+ 𝐵

₄

𝑏

₄₂

)𝑏

₁₁

𝑎

₁₁

𝑏

₁₂

− 𝑎

₁₂

𝑏

₁₁

𝐵

₁

= −(𝐴

₄

𝑎

₄₁

+ 𝐵

₄

𝑏

₄₁

)𝑎

₁₂

+ (𝐴

₄

𝑎

₄₂

+ 𝐵

₄

𝑏

₄₂

)𝑎

₁₁

𝑎

₁₁

𝑏

₁₂

− 𝑎

₁₂

𝑏

₁₁

𝐴

₂

= (𝐴

₁

𝑎

₁₁

+ 𝐵

₁

𝑏

₁₁

)𝑏

₂₂

− (𝐴

₁

𝑎

₁₂

+ 𝐵

₁

𝑏

₁₂

)𝑏

₂₁

𝑎

₂₁

𝑏

₂₂

− 𝑎

₂₂

𝑏

₂₁

𝐵

₂

= −(𝐴

₁

𝑎

₁₁

+ 𝐵

₁

𝑏

₁₁

)𝑎

₂₂

+ (𝐴

₁

𝑎

₁₂

+ 𝐵

₁

𝑏

₁₂

)𝑎

₂₁

𝑎

₂₁

𝑏

₂₂

− 𝑎

₂₂

𝑏

₂₁

𝐴

₃

= (𝐴

₂

𝑎

₂₁

+ 𝐵

₂

𝑏

₂₁

)𝑏

₃₂

− (𝐴

₂

𝑎

₂₂

+ 𝐵

₂

𝑏

₂₂

)𝑏

₃₁

𝑎

₃₁

𝑏

₃₂

− 𝑎

₃₂

𝑏

₃₁

𝐵

₃

= −(𝐴

2

𝑎

₂₁

+ 𝐵

₂

𝑏

₂₁

)𝑎

₃₂

+ (𝐴

2

𝑎

₂₂

+ 𝐵

₂

𝑏

₂₂

)𝑎

₃₁

𝑎

₃₁

𝑏

₃₂

− 𝑎

₃₂

𝑏

₃₁

𝐴

₄

= (𝐴

₃

𝑎

₃₁

+ 𝐵

₃

𝑏

₃₁

)𝑏

₄₂

− (𝐴

₃

𝑎

₃₂

+ 𝐵

₃

𝑏

₃₂

)𝑏

₄₁

𝑎

₄₁

𝑏

₄₂

− 𝑎

₄₂

𝑏

₄₁

𝐵

₄

= −(𝐴

3

𝑎

₃₁

+ 𝐵

₃

𝑏

₃₁

)𝑎

₄₂

+ (𝐴

3

𝑎

₃₂

+ 𝐵

₃

𝑏

₃₂

)𝑎

₄₁

𝑎

₄₁

𝑏

₄₂

− 𝑎

₄₂

𝑏

₄₁

Equation 67

This allows for the unknowns to be solved as long as 𝛼 is determined, as it is the only unknown in the equations above other than the subjects.

By changing notation, a more compact set of equations can be formed:

Finite Analytic Method Extension: Anisotropic Solution

46 | P a g e 𝑢

₁₁

= 𝑎

₄₁

𝑏

₁₂

− 𝑎

₄₂

𝑏

₁₁

𝑎

₁₁

𝑏

₁₂

− 𝑎

₁₂

𝑏

₁₁

, 𝑢

₁₂

= −𝑎

₄₁

𝑎

₁₂

+ 𝑎

₄₂

𝑎

₁₁

𝑎

₁₁

𝑏

₁₂

− 𝑎

₁₂

𝑏

₁₁

𝑣

₁₁

= 𝑏

₄₁

𝑏

₁₂

− 𝑏

₄₂

𝑏

₁₁

𝑎

₁₁

𝑏

₁₂

− 𝑎

₁₂

𝑏

₁₁

, 𝑣

₁₂

= −𝑏

₄₁

𝑎

₁₂

+ 𝑏

₄₂

𝑎

₁₁

𝑎

₁₁

𝑏

₁₂

− 𝑎

₁₂

𝑏

₁₁

𝑢

₂₁

= 𝑎

₁₁

𝑏

₂₂

− 𝑎

₁₂

𝑏

₂₁

𝑎

₂₁

𝑏

₂₂

− 𝑎

₂₂

𝑏

₂₁

, 𝑢

₂₂

= −𝑎

₁₁

𝑎

₂₂

+ 𝑎

₁₂

𝑎

₂₁

𝑎

₂₁

𝑏

₂₂

− 𝑎

₂₂

𝑏

₂₁

𝑣

₂₁

= 𝑏

₁₁

𝑏

₂₂

− 𝑏

₁₂

𝑏

₂₁

𝑎

₂₁

𝑏

₂₂

− 𝑎

₂₂

𝑏

₂₁

, 𝑣

₂₂

= −𝑏

₁₁

𝑎

₂₂

+ 𝑏

₁₂

𝑎

₂₁

𝑎

₂₁

𝑏

₂₂

− 𝑎

₂₂

𝑏

₂₁

𝑢

₃₁

= 𝑎

₂₁

𝑏

₃₂

− 𝑎

₂₂

𝑏

₃₁

𝑎

₃₁

𝑏

₃₂

− 𝑎

₃₂

𝑏

₃₁

, 𝑢

₃₂

= −𝑎

₂₁

𝑎

₃₂

+ 𝑎

₂₂

𝑎

₃₁

𝑎

₃₁

𝑏

₃₂

− 𝑎

₃₂

𝑏

₃₁

𝑣

₃₁

= 𝑏

₂₁

𝑏

₃₂

− 𝑏

₂₂

𝑏

₃₁

𝑎

₃₁

𝑏

₃₂

− 𝑎

₃₂

𝑏

₃₁

, 𝑣

₃₂

= −𝑏

₂₁

𝑎

₃₂

+ 𝑏

₂₂

𝑎

₃₁

𝑎

₃₁

𝑏

₃₂

− 𝑎

₃₂

𝑏

₃₁

𝑢

₄₁

= 𝑎

₃₁

𝑏

₄₂

− 𝑎

₃₂

𝑏

₄₁

𝑎

₄₁

𝑏

₄₂

− 𝑎

₄₂

𝑏

₄₁

, 𝑢

₄₂

= −𝑎

₃₁

𝑎

₄₂

+ 𝑎

₃₂

𝑎

₄₁

𝑎

₄₁

𝑏

₄₂

− 𝑎

₄₂

𝑏

₄₁

𝑣

₄₁

= 𝑏

₃₁

𝑏

₄₂

− 𝑏

₃₂

𝑏

₄₁

𝑎

₄₁

𝑏

₄₂

− 𝑎

₄₂

𝑏

₄₁

, 𝑣

₄₂

= −𝑏

₃₁

𝑎

₄₂

+ 𝑏

₃₂

𝑎

₄₁

𝑎

₄₁

𝑏

₄₂

− 𝑎

₄₂

𝑏

₄₁

Equation 68

Thus, the equations for the unknowns 𝐴

_𝑖

and 𝐵

_𝑖

are as follows:

𝐴

₁

= 𝐴

₄

𝑢

₁₁

+ 𝐵

₄

𝑣

₁₁

𝐵

₁

= 𝐴

₄

𝑢

₁₂

+ 𝐵

₄

𝑣

₁₂

𝐴

₂

= 𝐴

₁

𝑢

₂₁

+ 𝐵

₁

𝑣

₂₁

𝐵

₂

= 𝐴

₁

𝑢

₂₂

+ 𝐵

₁

𝑣

₂₂

𝐴

₃

= 𝐴

₂

𝑢

₃₁

+ 𝐵

₂

𝑣

₃₁

𝐵

₃

= 𝐴

₂

𝑢

₃₂

+ 𝐵

₂

𝑣

₃₂

𝐴

₄

= 𝐴

₃

𝑢

₄₁

+ 𝐵

₃

𝑣

₄₁

𝐵

₄

= 𝐴

₃

𝑢

₄₂

+ 𝐵

₃

𝑣

₄₂

Equation 69

Applying the backward substitution method here, similar to that seen in the

Finite Analytic Method section, a general solution for 𝛼 can be derived. Starting with the expressions

deduced from the substitution:

Finite Analytic Method Extension: Anisotropic Solution

47 | P a g e 𝐴

₁

= 𝐴

₃

[𝑢

₄₁

𝑢

₁₁

+ 𝑢

₄₂

𝑣

₁₁

] + 𝐵

₃

[𝑣

₄₁

𝑢

₁₁

+ 𝑣

₄₂

𝑣

₁₁

]

= 𝐴

₂

[𝑢

₃₁

𝑢

₄₁

𝑢

₁₁

+ 𝑢

₃₁

𝑢

₄₂

𝑣

₁₁

+ 𝑢

₃₂

𝑣

₄₁

𝑢

₁₁

+ 𝑢

₃₂

𝑣

₄₂

𝑣

₁₁

]

+ 𝐵

₂

[𝑣

₃₁

𝑢

₄₁

𝑢

₁₁

+ 𝑣

₃₁

𝑢

₄₂

𝑣

₁₁

+ 𝑣

₃₂

𝑣

₄₁

𝑢

₁₁

+ 𝑣

₃₂

𝑣

₄₂

𝑣

₁₁

]

= 𝐴

₁

[𝑢

₂₁

𝑢

₃₁

𝑢

₄₁

𝑢

₁₁

+ 𝑢

₂₁

𝑢

₃₁

𝑢

₄₂

𝑣

₁₁

+ 𝑢

₂₁

𝑢

₃₂

𝑣

₄₁

𝑢

₁₁

+ 𝑢

₂₁

𝑢

₃₂

𝑣

₄₂

𝑣

₁₁

+ 𝑢

₂₂

𝑣

₃₁

𝑢

₄₁

𝑢

₁₁

+ 𝑢

₂₂

𝑣

₃₁

𝑢

₄₂

𝑣

₁₁

+ 𝑢

₂₂

𝑣

₃₂

𝑣

₄₁

𝑢

₁₁

+ 𝑢

₂₂

𝑣

₃₂

𝑣

₄₂

𝑣

₁₁

]

+ 𝐵

₁

[𝑣

₂₁

𝑢

₃₁

𝑢

₄₁

𝑢

₁₁

+ 𝑣

₂₁

𝑢

₃₁

𝑢

₄₂

𝑣

₁₁

+ 𝑣

₂₁

𝑢

₃₂

𝑣

₄₁

𝑢

₁₁

+ 𝑣

₂₁

𝑢

₃₂

𝑣

₄₂

𝑣

₁₁

+ 𝑣

₂₂

𝑣

₃₁

𝑢

₄₁

𝑢

₁₁

+ 𝑣

₂₂

𝑣

₃₁

𝑢

₄₂

𝑣

₁₁

+ 𝑣

₂₂

𝑣

₃₂

𝑣

₄₁

𝑢

₁₁

+ 𝑣

₂₂

𝑣

₃₂

𝑣

₄₂

𝑣

₁₁

]

= 𝐴

₁

𝜑

₁₁

+ 𝐵

₁

𝜑

₂₁

𝐵

₁

= 𝐴

₃

[𝑢

₄₁

𝑢

₁₂

+ 𝑢

₄₂

𝑣

₁₂

] + 𝐵

₃

[𝑣

₄₁

𝑢

₁₂

+ 𝑣

₄₂

𝑣

₁₂

]

= 𝐴

₂

[𝑢

₃₁

𝑢

₄₁

𝑢

₁₂

+ 𝑢

₃₁

𝑢

₄₂

𝑣

₁₂

+ 𝑢

₃₂

𝑣

₄₁

𝑢

₁₂

+ 𝑢

₃₂

𝑣

₄₂

𝑣

₁₂

]

+ 𝐵

₂

[𝑣

₃₁

𝑢

₄₁

𝑢

₁₂

+ 𝑣

₃₁

𝑢

₄₂

𝑣

₁₂

+ 𝑣

₃₂

𝑣

₄₁

𝑢

₁₂

+ 𝑣

₃₂

𝑣

₄₂

𝑣

₁₂

]

= 𝐴

₁

[𝑢

₂₁

𝑢

₃₁

𝑢

₄₁

𝑢

₁₂

+ 𝑢

₂₁

𝑢

₃₁

𝑢

₄₂

𝑣

₁₂

+ 𝑢

₂₁

𝑢

₃₂

𝑣

₄₁

𝑢

₁₂

+ 𝑢

₂₁

𝑢

₃₂

𝑣

₄₂

𝑣

₁₂

+ 𝑢

₂₂

𝑣

₃₁

𝑢

₄₁

𝑢

₁₂

+ 𝑢

₂₂

𝑣

₃₁

𝑢

₄₂

𝑣

₁₂

+ 𝑢

₂₂

𝑣

₃₂

𝑣

₄₁

𝑢

₁₂

+ 𝑢

₂₂

𝑣

₃₂

𝑣

₄₂

𝑣

₁₂

]

+ 𝐵

₁

[𝑣

₂₁

𝑢

₃₁

𝑢

₄₁

𝑢

₁₂

+ 𝑣

₂₁

𝑢

₃₁

𝑢

₄₂

𝑣

₁₂

+ 𝑣

₂₁

𝑢

₃₂

𝑣

₄₁

𝑢

₁₂

+ 𝑣

₂₁

𝑢

₃₂

𝑣

₄₂

𝑣

₁₂

+ 𝑣

₂₂

𝑣

₃₁

𝑢

₄₁

𝑢

₁₂

+ 𝑣

₂₂

𝑣

₃₁

𝑢

₄₂

𝑣

₁₂

+ 𝑣

₂₂

𝑣

₃₂

𝑣

₄₁

𝑢

₁₂

+ 𝑣

₂₂

𝑣

₃₂

𝑣

₄₂

𝑣

₁₂

]

= 𝐴

₁

𝜑

₁₂

+ 𝐵

₁

𝜑

₂₂

Equation 70

The variables 𝜑

_𝑖𝑗

are simply a change of notation, used to make the equations shorter and more comprehensible.

Finally, the following can be obtained:

𝐴

₁

(𝜑

₁₁

− 1) + 𝐵

₁

𝜑

₂₁

= 0 𝐴

₁

𝜑

₁₂

+ 𝐵

₁

(𝜑

₂₂

− 1) = 0

Equation 71

Again, the corresponding determinant is taken is set to equal zero, as follows:

| (𝜑

₁₁

− 1) 𝜑

₂₁

𝜑

₁₂

(𝜑

₂₂

− 1) | = 𝜑

₁₁

𝜑

₂₂

− 𝜑

₁₂

𝜑

₂₁

− 𝜑

₁₁

− 𝜑

₂₂

+ 1 = 0

Equation 72

Note: 𝜑

_𝑖𝑗

are the components derived from the substitution method.

Now, a solution for Equation 72 must be found in order to determine the value of 𝛼. The Newton Secant

method would be applied to find the best fit for 𝛼. This is a numerical approach to the solution of 𝛼, as

opposed to the analytical solution that was determined in the previous section for the isotropic

Finite Analytic Method Extension: Anisotropic Solution

48 | P a g e permeability. Using the equation established above, the same solution for 𝛼 was achieved for the checkerboard test case used in this study, despite the two methods and algorithms varying so greatly in the approach. In addition to finding a solution for 𝛼, the constant 𝐶 =

^𝐵1

𝐴1

can be solved by taking either 𝐶 = −

^(𝜑_𝜑¹¹⁻¹⁾

21

, or 𝐶 = −

^𝜑12

(𝜑₂₂−1)

. This was also determined to give the same value as the theoretical approach described in the

Finite Analytic Method section.

Newton Secant, also referred to as the Newton-Raphson method, is a root seeking algorithm where better approximations for the roots are found in a successive fashion (Papakonstantinou, 2007). The approach starts with a function 𝑓 and its corresponding derivative 𝑓

^′

(approximated numerically in the approach used in this study, see Equation 73), together with an initial guess 𝑥

₀

for the root of the function. If the guess converges towards a better assumption, then the process is repeated until a sufficiently better approximation is found based on a criteria and/or margin of error. The general formula takes the form of:

𝑓

^′

(𝑥

₁

) = 𝑓(𝑥

₁

) − 𝑓(𝑥

₀

)

𝑥

₁

− 𝑥

₀

→ 𝑓

^′

(𝑥

_𝑛+1

) = 𝑓(𝑥

_𝑛+1

) − 𝑓(𝑥

_𝑛

) 𝑥

_𝑛+1

− 𝑥

_𝑛

𝑥

₁

= 𝑥

₀

− 𝑓(𝑥

₀

)

𝑓

^′

(𝑥

₀

) → 𝑥

_𝑛+1

= 𝑥

_𝑛

− 𝑓(𝑥

_𝑛

) 𝑓

^′

(𝑥

_𝑛

)

Equation 73

In this manner, the finite analytic method is reconfigured to support anisotropic permeability. This novel method differs greatly from what was recently developed by the same authors of the article that was explored in the

Finite Analytic Method section. The approach that they tackled was a more generalized form of the

Laplace equations (Liu & Wang, 2016), where an entirely new set of equations are derived. However, the

novel method developed during this study is entirely independent.

Results

49 | P a g e

In document Two-Dimensional Fluid Flow in Heterogeneous Porous Media Using Finite Analytic Method (sider 38-59)