Comment déterminez-vous la résistance effective d'une grille finie de résistances?

Avertissement: Je suis géophysicien avec une formation limitée en génie électrique. Je ne sais pas si ce problème est incroyablement facile, incroyablement complexe ou complètement insensé.

Mon objectif: déterminer la résistivité globale d'un échantillon de roche à l'aide de réseaux de résistances.

L'échantillon de roche doit être modélisé en utilisant un réseau de résistances avec certaines résistances ayant une résistance élevée (représentant la roche solide) et d'autres résistances ayant une faible résistance (représentant les voies de fluide dans la roche).

Supposons que j'ai un réseau de résistances sur une grille uniforme comme indiqué ci-dessous. Dans l'exemple illustré, chaque segment de ligne a une résistance associée étiquetée de 1 à 24 sur une grille 3 par 3. Les résistances de chaque segment de ligne sont connues.

La longueur totale de la grille est et la "zone" est A (dans ce cas, il s'agit d'un exemple 2D, donc la zone n'est également qu'une longueur). La résistivité globale de l'échantillon est alors donnée par:$L$$A$

$$\rho_{bulk} = \frac{R_{eff}A}{L}$$

Ma question: comment déterminer la résistance effective, du réseau?$R_{eff}$

J'ai regardé en ligne mais tout ce que je peux trouver, ce sont des discussions sur les réseaux infinis , les courants de source et de puits, etc. Je ne suis pas intéressé par le courant ou la tension.

Ce problème peut-il être résolu tel quel?

L'idée de base est assez simple. Vous disposez d'une matrice ( ) qui représente des "nœuds" ou sommets dans votre système. Chacun de ces nœuds est associé à une "tension" à valeur scalaire qui peut être modifiée ou mise à jour au fur et à mesure de l'algorithme. Il y aura également deux nœuds dont la tension ne peut pas être modifiée. Nous allons appliquer ici une sorte de "batterie", donc ces deux nœuds représentent les deux extrémités de cette batterie.$V$

Séparément, deux autres matrices ( et R h ) représentent les bords du système, horizontal et vertical. Ce sont vos valeurs de résistance, je suppose. Je ne sais pas comment vous comptez les remplir. Mais c'est votre problème. Cette technique suppose que vous pouvez également remplir ces matrices.$Rv$$Rh$

Selon le langage informatique que vous utilisez, vous pouvez ou non être en mesure d'utiliser des indices négatifs. Peu importe. C'est juste une question de garder à l'esprit ce à quoi vous êtes confronté.

Supposons que la longueur soit divisée en N sections L et que la "longueur" A soit divisée en N sections A. Ensuite, vous devrez construire une matrice avec ( N L + 1 ) ⋅ ( N A + 1 ) sommets pour les valeurs de tension scalaire. (ou plus.) Vous aurez également besoin de ces deux autres matrices avec N A ⋅ ( N L + 1 ) bords verticaux et N L ⋅ ( N A + $L$$N_L$$A$$N_A$$\left(N_L+1\right)\cdot\left(N_A+1\right)$$N_A\cdot\left(N_L+1\right)$ bords horizontaux entre ces sommets.$N_L\cdot\left(N_A+1\right)$

Maintenant. Initialise tous les sommets avec . Choisissez l'un des sommets à gauche (au milieu, de préférence) et notez-le comme un$0\:\textrm{V}$Valeur V qui n'est PAS autorisée à changer. Utilisez la méthode que vous souhaitez pour cela. Choisissez l'un des sommets à droite (au milieu, de préférence) et changez sa valeur en$0\:\textrm{V}$ , tout en prenant note à nouveau que sa valeur n'est PAS autorisée à changer. Une technique qui fonctionne ici consiste simplement à le laisser changer normalement, puis à remplacer la valeur à chaque étape. Mais peu importe comment vous y parvenez, tant que vous y parvenez.$1\:\textrm{V}$

(Il existe d'autres techniques pour des raisons d'efficacité. Mais cela ne vaut probablement pas la peine de les utiliser ici.)

Passons maintenant à l'algorithme, qui est parfois appelé algorithme en damier ou rouge-noir . En vous déplaçant dans la matrice de tension de votre nœud, traitez chaque nœud où la somme des deux indices, est paire, en effectuant l'affectation simple suivante:$i+j$

$$\begin{align*} V_{i,j}&=\frac{Rh_{i,j-1}\cdot Rh_{i,j}\cdot\left(V_{i-1,j}\cdot Rv_{i,j}+V_{i+1,j}\cdot Rv_{i-1,j}\right)}{Rh_{i,j-1}\cdot Rh_{i,j}\cdot \left(Rv_{i,j}+Rv_{i-1,j}\right)+Rv_{i-1,j}\cdot Rv_{i,j}\left(Rh_{i,j}+Rh_{i,j-1}\right)}\\\\ &+\frac{Rv_{i-1,j}\cdot Rv_{i,j}\cdot \left(V_{i,j-1}\cdot Rh_{i,j}+V_{i,j+1}\cdot Rh_{i,j-1}\right)}{Rh_{i,j-1}\cdot Rh_{i,j}\cdot \left(Rv_{i,j}+Rv_{i-1,j}\right)+Rv_{i-1,j}\cdot Rv_{i,j}\left(Rh_{i,j}+Rh_{i,j-1}\right)} \end{align*}$$

L'équation ci-dessus n'est rien de plus que le calcul de la tension d'un nœud central ayant quatre résistances qui s'y connectent, où les tensions aux autres extrémités des quatre résistances sont connues. La tension du nœud central est ensuite calculée à partir de l'équation ci-dessus. Étant donné que le diviseur est le même pour chaque terme, vous pouvez simplement calculer la somme des numérateurs, puis diviser une fois par le dénominateur.

Cela mettra à jour tous les nœuds où la somme est paire. Maintenant, vous effectuez la même procédure pour tous les nœuds où la somme i + $i+j$$i+j$ est impaire. Une fois ces deux étapes effectuées, vous avez terminé un cycle.

Si nécessaire, réinitialisez les deux nœuds spéciaux (pour $0\:\textrm{V}$ et pour $1\:\textrm{V}$ comme expliqué précédemment.) Ou, si vous avez protégé ces deux nœuds, il n'est pas nécessaire de les réinitialiser.

Vous êtes prêt pour le prochain cycle. Effectuez ces cycles autant de fois que vous le jugez nécessaire pour que l'état général se stabilise (et ce sera le cas.)

Lorsque vous arrêtez le processus, vous pouvez facilement déterminer la résistance en choisissant de regarder les nœuds entourant votre nœud protégé gauche ou bien de regarder les nœuds entourant votre nœud protégé droit. (Il peut être judicieux d'agrandir votre matrice juste assez [par 1 dans toutes les directions] pour avoir en fait quatre nœuds entourant l'un ou l'autre choix.) La différence de tension entre les nœuds environnants et le nœud spécial, divisée par le la résistance dans les bords entre eux vous indique la sortie / entrée actuelle de votre nœud spécial. Puisqu'il s'agit d'un nœud "batterie", ce courant doit être TOUT du courant. Puisque la tension est$1\:\textrm{V}$ , par définition, divisant 1 par la somme de ces quatre courants que vous trouvez vous indique la résistance totale.

Je regarde un code que j'ai écrit qui totalise, avec beaucoup de commentaires, seulement 67 lignes. Il n'est donc PAS difficile d'écrire.

$1\:\textrm{V}$

---

Pourquoi faut-il séparer le système en i + j = pair et i + j = impair?

$V_{5,5}=f\left(V_{4,5},V_{6,5},V_{5,4},V_{5,6}\right)$$V_{5,5}$$V_{5,6}=f\left(V_{4,6},V_{6,6},V_{5,5},V_{5,7}\right)$$V_{5,5}$$V_{5,7}=f\left(V_{4,7},V_{6,7},V_{5,6},V_{5,8}\right)$car aucune des entrées de la fonction n'est un nœud qui a été modifié au cours de cette étape. Ensuite, vous vous retournez et calculez les suppléants, en évitant le maculage, mais maintenant en mettant à jour les suppléants. Vous devez vraiment le faire de cette façon.

De plus, la formule est-elle identique pour les étapes paires et impaires?

Oui, c'est pareil.

Tout cela peut-il être résolu en une seule étape en utilisant une sorte de système linéaire Ax = b où A est un opérateur linéaire et b fournit les conditions aux limites? En le regardant, cela semble quelque peu analogue aux méthodes de différences finies pour résoudre des équations aux dérivées partielles.

Il y a une connexion. Je pense que cela s'appelle une implémentation «sans matrice».

---

Voici un exemple. L'ensemble de valeurs de résistance suivant a été placé dans LTSpice pour la simulation:

$1\:\textrm{V}$$30.225\:\textrm{mA}$$30.224552\:\textrm{mA}$

J'ai exécuté le programme VB.NET suivant:

Module GEOGRID

    Const NL As Integer = 2
    Const NA As Integer = 2
    Const INF As Double = 1.0E+32

    Sub Main()

        Static Rh As Double(,) = New Double(NL + 2, NA + 1) {
                    {INF, INF, INF, INF},
                    {INF, 5, 21, INF},
                    {INF, 76, 10, INF},
                    {INF, 32, 22, INF},
                    {INF, INF, INF, INF}}
        Static Rv As Double(,) = New Double(NA + 1, NL + 2) {
                    {INF, INF, INF, INF, INF},
                    {INF, 61, 50, 16, INF},
                    {INF, 56, 45, 18, INF},
                    {INF, INF, INF, INF, INF}}
        Dim V As Double(,) = New Double(NL + 2, NA + 2) {
                    {0, 0, 0, 0, 0},
                    {0, 0, 0, 0, 0},
                    {0, 0, 0, 1, 0},
                    {0, 0, 0, 0, 0},
                    {0, 0, 0, 0, 0}}
        Dim PDE As Func(Of Integer, Integer, Double) = Function(ByVal i As Integer, ByVal j As Integer) (
                    Rh(i, j - 1) * Rh(i, j) * (V(i - 1, j) * Rv(i, j) + V(i + 1, j) * Rv(i - 1, j)) +
                    Rv(i - 1, j) * Rv(i, j) * (V(i, j - 1) * Rh(i, j) + V(i, j + 1) * Rh(i, j - 1))
                  ) / (
                    Rh(i, j - 1) * Rh(i, j) * (Rv(i, j) + Rv(i - 1, j)) +
                    Rv(i - 1, j) * Rv(i, j) * (Rh(i, j) + Rh(i, j - 1))
                  )
        Dim IV As Func(Of Integer, Integer, Double) = Function(ByVal i As Integer, ByVal j As Integer) 0 +
                    (V(i, j) - V(i - 1, j)) / Rv(i - 1, j) + (V(i, j) - V(i + 1, j)) / Rv(i, j) +
                    (V(i, j) - V(i, j - 1)) / Rh(i, j - 1) + (V(i, j) - V(i, j + 1)) / Rh(i, j)
        Dim idx As Integer = NA \ 2 + 1
        Dim jdx1 As Integer = NL + 1
        Dim jdx2 As Integer = 1
        For x As Integer = 1 To 1000
            For k As Integer = 0 To (NA + 1) * (NL + 1) - 1 Step 2
                Dim i As Integer = k \ (NL + 1)
                Dim j As Integer = k - i * (NL + 1) + 1
                i += 1
                If Not (i = idx AndAlso (j = jdx1 OrElse j = jdx2)) Then V(i, j) = PDE(i, j)
            Next
            For k As Integer = 1 To (NA + 1) * (NL + 1) - 1 Step 2
                Dim i As Integer = k \ (NL + 1)
                Dim j As Integer = k - i * (NL + 1) + 1
                i += 1
                If Not (i = idx AndAlso (j = jdx1 OrElse j = jdx2)) Then V(i, j) = PDE(i, j)
            Next
        Next
        Console.WriteLine("R = " & (1.0 / IV(idx, jdx1)).ToString)
        Console.WriteLine("R = " & (-1.0 / IV(idx, jdx2)).ToString)
    End Sub

End Module

$R = 33.0856844038614\:\Omega$

Le programme ci-dessus montre un moyen de configurer les résistances, verticales et horizontales, ainsi que la matrice de tension, de sorte qu'il simplifie certains des tests pour les nœuds et / ou les valeurs de résistance inexistants. Le code est un peu plus propre de cette façon, bien qu'il nécessite quelques éléments de tableau supplémentaires. (J'ai simplement rendu les valeurs de résistance supplémentaires infinies.) Il suffit de comparer la façon dont j'ai configuré les tableaux avec la façon dont le schéma a été présenté également, et je pense que vous serez en mesure de déterminer exactement détails ici.

J'ai également piraté les résistances et les valeurs des nœuds, bien sûr, sans en faire un programme à usage général pour lire un tableau de valeurs. Mais cette généralité est assez facile à ajouter. Et ce code devrait rendre tout ce que j'ai écrit absolument sans ambiguïté.

$x$$x$

$0\:\textrm{V}$$1\:\textrm{V}$

(D'accord. Une dernière note. Ce serait beaucoup mieux ciblé sur F # ou tout compilateur décent ciblant un système informatique massivement parallèle. Chaque calcul en "rouge" ou "noir" peut être effectué en parallèle, complètement indépendamment les uns des autres. F # rend cela trivial. Donc codé en F #, vous pouvez l'exécuter sur tous vos cœurs disponibles sans rien de spécial à faire. Cela fonctionne juste. Juste une note au cas où vous collectez BEAUCOUP de données d'une manière et que vous pourriez vouloir prendre tous les avantages d'un système multicœur.)

---

NOTE DE FIN:

La dérivation est assez simple à partir de KCL. Placez quatre résistances dans la disposition suivante:

^{simuler ce circuit - Schéma créé à l'aide de CircuitLab}

Appliquer KCL:

$$\begin{align*} \frac{V}{R_1}+\frac{V}{R_2}+\frac{V}{R_3}+\frac{V}{R_4} &= \frac{V_1}{R_1}+\frac{V_2}{R_2}+\frac{V_3}{R_3}+\frac{V_4}{R_4}\\\\ &\therefore\\\\ V &=\left(\frac{V_1}{R_1}+\frac{V_2}{R_2}+\frac{V_3}{R_3}+\frac{V_4}{R_4}\right)\bigg(R_1 \mid\mid R_2 \mid\mid R_3 \mid\mid R_4\bigg) \end{align*}$$

Certains jouer avec l'algèbre obtient le résultat que j'ai utilisé dans le code.

Vous pouvez certainement adopter l'approche d'un réseau de résistances 2D pour modéliser un problème 2D mais cela peut devenir quelque peu délicat lors du passage à 3 dimensions. Vous voudrez peut-être envisager d'utiliser une approche plus traditionnelle (de nos jours) avec des conducteurs de volume définis dans vos domaines avec une conductivité appropriée attribuée à chacun. Le code freeware FEMM ( http://www.femm.info/wiki/HomePage ) est très performant et peut être utilisé pour la 2D, la symétrie axiale et la 3D. De là, vous pouvez envisager de passer à des codes beaucoup plus performants comme SCIrun ( https://www.sci.utah.edu/), qui est un code académique pour les problèmes de conducteurs de volume d'une complexité substantielle. Je l'utilise régulièrement pour des mailles de plus d'un million de tétraèdres. Même s'il a été principalement développé pour la modélisation biologique, il devrait fonctionner parfaitement pour ce que vous faites. Les exemples de problèmes avancés dans la boîte à outils avant / inverse devraient vous aider. Les problèmes inverses peuvent également être utiles pour la tomographie par impédance. J'utilise généralement la version 4 car la version 5 est toujours en cours de réalisation. Le logiciel possède également une interface avec tetgen qui est un excellent code de construction de maillage.

Enfin, si vous n'êtes pas opposé à dépenser de l'argent, il y a toujours COMSOL, qui est très facile à utiliser (et assez cher).